Working on docs

demo/ChannelFlow.py

@@ -62,7 +62,7 @@ def __init__(self,
                  modsave=1e8,
                  moderror=100,
                  checkpoint=1000,
-                 timestepper='PDEIRK3'):
+                 timestepper='IMEXRK3'):
         self.nu = nu
         self.dt = dt
         self.conv = conv
@@ -292,8 +292,9 @@ def prepare_step(self, rk):
     def assemble(self):
         for pde in self.pdes.values():
             pde.assemble()
-        for pde in self.pdes1d.values():
-            pde.assemble()
+        if comm.Get_rank() == 0:
+            for pde in self.pdes1d.values():
+                pde.assemble()
 
     def solve(self, t=0, tstep=0, end_time=1000):
         self.assemble()

demo/ChannelFlow2D.py

@@ -87,7 +87,7 @@ def __init__(self,
         self.TD = TensorProductSpace(comm, (self.D0, self.F1), collapse_fourier=False, modify_spaces_inplace=True) # Streamwise velocity
         self.TC = TensorProductSpace(comm, (self.C0, self.F1), collapse_fourier=False, modify_spaces_inplace=True) # No bc
         self.BD = VectorSpace([self.TB, self.TD])  # Velocity vector space
-        self.CD = VectorSpace(self.TD)                      # Convection vector space
+        self.CD = VectorSpace(self.TD)             # Convection vector space
         self.CC = VectorSpace([self.TD, self.TC])  # Curl vector space
 
         # Padded space for dealiasing
@@ -124,12 +124,12 @@ def __init__(self,
                                      data={'0': {'U': [self.u_]}})
 
         # set up equations
+        #v = TestFunction(self.TB.get_testspace(kind='PG'))
         v = TestFunction(self.TB)
-        h = TestFunction(self.TD)
 
         # Chebyshev matrices are not sparse, so need a tailored solver. Legendre has simply 5 nonzero diagonals and can use generic solvers.
         sol1 = chebyshev.la.Biharmonic if self.B0.family() == 'chebyshev' else la.SolverGeneric1ND
-        sol2 = chebyshev.la.Helmholtz if self.B0.family() == 'chebyshev' else la.SolverGeneric1ND
+        #sol1 = la.SolverGeneric1ND
 
         self.pdes = {
 
@@ -238,8 +238,9 @@ def prepare_step(self, rk):
     def assemble(self):
         for pde in self.pdes.values():
             pde.assemble()
-        for pde in self.pdes1d.values():
-            pde.assemble()
+        if comm.Get_rank() == 0:
+            for pde in self.pdes1d.values():
+                pde.assemble()
 
     def solve(self, t=0, tstep=0, end_time=1000):
         self.assemble()
@@ -250,7 +251,7 @@ def solve(self, t=0, tstep=0, end_time=1000):
                     eq.compute_rhs(rk)
                 for eq in self.pdes.values():
                     eq.solve_step(rk)
-                self.compute_vw(rk)
+                self.compute_v(rk)
             t += self.dt
             tstep += 1
             self.update(t, tstep)

demo/MKM_MicroPolar.py

@@ -23,7 +23,7 @@ def __init__(self,
                  family='C',
                  padding_factor=(1, 1.5, 1.5),
                  checkpoint=1000,
-                 timestepper='PDEIRK3',
+                 timestepper='IMEXRK3',
                  rand=1e-7):
         MicroPolar.__init__(self, N=N, domain=domain, Re=Re, J=J, m=m, NP=NP, dt=dt, conv=conv, modplot=modplot,
                             modsave=modsave, moderror=moderror, filename=filename, family=family,
@@ -282,7 +282,7 @@ def fromfile(self, filename="stats"):
         'checkpoint': 1,
         'sample_stats': 10,
         'padding_factor': (1.5, 1.5, 1.5),
-        'timestepper': 'PDEIRK3', # IMEXRK222, IMEXRK443
+        'timestepper': 'IMEXRK222', # IMEXRK222, IMEXRK443
         }
     c = MKM(**d)
     t, tstep = c.initialize(from_checkpoint=False)

demo/OrrSommerfeld.py

@@ -29,11 +29,6 @@ def initialize(self, from_checkpoint=False):
         # This is the convection at t=0
         self.e0 = 0.5*dx(ub[0]**2+(ub[1]-(1-self.X[0]**2))**2)
         self.acc = np.zeros(1)
-        self.convection()
-        for pde in self.pdes.values():
-            pde.initialize()
-        for pde in self.pdes1d.values():
-            pde.initialize()
         return 0, 0
 
     def initOS(self, OS, eigvals, eigvectors, U, t=0.):
@@ -63,7 +58,7 @@ def compute_error(self, t):
     def init_plots(self):
         self.ub = ub = self.u_.backward()
         self.im1 = 1
-        if comm.Get_rank() == 0:
+        if comm.Get_rank() == 0 and comm.Get_size() == 1:
             plt.figure(1, figsize=(6, 3))
             self.im1 = plt.contourf(self.X[1][:, :, 0], self.X[0][:, :, 0], ub[0, :, :, 0], 100)
             plt.colorbar(self.im1)
@@ -84,7 +79,7 @@ def plot(self, t, tstep):
             self.init_plots()
 
         if tstep % self.modplot == 0 and self.modplot > 0:
-            if comm.Get_rank() == 0:
+            if comm.Get_rank() == 0 and comm.Get_size() == 1:
                 ub = self.u_.backward(self.ub)
                 X = self.X
                 self.im1.axes.clear()

demo/RayleighBenard2D.py

@@ -2,6 +2,7 @@
 from ChannelFlow2D import KMM
 import matplotlib.pyplot as plt
 import sympy
+
 np.warnings.filterwarnings('ignore')
 
 # pylint: disable=attribute-defined-outside-init
@@ -153,7 +154,7 @@ def plot(self, t, tstep):
 
         if tstep % self.modplot == 0 and self.modplot > 0:
             ub = self.u_.backward(self.ub)
-            Tb = self.T_.backward(self.Tb)
+            self.Tb = self.T_.backward(self.Tb)
             if comm.Get_rank() == 0:
                 plt.figure(1)
                 #self.im1.set_UVC(ub[1, ::4, ::4], ub[0, ::4, ::4])
@@ -162,7 +163,7 @@ def plot(self, t, tstep):
                 plt.pause(1e-6)
                 plt.figure(2)
                 self.im2.axes.clear()
-                self.im2.axes.contourf(self.X[1][:, :], self.X[0][:, :], Tb[:, :], 100)
+                self.im2.axes.contourf(self.X[1][:, :], self.X[0][:, :], self.Tb[:, :], 100)
                 self.im2.autoscale()
                 plt.pause(1e-6)
 
@@ -188,28 +189,28 @@ def solve(self, t=0, tstep=0, end_time=1000):
 
 if __name__ == '__main__':
     from time import time
-    N = (256, 512)
+    N = (64, 128)
     d = {
         'N': N,
         'Ra': 1000000.,
         'Pr': 0.7,
-        'dt': 0.0025,
+        'dt': 0.025,
         'filename': f'RB_{N[0]}_{N[1]}',
         'conv': 1,
-        'modplot': 100,
-        'moderror': 100,
+        'modplot': 10,
+        'moderror': 10,
         'modsave': 100,
         #'bcT': (0.9+0.1*sp.sin(2*(y-tt)), 0),
         #'bcT': (0.9+0.1*sympy.sin(2*y), 0),
         'bcT': (1, 0),
         'family': 'C',
         'checkpoint': 100,
         #'padding_factor': 1,
-        'timestepper': 'IMEXRK222'
+        'timestepper': 'IMEXRK3'
         }
     c = RayleighBenard(**d)
     t, tstep = c.initialize(rand=0.001, from_checkpoint=False)
     t0 = time()
-    c.solve(t=t, tstep=tstep, end_time=100)
+    c.solve(t=t, tstep=tstep, end_time=15)
     print('Computing time %2.4f'%(time()-t0))
 

docs/demos/Functions/functions.do.txt

@@ -337,8 +337,8 @@ for xj in X[0]:
 We may alternatively plot on a uniform mesh
 
 !bc pycod
-X = T.local_mesh(broadcast=True, uniform=True)
-plt.contourf(X[0], X[1], u.backward(kind='uniform'))
+X = T.local_mesh(bcast=True, kind='uniform')
+plt.contourf(X[0], X[1], u.backward(mesh='uniform'))
 !ec
 
 ===== Curvilinear coordinates =====

docs/demos/Integration/surfaceintegration.do.txt

@@ -138,7 +138,7 @@ spiral below, the result looks reasonable.
 !bc pycod
 import matplotlib.pyplot as plt
 fig = plt.figure(figsize=(4, 3))
-X = C.cartesian_mesh(uniform=True)
+X = C.cartesian_mesh(kind='uniform')
 ax = fig.add_subplot(111, projection='3d')
 p = ax.plot(X[0], X[1], X[2], 'r')
 hx = ax.set_xticks(np.linspace(-1, 1, 5))

docs/demos/RayleighBenard/Makefile

@@ -22,6 +22,11 @@ ${name}.rst: ${name}.do.txt
 	doconce format sphinx ${name} --sphinx_preserve_bib_keys
 	doconce subst 'XXX' '   ' $(basename $@).rst
 
+${name}.ipynb: ${name}.do.txt
+	doconce format ipynb $(basename $@).do.txt
+	#python add_metadata.py $(basename $@).ipynb
+	#jupyter nbconvert --inplace --execute ${basename $@}.ipynb
+
 view:
 	open ${name}.pdf &