Fixed vortex example

pySDC/implementations/problem_classes/HeatEquation_1D_FEniCS_matrix_forced.py

@@ -407,7 +407,6 @@ def solve_system(self, rhs, factor, u0, t):
         b = self.dtype_u(rhs)
 
         self.bc.apply(T, b.values.vector())
-        self.bc.apply(b.values.vector())
 
         df.solve(T, u.values.vector(), b.values.vector())
 

pySDC/implementations/problem_classes/VorticityVelocity_2D_FEniCS_periodic.py

@@ -3,7 +3,6 @@
 import dolfin as df
 import numpy as np
 
-from pySDC.core.Errors import ParameterError
 from pySDC.core.Problem import ptype
 from pySDC.implementations.datatype_classes.fenics_mesh import fenics_mesh, rhs_fenics_mesh
 
@@ -78,9 +77,9 @@ def __init__(self, c_nvars=None, family='CG', order=4, refinements=None, nu=0.01
             c_nvars = [(32, 32)]
 
         if refinements is None:
-            refinements = [1, 0]
+            refinements = 1
 
-        # Sub domain for Periodic boundary condition
+        # Subdomain for Periodic boundary condition
         class PeriodicBoundary(df.SubDomain):
             # Left boundary is "target domain" G
             def inside(self, x, on_boundary):
@@ -103,8 +102,8 @@ def map(self, x, y):
                     y[1] = x[1] - 1.0
 
         # set logger level for FFC and dolfin
-        df.set_log_level(df.WARNING)
         logging.getLogger('FFC').setLevel(logging.WARNING)
+        logging.getLogger('UFL').setLevel(logging.WARNING)
 
         # set solver and form parameters
         df.parameters["form_compiler"]["optimize"] = True
@@ -120,7 +119,8 @@ def map(self, x, y):
         # define function space for future reference
         self.V = df.FunctionSpace(mesh, family, order, constrained_domain=PeriodicBoundary())
         tmp = df.Function(self.V)
-        print('DoFs on this level:', len(tmp.vector().vector()[:]))
+        print('DoFs on this level:', len(tmp.vector()[:]))
+        self.fix_bc_for_residual = False
 
         # invoke super init, passing number of dofs, dtype_u and dtype_f
         super(fenics_vortex_2d, self).__init__(self.V)
@@ -162,7 +162,8 @@ def solve_system(self, rhs, factor, u0, t):
         """
 
         A = self.M + self.nu * factor * self.K
-        b = self.__apply_mass_matrix(rhs)
+        # b = self.apply_mass_matrix(rhs)
+        b = self.dtype_u(rhs)
 
         u = self.dtype_u(u0)
         df.solve(A, u.values.vector(), b.values.vector())
@@ -186,15 +187,17 @@ def __eval_fexpl(self, u, t):
             Explicit part of the right-hand side.
         """
 
-        A = 1.0 * self.K
-        b = self.__apply_mass_matrix(u)
+        A = self.K
+        b = self.apply_mass_matrix(u)
+        # b = self.dtype_u(u)
         psi = self.dtype_u(self.V)
         df.solve(A, psi.values.vector(), b.values.vector())
 
         fexpl = self.dtype_u(self.V)
         fexpl.values = df.project(
             df.Dx(psi.values, 1) * df.Dx(u.values, 0) - df.Dx(psi.values, 0) * df.Dx(u.values, 1), self.V
         )
+        fexpl = self.apply_mass_matrix(fexpl)
 
         return fexpl
 
@@ -215,9 +218,10 @@ def __eval_fimpl(self, u, t):
             Implicit part of the right-hand side.
         """
 
-        tmp = self.dtype_u(self.V)
-        tmp.values = df.Function(self.V, -1.0 * self.nu * self.K * u.values.vector())
-        fimpl = self.__invert_mass_matrix(tmp)
+        A = -self.nu * self.K
+        fimpl = self.dtype_u(self.V)
+        A.mult(u.values.vector(), fimpl.values.vector())
+        # fimpl = self.__invert_mass_matrix(fimpl)
 
         return fimpl
 
@@ -243,7 +247,7 @@ def eval_f(self, u, t):
         f.expl = self.__eval_fexpl(u, t)
         return f
 
-    def __apply_mass_matrix(self, u):
+    def apply_mass_matrix(self, u):
         r"""
         Routine to apply mass matrix.
 
@@ -258,7 +262,7 @@ def __apply_mass_matrix(self, u):
         """
 
         me = self.dtype_u(self.V)
-        me.values = df.Function(self.V, self.M * u.values.vector())
+        self.M.mult(u.values.vector(), me.values.vector())
 
         return me
 
@@ -278,12 +282,7 @@ def __invert_mass_matrix(self, u):
         """
 
         me = self.dtype_u(self.V)
-
-        A = 1.0 * self.M
-        b = self.dtype_u(u)
-
-        df.solve(A, me.values.vector(), b.values.vector())
-
+        df.solve(self.M, me.values.vector(), u.values.vector())
         return me
 
     def u_exact(self, t):

pySDC/playgrounds/FEniCS/vortex.py

@@ -1,9 +1,11 @@
 from pySDC.implementations.controller_classes.controller_nonMPI import controller_nonMPI
 from pySDC.implementations.problem_classes.VorticityVelocity_2D_FEniCS_periodic import fenics_vortex_2d
-from pySDC.implementations.sweeper_classes.generic_implicit import generic_implicit
+from pySDC.implementations.sweeper_classes.imex_1st_order_mass import imex_1st_order_mass
+from pySDC.implementations.sweeper_classes.imex_1st_order import imex_1st_order
 from pySDC.implementations.transfer_classes.TransferFenicsMesh import mesh_to_mesh_fenics
 
-if __name__ == "__main__":
+def setup_and_run(variant='mass'):
+
     num_procs = 1
 
     t0 = 0
@@ -12,7 +14,12 @@
 
     # initialize level parameters
     level_params = dict()
-    level_params['restol'] = 5e-09
+    if variant == 'mass':
+        level_params['restol'] = 5e-09 / 500
+    elif variant == 'mass_inv':
+        level_params['restol'] = 5e-09
+    else:
+        raise NotImplementedError('variant unknown')
     level_params['dt'] = dt
 
     # initialize step parameters
@@ -27,7 +34,8 @@
     sweeper_params = dict()
     sweeper_params['quad_type'] = 'RADAU-RIGHT'
     sweeper_params['num_nodes'] = [3]
-    sweeper_params['QI'] = 'LU'
+    sweeper_params['QI'] = 'MIN-SR-S'
+    sweeper_params['QE'] = 'PIC'
 
     problem_params = dict()
     problem_params['nu'] = 0.01
@@ -36,7 +44,7 @@
     problem_params['c_nvars'] = [(32, 32)]
     problem_params['family'] = 'CG'
     problem_params['order'] = [4]
-    problem_params['refinements'] = [1, 0]
+    problem_params['refinements'] = [1]
 
     # initialize controller parameters
     controller_params = dict()
@@ -46,7 +54,12 @@
     description = dict()
     description['problem_class'] = fenics_vortex_2d
     description['problem_params'] = problem_params
-    description['sweeper_class'] = generic_implicit
+    if variant == 'mass_inv':
+        description['sweeper_class'] = imex_1st_order
+    elif variant == 'mass':
+        description['sweeper_class'] = imex_1st_order_mass
+    else:
+        raise NotImplementedError('variant unknown')
     description['sweeper_params'] = sweeper_params
     description['level_params'] = level_params
     description['step_params'] = step_params
@@ -63,7 +76,11 @@
     # call main function to get things done...
     uend, stats = controller.run(u0=uinit, t0=t0, Tend=Tend)
 
-    # compute exact solution and compare
-    uex = P.u_exact(Tend)
+    # # compute exact solution and compare
+    # uex = P.u_exact(Tend)
+    #
+    # print('(classical) error at time %s: %s' % (Tend, abs(uex - uend) / abs(uex)))
 
-    print('(classical) error at time %s: %s' % (Tend, abs(uex - uend) / abs(uex)))
+if __name__ == "__main__":
+    setup_and_run(variant='mass')
+    # setup_and_run(variant='mass_inv')