Forgot to add different versions of vortex problem class

brownbaerchen · Dec 27, 2023 · 9bd2c2c · 9bd2c2c
1 parent 0206dd9
commit 9bd2c2c
Show file tree

Hide file tree

Showing 2 changed files with 169 additions and 13 deletions.
diff --git a/pySDC/implementations/problem_classes/VorticityVelocity_2D_FEniCS_periodic.py b/pySDC/implementations/problem_classes/VorticityVelocity_2D_FEniCS_periodic.py
@@ -22,7 +22,8 @@ class fenics_vortex_2d(ptype):
     .. math::
         \int_\Omega w_t v\,dx = - \nu \int_\Omega \nabla w \nabla v\,dx
 
-    The nonlinear system is solved in a *fully-implicit* way using Dolfin's weak solver.
+    This problem class treats the PDE in an IMEX way, with diffusion being the implicit part and everything else the explicit one.
+    The mass matrix needs inversion for this type of problem class, see the derived one for the mass-matrix version without inversion.
 
     Parameters
     ----------
@@ -162,9 +163,7 @@ def solve_system(self, rhs, factor, u0, t):
         """
 
         A = self.M + self.nu * factor * self.K
-        # b = self.apply_mass_matrix(rhs)
-        b = self.dtype_u(rhs)
-
+        b = self.apply_mass_matrix(rhs)
         u = self.dtype_u(u0)
         df.solve(A, u.values.vector(), b.values.vector())
 
@@ -187,17 +186,14 @@ def __eval_fexpl(self, u, t):
             Explicit part of the right-hand side.
         """
 
-        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())
+        df.solve(self.K, 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
 
@@ -221,7 +217,7 @@ def __eval_fimpl(self, u, t):
         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)
+        fimpl = self.__invert_mass_matrix(fimpl)
 
         return fimpl
 
@@ -326,3 +322,162 @@ def u_exact(self, t):
         # exit()
 
         return me
