# --- # jupyter: # jupytext: # formats: py:percent # kernelspec: # display_name: Python 3 # name: python3 # language_info: # name: python # --- # %% [markdown] # # Part A: pySDC and FEniCS # # In this example, pySDC is coupled with the [FEniCS framework](https://fenicsproject.org/) for finite elements in # space. This implies significant changes to the algorithm, depending on whether or not the mass matrix should be # inverted. SDC, MLSDC and PFASST can be used without changes when the right-hand side of the ODE is defined with the # inverse of the mass matrix. Otherwise, the mass matrix has to be used, e.g. in the tau-correction. This example tests # different variants of this methodology for SDC, MLSDC and PFASST. # # ## The setup # # The forced heat equation in 1D, with continuous Lagrange elements of order 4, and for MLSDC a second level on a # mesh half as fine, with the same element order and the same collocation nodes: coarsening in the mesh is where the # multilevel gain is, while coarsening in the nodes or in the element order throws it away. The problem and sweeper # classes are set per variant below. # %% from pathlib import Path import numpy as np 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, 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 from pySDC.implementations.transfer_classes.BaseTransfer_mass import base_transfer_mass def setup(t0=None, ml=None): """ Helper routine to set up parameters Args: t0 (float): initial time ml (bool): use single or multiple levels Returns: description and controller_params parameter dictionaries """ # initialize level parameters level_params = dict() level_params['restol'] = 5e-10 level_params['dt'] = 0.2 # initialize step parameters step_params = dict() step_params['maxiter'] = 20 # initialize sweeper parameters sweeper_params = dict() sweeper_params['quad_type'] = 'RADAU-RIGHT' if ml: # Keep the collocation nodes on both levels. Coarsening them does not help: with M=1 the coarse # level is asymptotically inert, so MLSDC takes neither more nor fewer iterations than SDC. sweeper_params['num_nodes'] = [3, 3] else: sweeper_params['num_nodes'] = [3] problem_params = dict() problem_params['nu'] = 0.1 problem_params['t0'] = t0 # ugly, but necessary to set up this ProblemClass problem_params['c_nvars'] = [128] problem_params['family'] = 'CG' problem_params['c'] = 1.0 if ml: # Coarsen in the mesh only. This is where the multilevel gain actually is: it halves the number of # iterations (6 -> 3 here, and 11.6 -> 3.8 for PFASST below), for both the mass-inverse and the mass # formulation. Coarsening in the nodes and the element order instead throws that away. problem_params['order'] = [4, 4] problem_params['refinements'] = [1, 0] else: problem_params['order'] = [4] problem_params['refinements'] = [1] # initialize controller parameters controller_params = dict() controller_params['logger_level'] = 30 base_transfer_params = dict() base_transfer_params['finter'] = True # Fill description dictionary for easy hierarchy creation description = dict() description['problem_class'] = None description['problem_params'] = problem_params description['sweeper_class'] = None description['sweeper_params'] = sweeper_params description['level_params'] = level_params description['step_params'] = step_params description['space_transfer_class'] = mesh_to_mesh_fenics description['base_transfer_params'] = base_transfer_params return description, controller_params # %% [markdown] # ## The variants # # - `'mass_inv'`: the right-hand side includes the inverse of the mass matrix, with the problem class `fenics_heat`, # so that the standard IMEX sweeper works unchanged. # - `'mass'`: the mass matrix stays on the left, with `fenics_heat_mass` and the sweeper `imex_1st_order_mass`, which # applies it to the initial value and in the residual. Between the levels, `base_transfer_mass` restricts the # quantities that carry the mass matrix, the tau-correction and the initial value, as load vectors. The right-hand # side is one as well, so it is re-evaluated on the fine level instead of being interpolated (`finter=False`). # - `'mass_timebc'`: as `'mass'`, but with time-dependent boundary conditions, `fenics_heat_mass_timebc`. # # Each run prints the error, statistics of the iterations and the time to solution, and appends them to # `data/step_7_A_out.txt`. # %% def run_variants(variant=None, ml=None, num_procs=None): """ Main routine to run the different implementations of the heat equation with FEniCS Args: variant (str): specifies the variant ml (bool): use single or multiple levels num_procs (int): number of processors in time """ Tend = 1.0 t0 = 0.0 description, controller_params = setup(t0=t0, ml=ml) if variant == 'mass': # Note that we need to reduce the tolerance for the residual here, since otherwise the error will be too high description['level_params']['restol'] /= 500 description['problem_class'] = fenics_heat_mass description['sweeper_class'] = imex_1st_order_mass description['base_transfer_class'] = base_transfer_mass # prolong_f is not available for the mass formulation: f is a load vector there, so it cannot be # interpolated. base_transfer_mass falls back to prolong, which re-evaluates f on the fine level. description['base_transfer_params']['finter'] = False elif variant == 'mass_inv': description['problem_class'] = fenics_heat description['sweeper_class'] = imex_1st_order elif variant == 'mass_timebc': # Trades accuracy for iterations: converged this runs to 1.7e-07 in 9.4 iterations, and # stopping 20x earlier costs about a factor two in error to get back to 6. description['level_params']['restol'] *= 20 description['problem_class'] = fenics_heat_mass_timebc description['sweeper_class'] = imex_1st_order_mass description['base_transfer_class'] = base_transfer_mass description['base_transfer_params']['finter'] = False else: raise NotImplementedError('Variant %s is not implemented' % variant) # quickly generate block of steps controller = controller_nonMPI(num_procs=num_procs, controller_params=controller_params, description=description) # get initial values on finest level P = controller.MS[0].levels[0].prob uinit = P.u_exact(0.0) # 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) err = abs(uex - uend) / abs(uex) Path("data").mkdir(parents=True, exist_ok=True) f = open('data/step_7_A_out.txt', 'a') out = f'Variant {variant} with ml={ml} and num_procs={num_procs} -- error at time {Tend}: {err}' f.write(out + '\n') print(out) # filter statistics by type (number of iterations) iter_counts = get_sorted(stats, type='niter', sortby='time') niters = np.array([item[1] for item in iter_counts]) out = ' Mean number of iterations: %4.2f' % np.mean(niters) f.write(out + '\n') print(out) out = ' Range of values for number of iterations: %2i ' % np.ptp(niters) f.write(out + '\n') print(out) out = ' Position of max/min number of iterations: %2i -- %2i' % (int(np.argmax(niters)), int(np.argmin(niters))) f.write(out + '\n') print(out) out = ' Std and var for number of iterations: %4.2f -- %4.2f' % (float(np.std(niters)), float(np.var(niters))) f.write(out + '\n') print(out) timing = get_sorted(stats, type='timing_run', sortby='time') out = f'Time to solution: {timing[0][1]:6.4f} sec.' f.write(out + '\n') print(out) # Bounds are meant to catch a regression, not to pin the current numbers: at the committed # settings the errors are 1.14e-08, or 2.8-3.2e-07 for mass_timebc, whose time-dependent # boundary data makes it a harder problem; the iteration counts are 6.00 serial, 3.00-3.20 with # a coarse level and 3.80 on five parallel steps. max_err = 2e-08 if variant == 'mass_timebc': # the loosened tolerance stops the parallel run a little earlier still, at a larger error max_err = 5e-07 if num_procs == 1 else 2e-06 max_niter = (5.0 if ml else 8.0) if num_procs == 1 else 6.0 assert np.mean(niters) <= max_niter, 'Mean number of iterations is too high, got %s' % np.mean(niters) assert err <= max_err, 'Error is too high, got %s' % err f.write('\n') print() f.close() # %% [markdown] # SDC, MLSDC and PFASST with 5 steps in parallel, emulated in one process, each with all three variants. # %% 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) run_variants(variant='mass', ml=True, num_procs=5) run_variants(variant='mass_timebc', ml=True, num_procs=5) if __name__ == "__main__": main() # %% [markdown] # ## Results # # FEniCS does not run in the browser, nor in the environment this website is built in. These are the results of our # CI, which runs this part in an environment with FEniCS, in the run that built this page: # # :::{literalinclude} /../../data/step_7_A_out.txt # :language: text # ::: # # :::{admonition} Important things to note # - Even core routines can be replaced where a method needs it: for the mass-matrix formulation, pySDC also has # `base_transfer_mass`, which the mass variants here use for MLSDC and PFASST. # - The project [Finite elements, the mass-matrix route](../../projects/FEM_with_FEniCS) takes the mass-matrix # formulation further: nonlinear problems, discontinuous elements, three levels, and which coarsening pays. # - It is also valuable to check out the data type and transfer classes required to work with FEniCS. Both can be # found in the `implementations` folder. # :::