Coverage for pySDC/projects/FEM_with_FEniCS/setups.py: 100%
69 statements
« prev ^ index » next coverage.py v7.17.0, created at 2026-11-02 12:30 +0000
1"""
2Problem setups. Every one takes the mass-matrix route: ``eval_f`` returns the assembled weak form,
3``solve_system`` takes a right-hand side already in the dual space, and the mass matrix is never
4inverted. See the README for what the combinations are for.
5"""
7import math
9from pySDC.implementations.problem_classes.HeatEquation_1D_FEniCS_matrix_forced import fenics_heat_mass
10from pySDC.implementations.problem_classes.GrayScott_1D_FEniCS_implicit import fenics_grayscott_mass
11from pySDC.implementations.problem_classes.VorticityVelocity_2D_FEniCS_periodic import fenics_vortex_2d_mass
12from pySDC.projects.FEM_with_FEniCS.problem_classes.Burgers_1D_FEniCS import fenics_burgers_mass
13from pySDC.projects.FEM_with_FEniCS.problem_classes.DG_1D_FEniCS import (
14 fenics_heat_dg_mass,
15 fenics_burgers_dg_mass,
16 fenics_grayscott_dg_mass,
17)
18from pySDC.implementations.sweeper_classes.imex_1st_order_mass import imex_1st_order_mass
19from pySDC.implementations.sweeper_classes.generic_implicit_mass import generic_implicit_mass
20from pySDC.implementations.transfer_classes.BaseTransfer_mass import base_transfer_mass
21from pySDC.implementations.transfer_classes.TransferFenicsMesh import mesh_to_mesh_fenics
23EXAMPLES = ('heat', 'burgers', 'grayscott', 'vortex')
24FAMILIES = ('CG', 'DG')
25COARSENINGS = ('h', 'p')
27#: Interior penalty constant for the DG forms. Only has to clear the coercivity threshold; the
28#: multilevel behaviour is flat in it from 2 upwards.
29SIGMA = 10.0
31#: Per-example defaults. ``refinements`` (h-coarsening) and ``orders`` (p-coarsening) are indexed by
32#: level, so the fine level is always the first entry; everything else is shared across levels. Both
33#: ladders halve the dof count per level, so the two coarsening directions cost the same.
34_SETUPS = {
35 'heat': {
36 'problem_class': fenics_heat_mass,
37 'dg_problem_class': fenics_heat_dg_mass,
38 'sweeper_class': imex_1st_order_mass,
39 'problem_params': {'nu': 0.1, 't0': 0.0, 'c_nvars': 128, 'c': 1.0},
40 'sweeper_params': {'quad_type': 'RADAU-RIGHT', 'QI': 'LU'},
41 'refinements': [2, 1, 0],
42 'orders': [4, 2, 1],
43 'dt': 0.2,
44 'nsteps': 8,
45 'restol': 5e-10 / 500,
46 'maxiter': 20,
47 'num_nodes': 3,
48 'utol': 1e-8,
49 'pfasst_procs': (1, 2, 4, 8),
50 },
51 'grayscott': {
52 'problem_class': fenics_grayscott_mass,
53 'dg_problem_class': fenics_grayscott_dg_mass,
54 'sweeper_class': generic_implicit_mass,
55 'problem_params': {
56 'c_nvars': 64,
57 'newton_tol': 1e-12,
58 'newton_rtol': 1e-11,
59 },
60 'sweeper_params': {'quad_type': 'RADAU-RIGHT', 'QI': 'LU'},
61 'refinements': [2, 1, 0],
62 'orders': [4, 2, 1],
63 'dt': 1.0,
64 'nsteps': 8,
65 'restol': 1e-9,
66 'maxiter': 30,
67 'num_nodes': 3,
68 'utol': 1e-6,
69 'pfasst_procs': (1, 2, 4, 8),
70 },
71 'burgers': {
72 'problem_class': fenics_burgers_mass,
73 'dg_problem_class': fenics_burgers_dg_mass,
74 'sweeper_class': generic_implicit_mass,
75 'problem_params': {
76 'c_nvars': 64,
77 'nu': 0.02,
78 'newton_tol': 1e-12,
79 'newton_rtol': 1e-11,
80 },
81 'sweeper_params': {'quad_type': 'RADAU-RIGHT', 'QI': 'LU'},
82 'refinements': [2, 1, 0],
83 'orders': [4, 2, 1],
84 'dt': 0.02,
85 'nsteps': 8,
86 'restol': 1e-9,
87 'maxiter': 30,
88 'num_nodes': 3,
89 # PFASST on a nonlinear problem drifts from serial by more than the linear examples do
90 'utol': 1e-6,
91 'pfasst_procs': (1, 2, 4, 8),
92 },
93 # The vortex is the correctness example, not a savings example: MLSDC converges to the same
94 # answer but the coarse level does not pay for itself, 0.73x on two levels, and four checks
95 # failed to move that -- the preconditioner, the shear-layer width, the element order and a 16x
96 # refinement of the whole hierarchy. It earns its place for what it covers rather than what it
97 # saves: it is the only 2d and the only periodic hierarchy here, and while it was out of the
98 # suite the prolongation was rewritten in a way that extrapolated on periodic spaces, with
99 # nothing left to catch it. CG only, because the DG problem classes are 1d.
100 'vortex': {
101 'problem_class': fenics_vortex_2d_mass,
102 'dg_problem_class': None,
103 'sweeper_class': imex_1st_order_mass,
104 'problem_params': {'nu': 0.01, 'delta': 0.05, 'rho': 50, 'c_nvars': [(8, 8)]},
105 'sweeper_params': {'quad_type': 'RADAU-RIGHT', 'QI': 'LU'},
106 'families': ('CG',),
107 'coarsenings': ('h',),
108 'refinements': [2, 1, 0],
109 'orders': [2],
110 'dt': 0.001,
111 'nsteps': 8,
112 'restol': 5e-9 / 500,
113 'maxiter': 30,
114 'num_nodes': 3,
115 'utol': 1e-6,
116 # The order study does not apply: it asks whether a high-order coarse space buys
117 # savings, and this example has none to buy. Its base order is 2 anyway, so the CG4
118 # rung would need a negative refinement.
119 'order_study': (),
120 'pays_off': False,
121 # 8 steps used to report convergence and return an O(1) wrong answer here. That was the
122 # point sampling: with it the error against serial is 1.21e+00 at 22 iterations, with the
123 # L2 projection 3.9e-08 at 11.25.
124 'pfasst_procs': (1, 2, 4, 8),
125 },
126}
129#: Marks a key every example must define, as opposed to one with a default.
130_REQUIRED = object()
133def _setting(example, key, default=_REQUIRED):
134 """One per-example entry from :data:`_SETUPS`."""
135 if example not in _SETUPS:
136 raise ValueError(f'unknown example {example!r}, expected one of {EXAMPLES}')
137 return _SETUPS[example][key] if default is _REQUIRED else _SETUPS[example].get(key, default)
140def get_tolerance(example):
141 """Tolerance for comparing solutions of the same example across configurations."""
142 return _setting(example, 'utol')
145def get_pfasst_procs(example):
146 """Process counts PFASST is known to be reliable for on this example."""
147 return _setting(example, 'pfasst_procs')
150def get_families(example):
151 """Element families this example has a problem class for. Not every example has both."""
152 return _setting(example, 'families', FAMILIES)
155def get_coarsenings(example):
156 """Coarsening directions this example is set up for."""
157 return _setting(example, 'coarsenings', COARSENINGS)
160def get_order_study(example):
161 """Element orders to compare at equal dof counts, or empty where that has no meaning."""
162 return _setting(example, 'order_study', (1, 2, 4))
165def pays_off(example):
166 """Whether the coarse level is expected to cost less than it saves on this example."""
167 return _setting(example, 'pays_off', True)
170def get_description(
171 example, nlevels=1, family='CG', coarsening='h', order=None, dt=None, restol=None, maxiter=None, nsteps=None
172):
173 """
174 Build the pySDC description for one of the three examples.
176 Parameters
177 ----------
178 example : str
179 One of :data:`EXAMPLES`.
180 nlevels : int
181 1 for SDC, 2 or 3 for MLSDC/PFASST.
182 family : str
183 One of :data:`FAMILIES`: ``'CG'`` or ``'DG'``.
184 coarsening : str
185 One of :data:`COARSENINGS`: ``'h'`` coarsens the mesh at fixed element order, ``'p'`` drops
186 the element order on the fine mesh. Irrelevant for ``nlevels=1``.
187 order : int, optional
188 Element order for h-coarsening, overriding the default. The refinement ladder is shifted to
189 keep the dof count on every level, so that changing this isolates the *order* of the coarse
190 space from its size: in 1d a Lagrange space has ``order * cells + 1`` dofs, so halving the
191 order and doubling the cells leaves the count alone. Only meaningful with ``coarsening='h'``.
192 dt, restol, maxiter, nsteps : optional
193 Override the per-example defaults.
195 Returns
196 -------
197 description : dict
198 controller_params : dict
199 t0, Tend : float
200 """
201 if example not in _SETUPS:
202 raise ValueError(f'unknown example {example!r}, expected one of {EXAMPLES}')
203 if family not in FAMILIES:
204 raise ValueError(f'unknown family {family!r}, expected one of {FAMILIES}')
205 if coarsening not in COARSENINGS:
206 raise ValueError(f'unknown coarsening {coarsening!r}, expected one of {COARSENINGS}')
207 if family not in get_families(example):
208 raise ValueError(f'{example!r} has no {family} problem class, only {get_families(example)}')
209 if coarsening not in get_coarsenings(example):
210 raise ValueError(f'{example!r} is only set up for {get_coarsenings(example)}-coarsening')
211 if not 1 <= nlevels <= 3:
212 raise ValueError(f'nlevels must be 1, 2 or 3, got {nlevels}')
214 s = _SETUPS[example]
215 dt = s['dt'] if dt is None else dt
216 nsteps = s['nsteps'] if nsteps is None else nsteps
218 problem_params = dict(s['problem_params'], family=family)
219 if coarsening == 'h':
220 # per level: take the coarsest nlevels entries so the fine level stays fixed
221 refinements = s['refinements'][:nlevels]
222 element_order = s['orders'][0] if order is None else order
223 if order is not None:
224 # keep the dof count: halving the order and doubling the cells per dimension leaves it
225 shift = round(math.log2(s['orders'][0] / order))
226 refinements = [r + shift for r in refinements]
227 if min(refinements) < 0:
228 # dolfin just skips the refine loop, silently giving two identical levels
229 raise ValueError(
230 f'{example!r} cannot hold its dof count at order {order}: that needs refinements '
231 f'{refinements}, and a negative refinement is not a coarsening'
232 )
233 problem_params['refinements'] = refinements
234 problem_params['order'] = element_order
235 orders = [element_order] * nlevels
236 else:
237 if order is not None:
238 raise ValueError('the order override only applies to h-coarsening')
239 # same mesh everywhere, the fine one; the element order is what comes down
240 refinements = [s['refinements'][0]] * nlevels
241 problem_params['refinements'] = s['refinements'][0]
242 problem_params['order'] = s['orders'][:nlevels]
243 orders = s['orders'][:nlevels]
245 if family == 'DG':
246 # Pin the interior penalty to the fine level. The SIPG form depends on the mesh and order it
247 # is built on, so rediscretising it on a coarse level would change the penalty by h or p^2 --
248 # and the penalty is the dominant term, so the coarse operator would stop being the Galerkin
249 # operator P^T A_F P of the fine one. alpha_l = sigma p_0^2 h_l / h_0 keeps alpha_l / h_l
250 # equal on every level, which restores A_G = P^T A_F P exactly.
251 problem_params['penalty'] = [SIGMA * orders[0] ** 2 * 2 ** (refinements[0] - r) for r in refinements]
253 sweeper_params = dict(s['sweeper_params'])
254 sweeper_params['num_nodes'] = [s['num_nodes']] * nlevels
256 description = {
257 'problem_class': s['dg_problem_class'] if family == 'DG' else s['problem_class'],
258 'problem_params': problem_params,
259 'sweeper_class': s['sweeper_class'],
260 'sweeper_params': sweeper_params,
261 'level_params': {'restol': s['restol'] if restol is None else restol, 'dt': dt},
262 'step_params': {'maxiter': s['maxiter'] if maxiter is None else maxiter},
263 }
265 if nlevels > 1:
266 description['space_transfer_class'] = mesh_to_mesh_fenics
267 description['space_transfer_params'] = {}
268 description['base_transfer_class'] = base_transfer_mass
269 # prolong_f is not available for the mass formulation: f is a load vector, so it cannot be
270 # interpolated. base_transfer_mass falls back to prolong, re-evaluating f on the fine level.
271 description['base_transfer_params'] = {'finter': False}
273 controller_params = {'logger_level': 30}
275 return description, controller_params, 0.0, nsteps * dt