Rebuild docs

demo/ChannelFlow2D.py

@@ -63,7 +63,7 @@ def __init__(self,
                  modsave=1e8,
                  moderror=100,
                  checkpoint=1000,
-                 timestepper='IMEXRK3'):
+                 timestepper='IMEXRK222'):
         self.N = N
         self.nu = nu
         self.dt = dt
@@ -78,7 +78,7 @@ def __init__(self,
         self.im1 = None
 
         # Regular spaces
-        self.B0 = FunctionSpace(N[0], family, bc=(0, 0, 0, 0), domain=domain[0])
+        self.B0 = FunctionSpace(N[0], family, bc={'left': {'D': 0, 'N': 0}, 'right': {'D': 0, 'N': 0}}, domain=domain[0])
         self.D0 = FunctionSpace(N[0], family, bc=(0, 0), domain=domain[0])
         self.C0 = FunctionSpace(N[0], family, domain=domain[0])
         self.F1 = FunctionSpace(N[1], 'F', dtype='d', domain=domain[1])
@@ -115,7 +115,7 @@ def __init__(self,
             self.dvdy = Project(Dx(self.u_[1], 1, 1), self.TD)
 
         self.curl = Project(curl(self.u_), self.TC)
-        self.divu = Project(div(self.u_), self.TC)
+        self.divu = Project(div(self.u_), self.TD)
         self.solP = None # For computing pressure
 
         # File for storing the results
@@ -130,7 +130,8 @@ def __init__(self,
         v = TestFunction(self.TB)
 
         # 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
+        #sol1 = chebyshev.la.Biharmonic if self.B0.family() == 'chebyshev' else la.SolverGeneric1ND
+        sol1 = la.SolverGeneric1ND
 
         self.pdes = {
 
@@ -175,6 +176,7 @@ def convection(self):
             dvdyp = self.dvdy().backward(padding_factor=self.padding_factor).v
             H[0] = self.TDp.forward(up[0]*dudxp+up[1]*dudyp, H[0])
             H[1] = self.TDp.forward(up[0]*dvdxp+up[1]*dvdyp, H[1])
+            #H[1] = self.TDp.forward(up[0]*dvdxp-up[1]*dudxp, H[1])
 
         elif self.conv == 1:
             curl = self.curl().backward(padding_factor=self.padding_factor)
@@ -204,13 +206,14 @@ def compute_pressure(self):
             self.d2udx2 = Project(self.nu*Dx(self.u_[0], 0, 2), self.TC)
             N0 = self.N0 = FunctionSpace(self.N[0], self.B0.family(), bc={'left': {'N': self.d2udx2()}, 'right': {'N': self.d2udx2()}})
             TN = self.TN = TensorProductSpace(comm, (N0, self.F1), modify_spaces_inplace=True)
-            sol = chebyshev.la.Helmholtz if self.B0.family() == 'chebyshev' else la.SolverGeneric1ND
+            #sol = chebyshev.la.Helmholtz if self.B0.family() == 'chebyshev' else la.SolverGeneric1ND
+            sol = la.SolverGeneric1ND
             self.divH = Inner(TestFunction(TN), -div(self.H_))
             self.solP = sol(inner(TestFunction(TN), div(grad(TrialFunction(TN)))))
             self.p_ = Function(TN)
 
         self.d2udx2()
-        self.N0.bc.set_tensor_bcs(self.N0, self.TN)
+        self.N0.bc.update(update_tensor=True)
         p_ = self.solP(self.divH(), self.p_, constraints=((0, 0),))
         return p_
 

demo/OrrSommerfeld2D.py

@@ -118,7 +118,7 @@ def print_energy_and_divergence(self, t, tstep):
         'family': 'C',
         'checkpoint': 10000000,
         'padding_factor': 1,
-        'timestepper': 'IMEXRK443'
+        'timestepper': 'IMEXRK222'
         }
     OS = True
     c = OrrSommerfeld(**d)

demo/OrrSommerfeld_eigs.py

@@ -26,10 +26,12 @@
 x = sp.Symbol('x', real=True)
 
 class OrrSommerfeld:
-    def __init__(self, alfa=1., Re=8000., N=80, quad='GC', **kwargs):
-        kwargs.update(dict(alfa=alfa, Re=Re, N=N, quad=quad))
+    def __init__(self, alfa=1., Re=8000., N=80, quad='GC', test='G', trial='G', **kwargs):
+        kwargs.update(dict(alfa=alfa, Re=Re, N=N, quad=quad, test=test, trial=trial))
         vars(self).update(kwargs)
         self.x, self.w = None, None
+        assert self.trial in ('PG', 'G')
+        assert self.test in ('PG', 'G')
 
     def interp(self, y, eigvals, eigvectors, eigval=1, verbose=False):
         """Interpolate solution eigenvector and it's derivative onto y
@@ -50,20 +52,29 @@ def interp(self, y, eigvals, eigvectors, eigval=1, verbose=False):
                 Print information or not
         """
         nx, eigval = self.get_eigval(eigval, eigvals, verbose)
-        SB = FunctionSpace(self.N, 'C', bc=(0, 0, 0, 0), quad=self.quad, dtype='D')
+        if self.trial == 'G':
+            SB = FunctionSpace(self.N, 'C', bc=(0, 0, 0, 0), quad=self.quad, dtype='D')
+        else:
+            SB = FunctionSpace(self.N, 'C', basis='Phi2', quad=self.quad, dtype='D')
         phi_hat = Function(SB)
         phi_hat[:-4] = np.squeeze(eigvectors[:, nx])
         phi = phi_hat.eval(y)
         dphidy = Dx(phi_hat, 0, 1).eval(y)
         return eigval, phi, dphidy
 
-    def assemble(self):
-        N = self.N
-        SB = FunctionSpace(N, 'C', bc=(0, 0, 0, 0), quad=self.quad)
-        test = SB.get_testspace('PG')
-        #test = FunctionSpace(N+4, 'C', basis='Compact4', quad=self.quad)
+    def get_trialspace(self, trial, dtype='d'):
+        if trial == 'G':
+            return FunctionSpace(self.N, 'C', basis='ShenBiharmonic', quad=self.quad, dtype=dtype)
+        return FunctionSpace(self.N, 'C', basis='Phi2', quad=self.quad, dtype=dtype)
+
+    def get_testspace(self, trial):
+        return trial.get_testspace(self.test)
+
+    def assemble(self, scale=None):
+        trial = self.get_trialspace(self.trial)
+        test = trial.get_testspace(self.test)
         v = TestFunction(test)
-        u = TrialFunction(SB)
+        u = TrialFunction(trial)
 
         # (u'', v)_w
         K = inner(v, Dx(u, 0, 2))
@@ -84,19 +95,27 @@ def assemble(self):
         a = self.alfa
         B = -Re*a*1j*(K-a**2*M)
         A = Q-2*a**2*K+(a**4 - 2*a*Re*1j)*M - 1j*a*Re*(K2-a**2*K1)
-
-        return A.diags().toarray(), B.diags().toarray()
-
-    def solve(self, verbose=False):
+        A, B = A.diags().toarray(), B.diags().toarray()
+        if scale is not None and not (scale[0] == 0 and scale[1] == 0):
+            assert isinstance(scale, tuple)
+            assert len(scale) == 2
+            s0, s1 = scale
+            k = np.arange(self.N-4)
+            testp = 1/(k+1)**(-s0) if s0 < 0 else (k+1)**s0
+            trialp = 1/(k+1)**(-s1) if s1 < 0 else (k+1)**s1
+            d =  testp[:, None] * trialp[None, :]
+            ##d = (1/A.diagonal())[:, None]
+            A *= d # * A
+            B *= d # * B
+        return A, B
+
+    def solve(self, verbose=False, scale=None):
         """Solve the Orr-Sommerfeld eigenvalue problem
         """
         if verbose:
             print('Solving the Orr-Sommerfeld eigenvalue problem...')
             print('Re = '+str(self.Re)+' and alfa = '+str(self.alfa))
-        A, B = self.assemble()
-        #d = (1/A.diagonal())[:, None]
-        #A *= d
-        #B *= d
+        A, B = self.assemble(scale=scale)
         return eig(A, B)
         # return eig(np.dot(inv(B), A))
 
@@ -136,6 +155,10 @@ def get_eigval(nx, eigvals, verbose=False):
                         help='Parameter')
     parser.add_argument('--quad', default='GC', type=str, choices=('GC', 'GL', 'LG'),
                         help='Discretization points: GC: Gauss-Chebyshev, GL: Gauss-Lobatto')
+    parser.add_argument('--test', default='G', type=str,
+                        help='G or PG, Galerkin or Petrov-Galerkin')
+    parser.add_argument('--trial', default='G', type=str,
+                        help='G or PG, Galerkin or Petrov-Galerkin')
     parser.add_argument('--plot', dest='plot', action='store_true', help='Plot eigenvalues')
     parser.add_argument('--verbose', dest='verbose', action='store_true', help='Print results')
     parser.set_defaults(plot=False)

demo/RayleighBenard2D.py

@@ -197,9 +197,9 @@ def solve(self, t=0, tstep=0, end_time=1000):
         'modplot': 10,
         'moderror': 10,
         'modsave': 1,
-        #'bcT': (0.9+0.1*sympy.sin(2*(y-tt)), 0),
+        'bcT': (0.9+0.1*sympy.sin(2*(y-tt)), 0),
         #'bcT': (0.9+0.1*sympy.sin(2*y), 0),
-        'bcT': (1, 0),
+        #'bcT': (1, 0),
         'family': 'C',
         'checkpoint': 100,
         #'padding_factor': 1,
@@ -208,7 +208,7 @@ def solve(self, t=0, tstep=0, end_time=1000):
     c = RayleighBenard(**d)
     t, tstep = c.initialize(rand=0.001, from_checkpoint=False)
     t0 = time()
-    c.solve(t=t, tstep=tstep, end_time=0.05)
+    c.solve(t=t, tstep=tstep, end_time=5)
     if comm.Get_rank() == 0:
         generate_xdmf('_'.join((d['filename'], 'U'))+'.h5', order='visit')
         generate_xdmf('_'.join((d['filename'], 'T'))+'.h5', order='visit')

docs/environment.yml

@@ -15,6 +15,7 @@ dependencies:
   - scipy
   - mpi4py
   - mpi4py-fft
+  - plotly
   - Sphinx==4.2
   - nbsphinx
   - pip:

docs/source/conf.py

@@ -21,7 +21,7 @@
 # -- Project information -----------------------------------------------------
 
 project = u'shenfun'
-copyright = u'2023, Mikael Mortensen'
+copyright = u'2024, Mikael Mortensen'
 author = u'Mikael Mortensen'
 
 p = subprocess.Popen(["git describe --tags | cut -d'-' -f 1"], stdout=subprocess.PIPE, shell=True)