From c042452cd37d87db0da07ddbecf5206edc4124a3 Mon Sep 17 00:00:00 2001 From: Robert Speck Date: Tue, 26 Dec 2023 14:05:28 +0100 Subject: [PATCH] Adding time-dependent BC example to FEniCS --- .gitignore | 2 + .../HeatEquation_1D_FEniCS_matrix_forced.py | 164 +++++++++++++++++- pySDC/playgrounds/FEniCS/heat_equation_M.py | 4 +- pySDC/tutorial/step_7/A_pySDC_with_FEniCS.py | 21 ++- 4 files changed, 182 insertions(+), 9 deletions(-) diff --git a/.gitignore b/.gitignore index 7ab0ead3bb..98df6abec9 100644 --- a/.gitignore +++ b/.gitignore @@ -156,3 +156,5 @@ Temporary Items # VSCode .vscode +*.cpp +pySDC/playgrounds/FEniCS/jitfailure-dolfin_expression_fc28530d435fa2de045af3312fc07c3b/recompile.sh diff --git a/pySDC/implementations/problem_classes/HeatEquation_1D_FEniCS_matrix_forced.py b/pySDC/implementations/problem_classes/HeatEquation_1D_FEniCS_matrix_forced.py index 6b7e028fba..55fac3edf1 100644 --- a/pySDC/implementations/problem_classes/HeatEquation_1D_FEniCS_matrix_forced.py +++ b/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 diff --git a/pySDC/playgrounds/FEniCS/heat_equation_M.py b/pySDC/playgrounds/FEniCS/heat_equation_M.py index db8fcda92b..b02db161f7 100644 --- a/pySDC/playgrounds/FEniCS/heat_equation_M.py +++ b/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: diff --git a/pySDC/tutorial/step_7/A_pySDC_with_FEniCS.py b/pySDC/tutorial/step_7/A_pySDC_with_FEniCS.py index 08ebedd172..7371c6a94c 100644 --- a/pySDC/tutorial/step_7/A_pySDC_with_FEniCS.py +++ b/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).