+
+
+class fenics_vortex_2d_mass(fenics_vortex_2d):
+    r"""
+    This class implements the vorticity-velocity problem in two dimensions with periodic boundary conditions
+    in :math:`[0, 1]^2`
+
+    .. math::
+        \frac{\partial w}{\partial t} = \nu \Delta w
+
+    for some parameter :math:`\nu`. In this class the problem is implemented that the spatial part is solved
+    using ``FEniCS`` [1]_. Hence, the problem is reformulated to the *weak formulation*
+
+    .. math::
+        \int_\Omega w_t v\,dx = - \nu \int_\Omega \nabla w \nabla v\,dx
+
+    This problem class treats the PDE in an IMEX way, with diffusion being the implicit part and everything else the explicit one.
+    No mass matrix inversion is needed here, i.e. using this problem class requires the imex_1st_order_mass sweeper.
+
+    Parameters
+    ----------
+    c_nvars : List of int tuple, optional
+        Spatial resolution, i.e., numbers of degrees of freedom in space, e.g. ``c_nvars=[(128, 128)]``.
+    family : str, optional
+        Indicates the family of elements used to create the function space
+        for the trail and test functions. The default is ``'CG'``, which are the class
+        of Continuous Galerkin, a *synonym* for the Lagrange family of elements, see [2]_.
+    order : int, optional
+        Defines the order of the elements in the function space.
+    refinements : int, optional
+        Denotes the refinement of the mesh. ``refinements=2`` refines the mesh by factor :math:`2`.
+    nu : float, optional
+        Diffusion coefficient :math:`\nu`.
+    rho : int, optional
+        Problem parameter.
+    delta : float, optional
+        Problem parameter.
+
+    Attributes
+    ----------
+    V : FunctionSpace
+        Defines the function space of the trial and test functions.
+    M : scalar, vector, matrix or higher rank tensor
+        Mass matrix for FENiCS.
+    K : scalar, vector, matrix or higher rank tensor
+        Stiffness matrix including diffusion coefficient (and correct sign).
+
+    References
+    ----------
+    .. [1] The FEniCS Project Version 1.5. M. S. Alnaes, J. Blechta, J. Hake, A. Johansson, B. Kehlet, A. Logg,
+        C. Richardson, J. Ring, M. E. Rognes, G. N. Wells. Archive of Numerical Software (2015).
+    .. [2] Automated Solution of Differential Equations by the Finite Element Method. A. Logg, K.-A. Mardal, G. N.
+        Wells and others. Springer (2012).
+    """
+
+    def solve_system(self, rhs, factor, u0, t):
+        r"""
+        Dolfin's linear solver for :math:`(M - factor \cdot A)\vec{u} = \vec{rhs}`.
+
+        Parameters
+        ----------
+        rhs : dtype_f
+            Right-hand side for the nonlinear system.
+        factor : float
+            Abbrev. for the node-to-node stepsize (or any other factor required).
+        u0 : dtype_u
+            Initial guess for the iterative solver (not used here so far).
+        t : float
+            Current time.
+
+        Returns
+        -------
+        u : dtype_u
+            The solution as mesh.
+        """
+
+        A = self.M + self.nu * factor * self.K
+        b = self.dtype_u(rhs)
+        u = self.dtype_u(u0)
+        df.solve(A, u.values.vector(), b.values.vector())
+
+        return u
+
+    def __eval_fexpl(self, u, t):
+        """
+        Helper routine to evaluate the explicit part of the right-hand side.
+
+        Parameters
+        ----------
+        u : dtype_u
+            Current values of the numerical solution.
+        t : float
+            Current time at which the numerical solution is computed.
+
+        Returns
+        -------
+        fexpl : dtype_u
+            Explicit part of the right-hand side.
+        """
+
+        b = self.apply_mass_matrix(u)
+        psi = self.dtype_u(self.V)
+        df.solve(self.K, 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
+
+    def __eval_fimpl(self, u, t):
+        """
+        Helper routine to evaluate the implicit part of the right-hand side.
+
+        Parameters
+        ----------
+        u : dtype_u
+            Current values of the numerical solution.
+        t : float
+            Current time at which the numerical solution is computed.
+
+        Returns
+        -------
+        fimpl : dtype_u
+            Implicit part of the right-hand side.
+        """
+
+        A = -self.nu * self.K
+        fimpl = self.dtype_u(self.V)
+        A.mult(u.values.vector(), fimpl.values.vector())
+
+        return fimpl
+
+    def eval_f(self, u, t):
+        """
+        Routine to evaluate both parts of the right-hand side.
+
+        Note: Need to add this here, because otherwise the parent class will call the "local" functions __eval_*
+        and not the ones of the child class.
+
+        Parameters
+        ----------
+        u : dtype_u
+            Current values of the numerical solution.
+        t : float
+            Current time at which the numerical solution is computed.
+
+        Returns
+        -------
+        f : dtype_f
+            The right-hand side divided into two parts.
+        """
+
+        f = self.dtype_f(self.V)
+        f.impl = self.__eval_fimpl(u, t)
+        f.expl = self.__eval_fexpl(u, t)
+        return f
diff --git a/pySDC/playgrounds/FEniCS/vortex.py b/pySDC/playgrounds/FEniCS/vortex.py
@@ -1,5 +1,5 @@
 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.problem_classes.VorticityVelocity_2D_FEniCS_periodic import fenics_vortex_2d, fenics_vortex_2d_mass
 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
@@ -10,7 +10,7 @@ def setup_and_run(variant='mass'):
 
     t0 = 0
     dt = 0.001
-    Tend = 4 * dt
+    Tend = 1 * dt
 
     # initialize level parameters
     level_params = dict()
@@ -52,11 +52,12 @@ def setup_and_run(variant='mass'):
 
     # Fill description dictionary for easy hierarchy creation
     description = dict()
-    description['problem_class'] = fenics_vortex_2d
     description['problem_params'] = problem_params
     if variant == 'mass_inv':
+        description['problem_class'] = fenics_vortex_2d
         description['sweeper_class'] = imex_1st_order
     elif variant == 'mass':
+        description['problem_class'] = fenics_vortex_2d_mass
         description['sweeper_class'] = imex_1st_order_mass
     else:
         raise NotImplementedError('variant unknown')
@@ -83,4 +84,4 @@ def setup_and_run(variant='mass'):
 
 if __name__ == "__main__":
     setup_and_run(variant='mass')
-    # setup_and_run(variant='mass_inv')
+    setup_and_run(variant='mass_inv')