Adding time-dependent BC example to FEniCS

.gitignore

@@ -156,3 +156,5 @@ Temporary Items
 
 # VSCode
 .vscode
+*.cpp
+pySDC/playgrounds/FEniCS/jitfailure-dolfin_expression_fc28530d435fa2de045af3312fc07c3b/recompile.sh

pySDC/implementations/problem_classes/HeatEquation_1D_FEniCS_matrix_forced.py

@@ -18,7 +18,7 @@ class fenics_heat(ptype):
     for :math:`x \in \Omega:=[0,1]`, where the forcing term :math:`f` is defined by
 
     .. math::
-        f(x, t) = -\cos(\pi x) (\sin(t) - \nu \pi^2 \cos(t)).
+        f(x, t) = -\sin(\pi x) (\sin(t) - \nu \pi^2 \cos(t)).
 
     For initial conditions with constant c and
 
@@ -311,12 +311,17 @@ class fenics_heat_mass(fenics_heat):
     for :math:`x \in \Omega:=[0,1]`, where the forcing term :math:`f` is defined by
 
     .. math::
-        f(x, t) = -\cos(\pi x) (\sin(t) - \nu \pi^2 \cos(t)).
+        f(x, t) = -\sin(\pi x) (\sin(t) - \nu \pi^2 \cos(t)).
 
-    The exact solution of the problem is
+    For initial conditions with constant c and
 
     .. math::
-        u(x, t) = \cos(\pi x)\cos(t).
+        u(x, 0) = \sin(\pi x) + c
+
+    the exact solution of the problem is given by
+
+    .. math::
+        u(x, t) = \sin(\pi x)\cos(t) + c.
 
     In this class the problem is implemented in the way that the spatial part is solved using ``FEniCS`` [1]_. Hence, the problem
     is reformulated to the *weak formulation*
