Coverage for pySDC/projects/StroemungsRaum/run_Navier_Stokes_TaylorGreen_FEniCS.py: 81%
69 statements
« prev ^ index » next coverage.py v7.16.1, created at 2026-09-25 20:28 +0000
1import dolfin as df
2import numpy as np
4from pySDC.implementations.controller_classes.controller_nonMPI import controller_nonMPI
5from pySDC.projects.StroemungsRaum.problem_classes.NavierStokes_2D_TaylorGreen_monolithic_FEniCS import (
6 fenics_NSE_2D_TaylorGreen,
7)
8from pySDC.projects.StroemungsRaum.sweepers.generic_implicit_mass import (
9 generic_implicit_mass,
10 generic_implicit_mass_diffbc,
11)
14def setup(
15 t0=0.0,
16 dt=0.1,
17 periodic=False,
18 differentiated_bc=False,
19 nelems=24,
20 nu=0.1,
21 num_nodes=4,
22 maxiter=40,
23 restol=1e-12,
24 Sol_tol=None,
25):
26 """
27 Helper routine to set up parameters
29 Args:
30 t0: float,
31 initial time
32 dt: float,
33 time step size
34 periodic: bool,
35 use periodic instead of time-dependent Dirichlet conditions in x
36 differentiated_bc: bool,
37 impose the time-dependent boundary data in differentiated form, which recovers
38 the order it otherwise costs; requires periodic=False
39 nelems: int,
40 number of elements per spatial direction
41 nu: float,
42 kinematic viscosity
43 num_nodes: int,
44 number of collocation nodes
45 maxiter: int,
46 maximum number of SDC iterations
47 restol: float,
48 residual tolerance
49 Sol_tol: float,
50 absolute tolerance of the Newton solve at each node; defaults to one decade below
51 ``restol``, which SDC cannot converge past
53 Returns:
54 description: dict,
55 pySDC description dictionary containing problem and method parameters.
56 controller_params: dict,
57 Parameters for the pySDC controller.
58 """
59 # initialize level parameters
60 level_params = dict()
61 level_params['restol'] = restol
62 level_params['dt'] = dt
64 # initialize step parameters
65 step_params = dict()
66 step_params['maxiter'] = maxiter
68 # initialize sweeper parameters
69 sweeper_params = dict()
70 sweeper_params['quad_type'] = 'RADAU-RIGHT'
71 sweeper_params['num_nodes'] = num_nodes
72 sweeper_params['QI'] = 'LU'
74 # initialize problem parameters
75 problem_params = dict()
76 problem_params['nelems'] = nelems
77 problem_params['t0'] = t0
78 problem_params['order'] = 2
79 problem_params['nu'] = nu
80 problem_params['periodic'] = periodic
81 problem_params['differentiated_bc'] = differentiated_bc
82 problem_params['Sol_tol'] = restol / 10 if Sol_tol is None else Sol_tol
84 # initialize controller parameters
85 controller_params = dict()
86 controller_params['logger_level'] = 30
88 # Fill description dictionary
89 description = dict()
90 description['problem_class'] = fenics_NSE_2D_TaylorGreen
91 description['sweeper_class'] = generic_implicit_mass_diffbc if differentiated_bc else generic_implicit_mass
92 description['problem_params'] = problem_params
93 description['sweeper_params'] = sweeper_params
94 description['level_params'] = level_params
95 description['step_params'] = step_params
97 return description, controller_params
100def run_simulation(description, controller_params, Tend):
101 """
102 Run the time integration for the 2D Taylor-Green Navier-Stokes benchmark.
104 Args:
105 description: dict,
106 pySDC problem and method description.
107 controller_params: dict,
108 Parameters for the pySDC controller.
109 Tend: float,
110 Final simulation time.
112 Returns:
113 P: problem instance,
114 Problem instance holding the function spaces and the exact solution.
115 stats: dict,
116 Collected runtime statistics.
117 uend: dtype_u,
118 Final solution at time Tend.
119 """
120 t0 = description['problem_params']['t0']
122 controller = controller_nonMPI(num_procs=1, controller_params=controller_params, description=description)
124 P = controller.MS[0].levels[0].prob
125 uend, stats = controller.run(u0=P.u_exact(t0), t0=t0, Tend=Tend)
127 return P, stats, uend
130def relative_errors(u, uref):
131 """
132 Relative L2 errors in velocity and pressure between two solutions on the same space.
134 Args:
135 u: dtype_u,
136 Numerical solution.
137 uref: dtype_u,
138 Reference solution.
140 Returns:
141 tuple of float: relative L2 error in velocity and in pressure.
142 """
143 un, pn = u.values.split(deepcopy=True)
144 ur, pr = uref.values.split(deepcopy=True)
146 return (
147 df.errornorm(ur, un, 'L2', degree_rise=0) / df.norm(ur, 'L2'),
148 df.errornorm(pr, pn, 'L2', degree_rise=0) / df.norm(pr, 'L2'),
149 )
152def run_postprocessing(P, uend, Tend):
153 """
154 Compute relative L2 errors between the numerical and the exact solution at the final time.
156 Args:
157 P: problem instance,
158 Problem instance holding the exact solution.
159 uend: dtype_u,
160 Final solution at time Tend.
161 Tend: float,
162 Final simulation time.
164 Returns:
165 tuple of float: relative L2 error in velocity and in pressure.
166 """
167 return relative_errors(uend, P.u_exact(Tend))
170def order_study(dts, Tend, periodic=False, **kwargs):
171 r"""
172 Measure the observed temporal order of convergence.
174 Errors are *not* taken against the exact solution: the spatial discretization error
175 dominates it for any affordable mesh, which hides the temporal order completely. Instead
176 consecutive step sizes are compared with each other (Richardson), which cancels the spatial
177 error exactly because every run uses the same mesh and needs no reference run, at the cost
178 of one order estimate.
180 Args:
181 dts: list of float,
182 Step sizes to run, largest first, each one half of the previous.
183 Tend: float,
184 Final simulation time; must be an integer multiple of every step size.
185 periodic: bool,
186 Use periodic instead of time-dependent Dirichlet conditions in x.
187 kwargs:
188 Passed on to :func:`setup`.
190 Returns:
191 dts_out: list of float,
192 Step sizes the errors belong to; one shorter than ``dts``.
193 errors_u: list of float,
194 Relative L2 velocity error per step size.
195 errors_p: list of float,
196 Relative L2 pressure error per step size.
197 """
198 solutions = []
199 for dt in dts:
200 description, controller_params = setup(dt=dt, periodic=periodic, **kwargs)
201 solutions.append(run_simulation(description, controller_params, Tend)[2])
203 pairs = zip(solutions[:-1], solutions[1:], strict=True)
204 errors = [relative_errors(u, ref) for u, ref in pairs]
206 return dts[:-1], [e[0] for e in errors], [e[1] for e in errors]
209def observed_order(dts, errors):
210 """
211 Observed order of convergence between consecutive step sizes.
213 Args:
214 dts: list of float,
215 Step sizes.
216 errors: list of float,
217 Corresponding errors.
219 Returns:
220 list of float: observed orders, one shorter than the inputs.
221 """
222 return [np.log(errors[i] / errors[i + 1]) / np.log(dts[i] / dts[i + 1]) for i in range(len(dts) - 1)]
225def main():
226 r"""
227 Run the order study for both boundary condition variants and report the observed orders.
229 RADAU-RIGHT with M nodes has design order :math:`2M-1` and, on a stiff problem with
230 time-dependent boundary data, drops to the stiff order :math:`M+1`. The gap the benchmark
231 can show is therefore :math:`M-2`, and **nothing at all is visible for M = 2**, where the
232 two coincide at 3. M = 4 is used here because it is the cheapest setting that makes the
233 reduction unmistakable: order 7 against 5 in the pressure.
234 """
235 Tend = 0.2
236 dts = [0.2, 0.1, 0.05, 0.025]
238 cases = [
239 ('periodic', dict(periodic=True)),
240 ('time-dependent Dirichlet', dict(periodic=False)),
241 ('time-dependent Dirichlet, differentiated', dict(periodic=False, differentiated_bc=True)),
242 ]
244 for label, kwargs in cases:
245 dts_out, errors_u, errors_p = order_study(dts, Tend, **kwargs)
247 print(f'\n{label} boundary conditions in x:')
248 print(f'{"dt":>10} {"err(u)":>12} {"order(u)":>9} {"err(p)":>12} {"order(p)":>9}')
249 orders_u = [None] + observed_order(dts_out, errors_u)
250 orders_p = [None] + observed_order(dts_out, errors_p)
251 for dt, eu, ou, ep, op in zip(dts_out, errors_u, orders_u, errors_p, orders_p, strict=True):
252 su = ' --- ' if ou is None else f'{ou:9.2f}'
253 sp = ' --- ' if op is None else f'{op:9.2f}'
254 print(f'{dt:10.5f} {eu:12.4e} {su} {ep:12.4e} {sp}')
257if __name__ == "__main__":
258 main()