Peano
Loading...
Searching...
No Matches
fv-fd-coupling.py
Go to the documentation of this file.
1# This file is part of the ExaHyPE2 project. For conditions of distribution and
2# use, please see the copyright notice at www.peano-framework.org
3import peano4
4import exahype2
5
6parser = exahype2.ArgumentParser(
7 "ExaHyPE 2 - Finite Volumes/Finite Differences Coupling Testing Script"
8)
9
10parser.set_defaults(
11 min_depth=4,
12 degrees_of_freedom=16,
13 end_time=5000,
14 number_of_snapshots=0,
15 periodic_boundary_conditions_x=False,
16 periodic_boundary_conditions_y=False,
17 periodic_boundary_conditions_z=False,
18)
19args = parser.parse_args()
20
21size = [1e6, 1e6, 1e6]
22max_h = 1.1 * min(size) / 3.0**args.min_depth
23min_h = max_h * 3.0 ** (-args.amr_levels)
24
25
28
29euler_unknowns = {
30 "rho": 1,
31 "rhoU": 1,
32 "rhoV": 1,
33}
34if args.dimensions == 3:
35 euler_unknowns["rhoW"] = 1
36euler_unknowns["rhoE"] = 1
37
38euler_initial_conditions = f"""
39 Q[Shortcuts::rho] = 0.1;
40 Q[Shortcuts::rhoU] = 0.0;
41 Q[Shortcuts::rhoV] = 0.0;
42#if DIMENSIONS == 3
43 Q[Shortcuts::rhoW] = 0.0;
44#endif
45 Q[Shortcuts::rhoE] = ((x(0) < {size[0] / 2}) ? (1.01) : (1.0));
46"""
47
48euler_boundary_conditions = """
49 // Reflective boundary conditions
50 Qoutside[Shortcuts::rho] = Qinside[Shortcuts::rho];
51 Qoutside[Shortcuts::rhoU] = -Qinside[Shortcuts::rhoU];
52 Qoutside[Shortcuts::rhoV] = -Qinside[Shortcuts::rhoV];
53#if DIMENSIONS == 3
54 Qoutside[Shortcuts::rhoW] = -Qinside[Shortcuts::rhoW];
55#endif
56 Qoutside[Shortcuts::rhoE] = Qinside[Shortcuts::rhoE];
57"""
58
59compute_primitive_variables = r"""
60 const auto irho = 1.0 / Q[Shortcuts::rho];
61 const auto u0 = Q[Shortcuts::rhoU + 0] * irho;
62 const auto u1 = Q[Shortcuts::rhoU + 1] * irho;
63#if DIMENSIONS == 3
64 const auto u2 = Q[Shortcuts::rhoU + 2] * irho;
65#endif
66
67#if DIMENSIONS == 3
68 const auto uSq = u0 * u0 + u1 * u1 + u2 * u2;
69 const auto u_n = (normal == 0) ? u0 : (normal == 1) ? u1 : u2;
70#else
71 const auto uSq = u0 * u0 + u1 * u1;
72 const auto u_n = (normal == 0) ? u0 : u1;
73#endif
74
75 const auto internalE = Q[Shortcuts::rhoE] - 0.5 * Q[Shortcuts::rho] * uSq;
76 const auto p = (GAMMA - 1.0) * internalE;
77"""
78
79euler_max_eigenvalue = r"""
80 {compute_primitive_variables}
81 const auto speedOfSound = std::sqrt(GAMMA * p * irho);
82 auto result = std::fmax(0.0, std::fabs(u_n - speedOfSound));
83 result = std::fmax(result, std::fabs(u_n + speedOfSound));
84 return result;
85""".format(
86 compute_primitive_variables=compute_primitive_variables
87)
88
89euler_flux = r"""
90 {compute_primitive_variables}
91
92 F[Shortcuts::rho] = Q[Shortcuts::rhoU + normal];
93
94 F[Shortcuts::rhoU + 0] = Q[Shortcuts::rhoU + 0] * u_n;
95 F[Shortcuts::rhoU + 1] = Q[Shortcuts::rhoU + 1] * u_n;
96#if DIMENSIONS == 3
97 F[Shortcuts::rhoU + 2] = Q[Shortcuts::rhoU + 2] * u_n;
98#endif
99
100 F[Shortcuts::rhoU + normal] += p;
101
102 F[Shortcuts::rhoE] = (Q[Shortcuts::rhoE] + p) * u_n;
103""".format(
104 compute_primitive_variables=compute_primitive_variables
105)
106
107euler_solver = exahype2.solvers.fv.godunov.GlobalAdaptiveTimeStep(
108 name="EulerFVSolver",
109 patch_size=args.degrees_of_freedom,
110 unknowns=euler_unknowns,
111 auxiliary_variables=0,
112 min_volume_h=min_h,
113 max_volume_h=max_h,
114 time_step_relaxation=0.5,
115)
116
117euler_solver.set_implementation(
118 initial_conditions=euler_initial_conditions,
119 boundary_conditions=euler_boundary_conditions,
120 flux=euler_flux,
121 max_eigenvalue=euler_max_eigenvalue,
122)
123
124
127
128diffusion_unknowns = {
129 "rhoU": 1,
130 "rhoV": 1,
131}
132if args.dimensions == 3:
133 diffusion_unknowns["rhoW"] = 1
134diffusion_unknowns["rhoE"] = 1
135diffusion_auxiliary_variables = {
136 "rho": 1,
137}
138
139diffusion_initial_conditions = f"""
140 Q[Shortcuts::rhoU] = 0.0;
141 Q[Shortcuts::rhoV] = 0.0;
142#if DIMENSIONS == 3
143 Q[Shortcuts::rhoW] = 0.0;
144#endif
145 Q[Shortcuts::rhoE] = ((x(0) < {size[0] / 2}) ? (1.01) : (1.0));
146 Q[Shortcuts::rho] = 0.1;
147"""
148
149diffusion_boundary_conditions = """
150 // Reflective boundary conditions
151 Qoutside0[Shortcuts::rhoU] = -Qinside[Shortcuts::rhoU];
152 Qoutside0[Shortcuts::rhoV] = -Qinside[Shortcuts::rhoV];
153#if DIMENSIONS == 3
154 Qoutside0[Shortcuts::rhoW] = -Qinside[Shortcuts::rhoW];
155#endif
156 Qoutside0[Shortcuts::rhoE] = Qinside[Shortcuts::rhoE];
157 Qoutside0[Shortcuts::rho] = Qinside[Shortcuts::rho];
158"""
159
160diffusion_stencil = """
161 {{COMPUTE_PRECISION}} velocityLaplacianU = 0.0;
162 {{COMPUTE_PRECISION}} velocityLaplacianV = 0.0;
163#if DIMENSIONS == 3
164 {{COMPUTE_PRECISION}} velocityLaplacianW = 0.0;
165#endif
166
167 const auto irho = 1.0 / Q[Shortcuts::rho];
168 const auto uSquared = (Q[Shortcuts::rhoU] * irho) * (Q[Shortcuts::rhoU] * irho);
169 const auto vSquared = (Q[Shortcuts::rhoV] * irho) * (Q[Shortcuts::rhoV] * irho);
170#if DIMENSIONS == 3
171 const auto wSquared = (Q[Shortcuts::rhoW] * irho) * (Q[Shortcuts::rhoW] * irho);
172 const auto squaredVelocities = uSquared + vSquared + wSquared;
173#else
174 const auto squaredVelocities = uSquared + vSquared;
175#endif
176
177 // x-direction
178 velocityLaplacianU += ((QRight0[Shortcuts::rhoU] / QRight0[Shortcuts::rho]) - 2 * (Q[Shortcuts::rhoU] * irho) + (QLeft0[Shortcuts::rhoU] / QLeft0[Shortcuts::rho])) / (h(0) * h(0));
179 velocityLaplacianV += ((QRight0[Shortcuts::rhoV] / QRight0[Shortcuts::rho]) - 2 * (Q[Shortcuts::rhoV] * irho) + (QLeft0[Shortcuts::rhoV] / QLeft0[Shortcuts::rho])) / (h(0) * h(0));
180#if DIMENSIONS == 3
181 velocityLaplacianW += ((QRight0[Shortcuts::rhoW] / QRight0[Shortcuts::rho]) - 2 * (Q[Shortcuts::rhoW] * irho) + (QLeft0[Shortcuts::rhoW] / QLeft0[Shortcuts::rho])) / (h(0) * h(0));
182#endif
183
184 // y-direction
185 velocityLaplacianU += ((QTop0[Shortcuts::rhoU] / QTop0[Shortcuts::rho]) - 2 * (Q[Shortcuts::rhoU] * irho) + (QDown0[Shortcuts::rhoU] / QDown0[Shortcuts::rho])) / (h(0) * h(0));
186 velocityLaplacianV += ((QTop0[Shortcuts::rhoV] / QTop0[Shortcuts::rho]) - 2 * (Q[Shortcuts::rhoV] * irho) + (QDown0[Shortcuts::rhoV] / QDown0[Shortcuts::rho])) / (h(0) * h(0));
187#if DIMENSIONS == 3
188 velocityLaplacianW += ((QTop0[Shortcuts::rhoW] / QTop0[Shortcuts::rho]) - 2 * (Q[Shortcuts::rhoW] * irho) + (QDown0[Shortcuts::rhoW] / QDown0[Shortcuts::rho])) / (h(0) * h(0));
189#endif
190
191#if DIMENSIONS == 3
192 // z-direction
193 velocityLaplacianU += ((QFront0[Shortcuts::rhoU] / QFront0[Shortcuts::rho]) - 2 * (Q[Shortcuts::rhoU] * irho) + (QBack0[Shortcuts::rhoU] / QBack0[Shortcuts::rho])) / (h(0) * h(0));
194 velocityLaplacianV += ((QFront0[Shortcuts::rhoV] / QFront0[Shortcuts::rho]) - 2 * (Q[Shortcuts::rhoV] * irho) + (QBack0[Shortcuts::rhoV] / QBack0[Shortcuts::rho])) / (h(0) * h(0));
195 velocityLaplacianW += ((QFront0[Shortcuts::rhoW] / QFront0[Shortcuts::rho]) - 2 * (Q[Shortcuts::rhoW] * irho) + (QBack0[Shortcuts::rhoW] / QBack0[Shortcuts::rho])) / (h(0) * h(0));
196#endif
197
198 QNew[Shortcuts::rhoU] = Q[Shortcuts::rhoU] + dt * Q[Shortcuts::rho] * KINEMATIC_VISCOSITY * velocityLaplacianU;
199 QNew[Shortcuts::rhoV] = Q[Shortcuts::rhoV] + dt * Q[Shortcuts::rho] * KINEMATIC_VISCOSITY * velocityLaplacianV;
200#if DIMENSIONS == 3
201 QNew[Shortcuts::rhoW] = Q[Shortcuts::rhoW] + dt * Q[Shortcuts::rho] * KINEMATIC_VISCOSITY * velocityLaplacianW;
202#endif
203
204 auto temperature = [](const double* Q) {
205 const auto irho = 1 / Q[Shortcuts::rho];
206 // Calculate temperature based on energy density
207 const auto uSquared = (Q[Shortcuts::rhoU] * irho) * (Q[Shortcuts::rhoU] * irho);
208 const auto vSquared = (Q[Shortcuts::rhoV] * irho) * (Q[Shortcuts::rhoV] * irho);
209#if DIMENSIONS == 3
210 const auto wSquared = (Q[Shortcuts::rhoW] * irho) * (Q[Shortcuts::rhoW] * irho);
211 const auto squaredVelocities = uSquared + vSquared + wSquared;
212#else
213 const auto squaredVelocities = uSquared + vSquared;
214#endif
215 const auto internalEnergy = Q[Shortcuts::rhoE] - 0.5 * (Q[Shortcuts::rho] * squaredVelocities);
216 const auto temperature = (internalEnergy * (GAMMA - 1.0)) / (Q[Shortcuts::rho] * UNIVERSAL_GAS_CONSTANT);
217 return temperature;
218};
219
220 {{COMPUTE_PRECISION}} temperatureLaplacian = 0.0;
221 const auto temperatureQ = temperature(Q);
222 const auto doubleTemperatureQ = 2 * temperatureQ;
223
224 // x-direction
225 temperatureLaplacian += (temperature(QRight0) - doubleTemperatureQ + temperature(QLeft0)) / (h(0) * h(0));
226 // y-direction
227 temperatureLaplacian += (temperature(QTop0) - doubleTemperatureQ + temperature(QDown0)) / (h(1) * h(1));
228#if DIMENSIONS == 3
229 // z-direction temperatureLaplacian += (temperature(QFront0) - doubleTemperatureQ + temperature(QBack0)) / (h(2) * h(2));
230#endif
231
232 const auto newPressure = Q[Shortcuts::rho] * UNIVERSAL_GAS_CONSTANT * (temperatureQ + dt * THERMAL_DIFFUSIVITY * temperatureLaplacian);
233 QNew[Shortcuts::rhoE] = newPressure / (GAMMA - 1.0) + 0.5 * Q[Shortcuts::rho] * vSquared;
234"""
235
236diffusion_solver = exahype2.solvers.fd.GlobalAdaptiveTimeStep(
237 name="DiffusionSolver",
238 patch_size=args.degrees_of_freedom,
239 unknowns=diffusion_unknowns,
240 auxiliary_variables=diffusion_auxiliary_variables,
241 min_h=min_h,
242 max_h=max_h,
243 time_step_relaxation=0.5,
244 solvers_before_in_coupling=[euler_solver],
245 order=1,
246)
247
248diffusion_solver.set_implementation(
249 initial_conditions=diffusion_initial_conditions,
250 boundary_conditions=diffusion_boundary_conditions,
251 solver=diffusion_stencil,
252 max_eigenvalue="""
253 return 1.0;
254""",
255)
256
257
260
261diffusion_solver._compute_time_step_size = (
262 """
263 const double timeStepSize = repositories::"""
264 + euler_solver.get_name_of_global_instance()
265 + """.getAdmissibleTimeStepSize(false);
266"""
267)
268
269number_of_variables_euler = euler_solver.auxiliary_variables + euler_solver.unknowns
270number_of_variables_diffusion = (
271 diffusion_solver.auxiliary_variables + diffusion_solver.unknowns
272)
273
274# EULER to DIFFUSION
275
276euler_solver.postprocess_updated_patch += f"""// Copy density, energy, and momenta
277 ::toolbox::blockstructured::copyUnknown(
278 {euler_solver.patch_size}, // no of unknowns per dimension
279 fineGridCell{euler_solver.name}Q.value, // source
280 {euler_solver._variable_names.index("rho")}, // density in source
281 {number_of_variables_euler}, // unknowns in source (incl material parameters)
282 0, // no overlap/halo here
283 fineGridCell{diffusion_solver.name}Q.value, // dest
284 {diffusion_solver._variable_names.index("rho")}, // density in dest
285 {number_of_variables_diffusion}, // unknowns in dest (incl material parameters)
286 0 // no overlap/halo here
287 );
288
289 ::toolbox::blockstructured::copyUnknown(
290 {euler_solver.patch_size}, // no of unknowns per dimension
291 fineGridCell{euler_solver.name}Q.value, // source
292 {euler_solver._variable_names.index("rhoU")}, // x-momentum in source
293 {number_of_variables_euler}, // unknowns in source (incl material parameters)
294 0, // no overlap/halo here
295 fineGridCell{diffusion_solver.name}Q.value, // dest
296 {diffusion_solver._variable_names.index("rhoU")}, // x-momentum in dest
297 {number_of_variables_diffusion}, // unknowns in dest (incl material parameters)
298 0 // no overlap/halo here
299 );
300
301 ::toolbox::blockstructured::copyUnknown(
302 {euler_solver.patch_size}, // no of unknowns per dimension
303 fineGridCell{euler_solver.name}Q.value, // source
304 {euler_solver._variable_names.index("rhoV")}, // y-momentum in source
305 {number_of_variables_euler}, // unknowns in source (incl material parameters)
306 0, // no overlap/halo here
307 fineGridCell{diffusion_solver.name}Q.value, // dest
308 {diffusion_solver._variable_names.index("rhoV")}, // y-momentum in dest
309 {number_of_variables_diffusion}, // unknowns in dest (incl material parameters)
310 0 // no overlap/halo here
311 );
312
313 ::toolbox::blockstructured::copyUnknown(
314 {euler_solver.patch_size}, // no of unknowns per dimension
315 fineGridCell{euler_solver.name}Q.value, // source
316 {euler_solver._variable_names.index("rhoE")}, // energy in source
317 {number_of_variables_euler}, // unknowns in source (incl material parameters)
318 0, // no overlap/halo here
319 fineGridCell{diffusion_solver.name}Q.value, // dest
320 {diffusion_solver._variable_names.index("rhoE")}, // energy in dest
321 {number_of_variables_diffusion}, // unknowns in dest (incl material parameters)
322 0 // no overlap/halo here
323 );
324"""
325
326if args.dimensions == 3:
327 euler_solver.postprocess_updated_patch += f"""
328#if DIMENSIONS == 3
329 ::toolbox::blockstructured::copyUnknown(
330 {euler_solver.patch_size}, // no of unknowns per dimension
331 fineGridCell{euler_solver.name}Q.value, // source
332 {euler_solver._variable_names.index("rhoW")}, // z-momentum in source
333 {number_of_variables_euler}, // unknowns in source (incl material parameters)
334 0, // no overlap/halo here
335 fineGridCell{diffusion_solver.name}Q.value, // dest
336 {diffusion_solver._variable_names.index("rhoW")}, // z-momentum in dest
337 {number_of_variables_diffusion}, // unknowns in dest (incl material parameters)
338 0 // no overlap/halo here
339 );
340#endif
341"""
342
343euler_solver.add_user_action_set_includes('#include "toolbox/blockstructured/Copy.h"\n')
344
345# DIFFUSION to EULER
346
347diffusion_solver.postprocess_updated_patch += f"""// Copy energy and momenta
348 ::toolbox::blockstructured::copyUnknown(
349 {diffusion_solver.patch_size}, // no of unknowns per dimension
350 fineGridCell{diffusion_solver.name}Q.value, // source
351 {diffusion_solver._variable_names.index("rhoU")}, // x-momentum in source
352 {number_of_variables_diffusion}, // unknowns in source (incl material parameters)
353 0, // no overlap/halo here
354 fineGridCell{euler_solver.name}Q.value, // dest
355 {euler_solver._variable_names.index("rhoU")}, // x-momentum in dest
356 {number_of_variables_euler}, // unknowns in dest (incl material parameters)
357 0 // no overlap/halo here
358 );
359
360 ::toolbox::blockstructured::copyUnknown(
361 {diffusion_solver.patch_size}, // no of unknowns per dimension
362 fineGridCell{diffusion_solver.name}Q.value, // source
363 {diffusion_solver._variable_names.index("rhoV")}, // y-momentum in source
364 {number_of_variables_diffusion}, // unknowns in source (incl material parameters)
365 0, // no overlap/halo here
366 fineGridCell{euler_solver.name}Q.value, // dest
367 {euler_solver._variable_names.index("rhoV")}, // y-momentum in dest
368 {number_of_variables_euler}, // unknowns in dest (incl material parameters)
369 0 // no overlap/halo here
370 );
371
372 ::toolbox::blockstructured::copyUnknown(
373 {diffusion_solver.patch_size}, // no of unknowns per dimension
374 fineGridCell{diffusion_solver.name}Q.value, // source
375 {diffusion_solver._variable_names.index("rhoE")}, // energy in source
376 {number_of_variables_diffusion}, // unknowns in source (incl material parameters)
377 0, // no overlap/halo here
378 fineGridCell{euler_solver.name}Q.value, // dest
379 {euler_solver._variable_names.index("rhoE")}, // energy in dest
380 {number_of_variables_euler}, // unknowns in dest (incl material parameters)
381 0 // no overlap/halo here
382 );
383"""
384
385if args.dimensions == 3:
386 diffusion_solver.postprocess_updated_patch += f"""
387#if DIMENSIONS == 3
388 ::toolbox::blockstructured::copyUnknown(
389 {diffusion_solver.patch_size}, // no of unknowns per dimension
390 fineGridCell{diffusion_solver.name}Q.value, // source
391 {diffusion_solver._variable_names.index("rhoW")}, // z-momentum in source
392 {number_of_variables_diffusion}, // unknowns in source (incl material parameters)
393 0, // no overlap/halo here
394 fineGridCell{euler_solver.name}Q.value, // dest
395 {euler_solver._variable_names.index("rhoW")}, // z-momentum in dest
396 {number_of_variables_euler}, // unknowns in dest (incl material parameters)
397 0 // no overlap/halo here
398 );
399#endif
400"""
401
402diffusion_solver.add_user_action_set_includes(
403 '#include "toolbox/blockstructured/Copy.h"\n'
404)
405
406euler_solver.last_solver = diffusion_solver
407
408
411
412project = exahype2.Project(
413 namespace=["tests", "exahype2", "fvfd"],
414 project_name=".",
415 directory=".",
416 executable="ExaHyPE",
417)
418
419project.add_solver(euler_solver)
420project.add_solver(diffusion_solver)
421
422if args.number_of_snapshots <= 0:
423 time_in_between_plots = 0.0
424else:
425 time_in_between_plots = args.end_time / args.number_of_snapshots
426 project.set_output_path(args.output)
427
428project.set_global_simulation_parameters(
429 dimensions=args.dimensions,
430 size=size[0 : args.dimensions],
431 offset=[0.0, 0.0, 0.0][0 : args.dimensions],
432 min_end_time=args.end_time,
433 max_end_time=args.end_time,
434 first_plot_time_stamp=0.0,
435 time_in_between_plots=time_in_between_plots,
436 periodic_BC=[
437 args.periodic_boundary_conditions_x,
438 args.periodic_boundary_conditions_y,
439 args.periodic_boundary_conditions_z,
440 ],
441)
442
443project.set_load_balancer("new ::exahype2::LoadBalancingConfiguration")
444project.set_Peano4_installation(
445 "../../../", mode=peano4.output.string_to_mode(args.build_mode)
446)
447project = project.generate_Peano4_project(verbose=False)
448project.output.makefile.add_CXX_flag("-DGAMMA=1.4")
449project.output.makefile.add_CXX_flag("-DUNIVERSAL_GAS_CONSTANT=8.31")
450project.output.makefile.add_CXX_flag("-DKINEMATIC_VISCOSITY=100")
451project.output.makefile.add_CXX_flag("-DTHERMAL_DIFFUSIVITY=150")
452project.set_fenv_handler(args.fpe)
453project.build(make=True, make_clean_first=True, throw_away_data_after_build=True)