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 55fac3edf1..2ee5fc990d 100644 --- a/pySDC/implementations/problem_classes/HeatEquation_1D_FEniCS_matrix_forced.py +++ b/pySDC/implementations/problem_classes/HeatEquation_1D_FEniCS_matrix_forced.py @@ -407,7 +407,6 @@ def solve_system(self, rhs, factor, u0, t): b = self.dtype_u(rhs) self.bc.apply(T, b.values.vector()) - self.bc.apply(b.values.vector()) df.solve(T, u.values.vector(), b.values.vector()) diff --git a/pySDC/implementations/problem_classes/VorticityVelocity_2D_FEniCS_periodic.py b/pySDC/implementations/problem_classes/VorticityVelocity_2D_FEniCS_periodic.py index 969c7a0fea..689ce0b2e1 100644 --- a/pySDC/implementations/problem_classes/VorticityVelocity_2D_FEniCS_periodic.py +++ b/pySDC/implementations/problem_classes/VorticityVelocity_2D_FEniCS_periodic.py @@ -3,7 +3,6 @@ import dolfin as df import numpy as np -from pySDC.core.Errors import ParameterError from pySDC.core.Problem import ptype from pySDC.implementations.datatype_classes.fenics_mesh import fenics_mesh, rhs_fenics_mesh @@ -78,9 +77,9 @@ def __init__(self, c_nvars=None, family='CG', order=4, refinements=None, nu=0.01 c_nvars = [(32, 32)] if refinements is None: - refinements = [1, 0] + refinements = 1 - # Sub domain for Periodic boundary condition + # Subdomain for Periodic boundary condition class PeriodicBoundary(df.SubDomain): # Left boundary is "target domain" G def inside(self, x, on_boundary): @@ -103,8 +102,8 @@ def map(self, x, y): y[1] = x[1] - 1.0 # set logger level for FFC and dolfin - df.set_log_level(df.WARNING) logging.getLogger('FFC').setLevel(logging.WARNING) + logging.getLogger('UFL').setLevel(logging.WARNING) # set solver and form parameters df.parameters["form_compiler"]["optimize"] = True @@ -120,7 +119,8 @@ def map(self, x, y): # define function space for future reference self.V = df.FunctionSpace(mesh, family, order, constrained_domain=PeriodicBoundary()) tmp = df.Function(self.V) - print('DoFs on this level:', len(tmp.vector().vector()[:])) + print('DoFs on this level:', len(tmp.vector()[:])) + self.fix_bc_for_residual = False # invoke super init, passing number of dofs, dtype_u and dtype_f super(fenics_vortex_2d, self).__init__(self.V) @@ -162,7 +162,8 @@ def solve_system(self, rhs, factor, u0, t): """ A = self.M + self.nu * factor * self.K - b = self.__apply_mass_matrix(rhs) + # b = self.apply_mass_matrix(rhs) + b = self.dtype_u(rhs) u = self.dtype_u(u0) df.solve(A, u.values.vector(), b.values.vector()) @@ -186,8 +187,9 @@ def __eval_fexpl(self, u, t): Explicit part of the right-hand side. """ - A = 1.0 * self.K - b = self.__apply_mass_matrix(u) + 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()) @@ -195,6 +197,7 @@ def __eval_fexpl(self, u, t): 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 @@ -215,9 +218,10 @@ def __eval_fimpl(self, u, t): Implicit part of the right-hand side. """ - tmp = self.dtype_u(self.V) - tmp.values = df.Function(self.V, -1.0 * self.nu * self.K * u.values.vector()) - fimpl = self.__invert_mass_matrix(tmp) + 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) return fimpl @@ -243,7 +247,7 @@ def eval_f(self, u, t): f.expl = self.__eval_fexpl(u, t) return f - def __apply_mass_matrix(self, u): + def apply_mass_matrix(self, u): r""" Routine to apply mass matrix. @@ -258,7 +262,7 @@ def __apply_mass_matrix(self, u): """ me = self.dtype_u(self.V) - me.values = df.Function(self.V, self.M * u.values.vector()) + self.M.mult(u.values.vector(), me.values.vector()) return me @@ -278,12 +282,7 @@ def __invert_mass_matrix(self, u): """ me = self.dtype_u(self.V) - - A = 1.0 * self.M - b = self.dtype_u(u) - - df.solve(A, me.values.vector(), b.values.vector()) - + df.solve(self.M, me.values.vector(), u.values.vector()) return me def u_exact(self, t): diff --git a/pySDC/playgrounds/FEniCS/vortex.py b/pySDC/playgrounds/FEniCS/vortex.py index f0cabfe4dc..70f1dca888 100644 --- a/pySDC/playgrounds/FEniCS/vortex.py +++ b/pySDC/playgrounds/FEniCS/vortex.py @@ -1,9 +1,11 @@ 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.sweeper_classes.generic_implicit import generic_implicit +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 -if __name__ == "__main__": +def setup_and_run(variant='mass'): + num_procs = 1 t0 = 0 @@ -12,7 +14,12 @@ # initialize level parameters level_params = dict() - level_params['restol'] = 5e-09 + if variant == 'mass': + level_params['restol'] = 5e-09 / 500 + elif variant == 'mass_inv': + level_params['restol'] = 5e-09 + else: + raise NotImplementedError('variant unknown') level_params['dt'] = dt # initialize step parameters @@ -27,7 +34,8 @@ sweeper_params = dict() sweeper_params['quad_type'] = 'RADAU-RIGHT' sweeper_params['num_nodes'] = [3] - sweeper_params['QI'] = 'LU' + sweeper_params['QI'] = 'MIN-SR-S' + sweeper_params['QE'] = 'PIC' problem_params = dict() problem_params['nu'] = 0.01 @@ -36,7 +44,7 @@ problem_params['c_nvars'] = [(32, 32)] problem_params['family'] = 'CG' problem_params['order'] = [4] - problem_params['refinements'] = [1, 0] + problem_params['refinements'] = [1] # initialize controller parameters controller_params = dict() @@ -46,7 +54,12 @@ description = dict() description['problem_class'] = fenics_vortex_2d description['problem_params'] = problem_params - description['sweeper_class'] = generic_implicit + if variant == 'mass_inv': + description['sweeper_class'] = imex_1st_order + elif variant == 'mass': + description['sweeper_class'] = imex_1st_order_mass + else: + raise NotImplementedError('variant unknown') description['sweeper_params'] = sweeper_params description['level_params'] = level_params description['step_params'] = step_params @@ -63,7 +76,11 @@ # call main function to get things done... uend, stats = controller.run(u0=uinit, t0=t0, Tend=Tend) - # compute exact solution and compare - uex = P.u_exact(Tend) + # # compute exact solution and compare + # uex = P.u_exact(Tend) + # + # print('(classical) error at time %s: %s' % (Tend, abs(uex - uend) / abs(uex))) - print('(classical) error at time %s: %s' % (Tend, abs(uex - uend) / abs(uex))) +if __name__ == "__main__": + setup_and_run(variant='mass') + # setup_and_run(variant='mass_inv') \ No newline at end of file