Moving residual fix from datatype to problem class

Parallel-in-Time · Dec 26, 2023 · 61145dc · 61145dc
1 parent d337dc9
commit 61145dc
Show file tree

Hide file tree

Showing 3 changed files with 27 additions and 13 deletions.
diff --git a/pySDC/implementations/datatype_classes/fenics_mesh.py b/pySDC/implementations/datatype_classes/fenics_mesh.py
@@ -108,11 +108,6 @@ def __abs__(self):
         # return maximum
         return absval
 
-    def fix_bc(self, bc):
-        bc.apply(self.values.vector())
-
-        return None
-
 
 class rhs_fenics_mesh(object):
     """

diff --git a/pySDC/implementations/problem_classes/HeatEquation_1D_FEniCS_matrix_forced.py b/pySDC/implementations/problem_classes/HeatEquation_1D_FEniCS_matrix_forced.py
@@ -70,10 +70,6 @@ class fenics_heat(ptype):
         The forcing term :math:`f` in the heat equation.
     bc : DirichletBC
         Denotes the 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
     ----------
@@ -129,10 +125,9 @@ def Boundary(x, on_boundary):
         self.M = df.assemble(a_M)
         self.K = df.assemble(a_K)
 
-        # set boundary values, incl. homogeneous ones for the residual
+        # set boundary values
         self.bc = df.DirichletBC(self.V, df.Constant(c), Boundary)
         self.bc_hom = df.DirichletBC(self.V, df.Constant(0), Boundary)
-        self.fix_bc_for_residual = True
 
         # set forcing term as expression
         self.g = df.Expression(
@@ -360,7 +355,11 @@ class fenics_heat_mass(fenics_heat):
     g : Expression
         The forcing term :math:`f` in the heat equation.
     bc : DirichletBC
-        Denotes the Dirichlet boundary conditions (currently not used here).
+        Denotes the 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
     ----------
@@ -370,6 +369,13 @@ class fenics_heat_mass(fenics_heat):
         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"""
+
+        super().__init__(c_nvars, t0, family, order, refinements, nu, c)
+
+        self.fix_bc_for_residual = True
+
     def solve_system(self, rhs, factor, u0, t):
         r"""
         Dolfin's linear solver for :math:`(M - factor A) \vec{u} = \vec{rhs}`.
@@ -428,3 +434,15 @@ def eval_f(self, u, t):
         f.expl = self.apply_mass_matrix(f.expl)
 
         return f
+
+    def fix_residual(self, res):
+        """
+        Applies homogeneous Dirichlet boundary conditions to the residual
+
+        Parameters
+        ----------
+        res : dtype_u
+              Residual
+        """
+        self.bc_hom.apply(res.values.vector())
+        return None
diff --git a/pySDC/implementations/sweeper_classes/imex_1st_order_mass.py b/pySDC/implementations/sweeper_classes/imex_1st_order_mass.py
@@ -125,8 +125,9 @@ def compute_residual(self, stage=None):
             # add tau if associated
             if L.tau[m] is not None:
                 res[m] += L.tau[m]
+            # Due to different boundary conditions we might have to fix the residual
             if L.prob.fix_bc_for_residual:
-                res[m].fix_bc(L.prob.bc_hom)
+                L.prob.fix_residual(res[m])
             # use abs function from data type here
             res_norm.append(abs(res[m]))