diff --git a/pySDC/implementations/problem_classes/VorticityVelocity_2D_FEniCS_periodic.py b/pySDC/implementations/problem_classes/VorticityVelocity_2D_FEniCS_periodic.py index 689ce0b2e1..f48f865ec2 100644 --- 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 index 70f1dca888..7b955d5c47 100644 --- 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') \ No newline at end of file + setup_and_run(variant='mass_inv') \ No newline at end of file