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