@@ -446,3 +451,154 @@ def fix_residual(self, res):
         """
         self.bc_hom.apply(res.values.vector())
         return None
+
+
+# noinspection PyUnusedLocal
+class fenics_heat_mass_timebc(fenics_heat_mass):
+    r"""
+    Example implementing the forced one-dimensional heat equation with time-dependent Dirichlet boundary conditions
+
+    .. math::
+        \frac{d u}{d t} = \nu \frac{d^2 u}{d x^2} + f
+
+    for :math:`x \in \Omega:=[0,1]`, where the forcing term :math:`f` is defined by
+
+    .. math::
+        f(x, t) = -\cos(\pi x) (\sin(t) - \nu \pi^2 \cos(t)).
+
+    and the boundary conditions are given by
+
+    .. math::
+        u(x, t) = \cos(\pi x)\cos(t).
+
+    The exact solution of the problem is given by
+
+    .. math::
+        u(x, t) = \cos(\pi x)\cos(t) + c.
+
+    In this class the problem is implemented in the way that the spatial part is solved using ``FEniCS`` [1]_. Hence, the problem
+    is reformulated to the *weak formulation*
+
+    .. math:
+        \int_\Omega u_t v\,dx = - \nu \int_\Omega \nabla u \nabla v\,dx + \int_\Omega f v\,dx.
+
+    The forcing term is treated explicitly, and is expressed via the mass matrix resulting from the left-hand side term
+    :math:`\int_\Omega u_t v\,dx`, and the other part will be treated in an implicit way.
+
+    Parameters
+    ----------
+    c_nvars : int, optional
+        Spatial resolution, i.e., numbers of degrees of freedom in space.
+    t0 : float, optional
+        Starting time.
+    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`.
+
+    Attributes
+    ----------
+    V : FunctionSpace
+        Defines the function space of the trial and test functions.
+    M : scalar, vector, matrix or higher rank tensor
+        Denotes the expression :math:`\int_\Omega u_t v\,dx`.
+    K : scalar, vector, matrix or higher rank tensor
+        Denotes the expression :math:`- \nu \int_\Omega \nabla u \nabla v\,dx`.
+    g : Expression
+        The forcing term :math:`f` in the heat equation.
+    bc : DirichletBC
+        Denotes the time-dependent Dirichlet boundary conditions.
+    bc_hom : DirichletBC
+        Denotes the homogeneous Dirichlet boundary conditions, potentially required for fixing the residual
+    fix_bc_for_residual: boolean
+        flag to indicate that the residual requires special treatment due to boundary conditions
+
+    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 __init__(self, c_nvars=128, t0=0.0, family='CG', order=4, refinements=1, nu=0.1, c=0.0):
+        """Initialization routine"""
+
+        # define the Dirichlet boundary
+        def Boundary(x, on_boundary):
+            return on_boundary
+
+        super().__init__(c_nvars, t0, family, order, refinements, nu, c)
+
+        self.u_D = df.Expression('cos(a*x[0]) * cos(t) + c', c=self.c, a=np.pi, t=t0, degree=self.order)
+        self.bc = df.DirichletBC(self.V, self.u_D, Boundary)
+        self.bc_hom = df.DirichletBC(self.V, df.Constant(0), Boundary)
+
+        # set forcing term as expression
+        self.g = df.Expression(
+            '-cos(a*x[0]) * (sin(t) - b*a*a*cos(t))',
+            a=np.pi,
+            b=self.nu,
+            t=self.t0,
+            degree=self.order,
+        )
+
+    def solve_system(self, rhs, factor, u0, t):
+        r"""
+        Dolfin's linear solver for :math:`(M - factor 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
+            Solution.
+        """
+
+        u = self.dtype_u(u0)
+        T = self.M - factor * self.K
+        b = self.dtype_u(rhs)
+
+        self.u_D.t = t
+
+        self.bc.apply(T, b.values.vector())
+        self.bc.apply(b.values.vector())
+
+        df.solve(T, u.values.vector(), b.values.vector())
+
+        return u
+
+    def u_exact(self, t):
+        r"""
+        Routine to compute the exact solution at time :math:`t`.
+
+        Parameters
+        ----------
+        t : float
+            Time of the exact solution.
+
+        Returns
+        -------
+        me : dtype_u
+            Exact solution.
+        """
+
+        u0 = df.Expression('cos(a*x[0]) * cos(t) + c', c=self.c, a=np.pi, t=t, degree=self.order)
+        me = self.dtype_u(df.interpolate(u0, self.V), val=self.V)
+
+        return me

pySDC/playgrounds/FEniCS/heat_equation_M.py

@@ -2,7 +2,7 @@
 
 
 from pySDC.implementations.controller_classes.controller_nonMPI import controller_nonMPI
-from pySDC.implementations.problem_classes.HeatEquation_1D_FEniCS_matrix_forced import fenics_heat_mass
+from pySDC.implementations.problem_classes.HeatEquation_1D_FEniCS_matrix_forced import fenics_heat_mass, fenics_heat_mass_timebc
 from pySDC.implementations.sweeper_classes.imex_1st_order_mass import imex_1st_order_mass
 from pySDC.implementations.transfer_classes.BaseTransfer_mass import base_transfer_mass
 from pySDC.implementations.transfer_classes.TransferFenicsMesh import mesh_to_mesh_fenics
@@ -55,7 +55,7 @@ def run_simulation(ml=None, mass=None):
     # Fill description dictionary for easy hierarchy creation
     description = dict()
     if mass:
-        description['problem_class'] = fenics_heat_mass
+        description['problem_class'] = fenics_heat_mass_timebc
         description['sweeper_class'] = imex_1st_order_mass
         description['base_transfer_class'] = base_transfer_mass
     else:

pySDC/tutorial/step_7/A_pySDC_with_FEniCS.py

@@ -4,7 +4,11 @@
 from pySDC.helpers.stats_helper import get_sorted
 
 from pySDC.implementations.controller_classes.controller_nonMPI import controller_nonMPI
-from pySDC.implementations.problem_classes.HeatEquation_1D_FEniCS_matrix_forced import fenics_heat_mass, fenics_heat
+from pySDC.implementations.problem_classes.HeatEquation_1D_FEniCS_matrix_forced import (
+    fenics_heat_mass,
+    fenics_heat,
+    fenics_heat_mass_timebc,
+)
 from pySDC.implementations.sweeper_classes.imex_1st_order_mass import imex_1st_order_mass, imex_1st_order
 from pySDC.implementations.transfer_classes.TransferFenicsMesh import mesh_to_mesh_fenics
 
@@ -98,6 +102,11 @@ def run_variants(variant=None, ml=None, num_procs=None):
     elif variant == 'mass_inv':
         description['problem_class'] = fenics_heat
         description['sweeper_class'] = imex_1st_order
+    elif variant == 'mass_timebc':
+        # Can increase the tolerance here, errors are higher anyway
+        description['level_params']['restol'] *= 20
+        description['problem_class'] = fenics_heat_mass_timebc
+        description['sweeper_class'] = imex_1st_order_mass
     else:
         raise NotImplementedError('Variant %s is not implemented' % variant)
 
@@ -146,10 +155,13 @@ def run_variants(variant=None, ml=None, num_procs=None):
 
     if num_procs == 1:
         assert np.mean(niters) <= 6.0, 'Mean number of iterations is too high, got %s' % np.mean(niters)
-        assert err <= 4.1e-08, 'Error is too high, got %s' % err
+        if variant == 'mass' or variant == 'mass_inv':
+            assert err <= 1.14e-08, 'Error is too high, got %s' % err
+        else:
+            assert err <= 3.25e-07, 'Error is too high, got %s' % err
     else:
         assert np.mean(niters) <= 11.6, 'Mean number of iterations is too high, got %s' % np.mean(niters)
-        assert err <= 4.0e-08, 'Error is too high, got %s' % err
+        assert err <= 1.14e-08, 'Error is too high, got %s' % err
 
     f.write('\n')
     print()
@@ -159,9 +171,12 @@ def run_variants(variant=None, ml=None, num_procs=None):
 def main():
     run_variants(variant='mass_inv', ml=False, num_procs=1)
     run_variants(variant='mass', ml=False, num_procs=1)
+    run_variants(variant='mass_timebc', ml=False, num_procs=1)
     run_variants(variant='mass_inv', ml=True, num_procs=1)
     run_variants(variant='mass', ml=True, num_procs=1)
+    run_variants(variant='mass_timebc', ml=True, num_procs=1)
     run_variants(variant='mass_inv', ml=True, num_procs=5)
+
     # WARNING: all other variants do NOT work, either because of FEniCS restrictions (weak forms with different meshes
     # will not work together) or because of inconsistent use of the mass matrix (locality condition for the tau
     # correction is not satisfied, mass matrix does not permute with restriction).