Commits

J.S. Oishi  committed edc970b

more changes to get Alfven waves working in parallel

  • Participants
  • Parent commits 3941576

Comments (0)

Files changed (3)

File dedalus/data_objects/fourier_data.py

 import numpy.fft as fpack
 from dedalus.utils.parallelism import comm, MPI
 from dedalus.utils.fftw import fftw
-
+from dedalus.funcs import insert_ipython
 try:
     import pycuda.autoinit
     import pycuda.gpuarray as gpuarray
             ki = fpack.fftfreq(S) * 2. * self.kny[i]
             ki.resize(kshape)
             self.k.append(ki)
+        insert_ipython()
         self.k = dict(zip(['z','y','x'][3-self.ndim:], self.k))
 
     def set_fft(self, method):
         self.data.shape = self._shape['kspace']
         self.data[:] = a
         self.fplan_x()
-        self.data /= na.sqrt(self.data.size)
+        self.data /= self.data.size
 
     def rev_fftw(self):
         self.rplan_x()
         self.data.shape = self._shape['xspace']
         self.data[:] = a
         self.rplan_yz()
-        self.data /= na.sqrt(self.data.size)
         self.data.imag = 0.
 
 class ParallelFourierShearRepresentation(ParallelFourierRepresentation, FourierShearRepresentation):

File dedalus/init_cond/init_cond.py

     f[tuple(-1*na.array(kindex))] = f[tuple(kindex)].conjugate()
 
 
-def alfven(data, k=(0, 0, 1), B0mag=5.0, u1mag=5e-6):
+def alfven(data, k=(1, 0, 0), B0mag=5.0, u1mag=5e-6):
     """
     Generate conditions for simulating Alfven waves in MHD.
     For 2d, must have k and B0 in same direction.
     kmag = na.linalg.norm(k)
     
     # Field setup and calculation
-    B0 = na.array([0., 0., 1.])[:N] * B0mag
+    B0 = na.array([1., 0., 0.])[:N] * B0mag
     
     # Background magnetic field
     if (0 in data['u']['x'].k['z'] and
             data['B'][i]['kspace'][k0] = B0[i]
     
     # Assign modes for k and -k
-    if (k[0] in data['u']['x'].k['z'] and
+    if (k[2] in data['u']['x'].k['z'] and
         k[1] in data['u']['x'].k['y'] and
-        k[2] in data['u']['x'].k['x']):
+        k[0] in data['u']['x'].k['x']):
         
-        kiz = na.array(na.where(data['u']['x'].k['z'] == k[0]))
+        kiz = na.array(na.where(data['u']['x'].k['z'] == k[2]))
         kiy = na.array(na.where(data['u']['x'].k['y'] == k[1]))
-        kix = na.array(na.where(data['u']['x'].k['x'] == k[2]))
+        kix = na.array(na.where(data['u']['x'].k['x'] == k[0]))
         
         kindex = tuple(kiz + kiy + kix)
         # Alfven speed and wave frequency
             data['u'][i]['kspace']
             data['B'][i]['kspace']
         
-            data['u'][i].data[kindex] = -u1[i] * 1j / 2.
+            data['u'][i].data[kindex] = u1[i] * 1j / 2.
             data['B'][i].data[kindex] = B1[i] * 1j / 2.
             
-    if (-k[0] in data['u']['x'].k['z'] and
+    if (-k[2] in data['u']['x'].k['z'] and
          -k[1] in data['u']['x'].k['y'] and
-         -k[2] in data['u']['x'].k['x']):
+         -k[0] in data['u']['x'].k['x']):
           
         # Find correct mode
         # dkx = na.abs(data['u']['x'].k['x'][1] - data['u']['x'].k['x'][0])
 #         dkz = na.abs(data['u']['x'].k['z'][1] - data['u']['x'].k['z'][0])
 #         
 #         
-        kiz = na.array(na.where(data['u']['x'].k['y'] == -k[0]))
+        kiz = na.array(na.where(data['u']['x'].k['z'] == -k[2]))
         kiy = na.array(na.where(data['u']['x'].k['y'] == -k[1]))
-        kix = na.array(na.where(data['u']['x'].k['x'] == -k[2]))
+        kix = na.array(na.where(data['u']['x'].k['x'] == -k[0]))
         
         kindex = tuple(kiz + kiy + kix)
         # Alfven speed and wave frequency
         omega = na.abs(cA * na.dot(k, B0) / B0mag)
     
         # u and B perturbations
-        u1 = na.array([0., 1., 0.])[3-N:] * u1mag
+        u1 = na.array([0, 1., 0.])[3-N:] * u1mag
         
         B1 = (na.dot(k, u1) * B0 - na.dot(k, B0) * u1) / omega
         B1mag = na.linalg.norm(B1)
             data['u'][i]['kspace']
             data['B'][i]['kspace']
         
-            data['u'][i].data[kindex] = - u1[i] * 1j / 2.
-            data['B'][i].data[kindex] = - B1[i] * 1j / 2.
+            data['u'][i].data[kindex] = -u1[i] * 1j / 2.
+            data['B'][i].data[kindex] = -B1[i] * 1j / 2.
 
 
+def alfven_old(data, k=(1, 0, 0), B0mag=5.0, u1mag=5e-6):
+    """
+    Generate conditions for simulating Alfven waves in MHD.
+    For 2d, must have k and B0 in same direction.
+    
+    Inputs:
+        data        StateData object
+        k           k-index tuple, (kx, ky, kz), adds sine wave here in u1 and B1
+        B0mag       B0 magnitude along x-direction
+        u1mag       u1 magnitude for k-wave being added
+    """
+    
+    N = len(k)
+    kmag = na.linalg.norm(k)
+    
+    # Field setup and calculation
+    B0 = na.array([1., 0., 0.])[:N] * B0mag
+        
+    # Alfven speed and wave frequency
+    cA = B0mag / na.sqrt(4 * na.pi * data.parameters['rho0'])
+    omega = na.abs(cA * na.dot(k, B0) / B0mag)
+    print '-' * 20
+    print 'cA = ', cA
+    print 'cA * cos(theta) = ', cA * na.dot(k, B0) / B0mag / kmag
+    print 'Max dt (CFL=1) = ', na.min(na.array(data['u']['x'].length) / 
+                                      na.array(data['u']['x'].shape)) / cA
+    print '-' * 20
+    
+    # Background magnetic field
+    for i in xrange(data['B'].ndim):
+        k0 = (0,) * N
+        data['B'][i]['kspace'][k0] = B0[i]
+
+    # u and B perturbations
+    u1 = na.array([0., 1., 0.])[3-N:] * u1mag
+    
+    B1 = (na.dot(k, u1) * B0 - na.dot(k, B0) * u1) / omega
+    B1mag = na.linalg.norm(B1)
+    
+    for i in xrange(data['u'].ndim):
+        data['u'][i]['kspace']
+        data['B'][i]['kspace']
+    
+        sin_k(data['u'][i].data, k[::-1], ampl=u1[i])
+        sin_k(data['B'][i].data, k[::-1], ampl=B1[i])
+
 def turb(ux, uy, spec, tot_en=0.5, **kwargs):
     """generate noise with a random phase and a spectrum given by
     the spec function.

File dedalus/utils/parallelism.py

 
     data_file = os.path.join(snap_dir, 'data.cpu%04i')
     fi = h5py.File(data_file % 0)
+    time = fi['/time'].value
     space = fi['/fields/u/0'].attrs['space']
     local_size = fi['/fields/u/0'].shape
     dtype = fi['/fields/u/0'].dtype
         concat_ax = 2
     else:
         concat_ax = 0
+
+    print "loaded %s at time = %f" % (field, time)
     return np.concatenate(data,axis=concat_ax), space