Peano 4
Loading...
Searching...
No Matches
elastic-wave.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 os
4import sys
5import argparse
6import subprocess
7import peano4
8import exahype2
9
10# Some reasonable default values
11dimensions = 2
12cell_size = 0.1
13amr_levels = 2
14patch_size = 4
15rk_order = 2
16dg_order = 1
17end_time = 3.0
18plot_interval = end_time / 100
19# plot_interval = 0
20rho = 2.7
21p_wave_speed = 6.0
22s_wave_speed = 3.464
23
24scenarios = [
25 "PointExplosion",
26]
27scenario = "PointExplosion"
28
29fv_rusanov_solvers = {
34}
35
36rkdg_rusanov_solvers = {
40}
41
42aderdg_rusanov_solvers = {
44 "ADERDGRusanovGlobalAdaptive": exahype2.solvers.aderdg.rusanov.GlobalAdaptiveTimeStep,
45}
46
47solvers = {}
48solvers.update(fv_rusanov_solvers)
49solvers.update(rkdg_rusanov_solvers)
50solvers.update(aderdg_rusanov_solvers)
51solver = "ADERDGRusanovGlobalAdaptive"
52
53build_modes = {
54 "Release": peano4.output.CompileMode.Release,
55 "Trace": peano4.output.CompileMode.Trace,
56 "Assert": peano4.output.CompileMode.Asserts,
57 "Stats": peano4.output.CompileMode.Stats,
58 "Debug": peano4.output.CompileMode.Debug,
59}
60build_mode = "Release"
61
62storage_types = {
63 "CallStack": exahype2.solvers.fv.Storage.CallStack,
64 "Heap": exahype2.solvers.fv.Storage.Heap,
65 "SmartPointers": exahype2.solvers.fv.Storage.SmartPointers,
66}
67storage_type = "CallStack"
68
69parser = argparse.ArgumentParser(
70 description="ExaHyPE 2 - Elastic Wave Application Script"
71)
72
73parser.add_argument(
74 "-s",
75 "--solver",
76 dest="solver",
77 choices=solvers.keys(),
78 default=solver,
79 help="|".join(solvers.keys()),
80)
81parser.add_argument(
82 "-o",
83 "--output",
84 dest="output_path",
85 type=str,
86 default="solution",
87 help="Output path for project solution output. The project will create a new folder at the given path. Default is 'solution'.",
88)
89parser.add_argument(
90 "-m",
91 "--build-mode",
92 dest="build_mode",
93 choices=build_modes.keys(),
94 default=build_mode,
95 help="|".join(build_modes.keys()),
96)
97parser.add_argument(
98 "-d",
99 "--dim",
100 dest="dim",
101 type=int,
102 default=dimensions,
103 help="Number of space dimensions",
104)
105
106parser.add_argument(
107 "-st",
108 "--storage",
109 dest="storage",
110 choices=storage_types.keys(),
111 default=storage_type,
112 help="|".join(storage_types.keys()),
113)
114parser.add_argument(
115 "-sc",
116 "--scenario",
117 dest="scenario",
118 choices=scenarios,
119 default=scenario,
120 help="|".join(scenarios),
121)
122parser.add_argument(
123 "-et",
124 "--end-time",
125 dest="end_time",
126 type=float,
127 default=end_time,
128 help="End time of simulation",
129)
130parser.add_argument(
131 "-pi",
132 "--plot-interval",
133 dest="plot_interval",
134 type=float,
135 default=plot_interval,
136 help="Time between plots and also between database entries when probes are used.",
137)
138
139parser.add_argument(
140 "-amr",
141 "--amr-levels",
142 dest="amr_levels",
143 type=int,
144 default=amr_levels,
145 help="Number of AMR grid levels on top of max. size of a finite volume cell.",
146)
147parser.add_argument(
148 "-dynamic",
149 "--dynamic-amr",
150 dest="dynamic_amr",
151 action="store_true",
152 default=False,
153 help="Use dynamic AMR.",
154)
155parser.add_argument(
156 "-stateless",
157 "--stateless",
158 dest="use_stateless_pde_terms",
159 action="store_true",
160 default=True,
161 help="Use stateless PDE terms (GPU offloading requires a GPU enabled Peano build and an enclave solver).",
162)
163parser.add_argument(
164 "-f",
165 "--fuse",
166 dest="fused_tasks",
167 type=int,
168 default=16384,
169 help="Number of enclave tasks to fuse into one meta task",
170)
171parser.add_argument(
172 "-vis",
173 "--visualise",
174 dest="visualise",
175 action="store_true",
176 default=False,
177 help="Visualise the output data using the postprocessor.",
178)
179
180parser.add_argument(
181 "-pbc",
182 "--periodic-boundary-conditions",
183 dest="periodic_boundary_conditions",
184 action="store_true",
185 help="Use periodic boundary conditions",
186)
187
188parser.add_argument(
189 "-maxh",
190 "--max-cell-size",
191 dest="max_cell_size",
192 type=float,
193 default=cell_size,
194 help="Maximum size of a single finite volume cell on each axis.",
195)
196parser.add_argument(
197 "-ps",
198 "--patch-size",
199 dest="patch_size",
200 type=int,
201 default=patch_size,
202 help="Number of finite volumes per axis (dimension) per patch.",
203)
204
205parser.add_argument(
206 "-rk-order",
207 "--rk-order",
208 dest="rk_order",
209 type=int,
210 default=rk_order,
211 help="Order of time discretisation for Runge-Kutta scheme.",
212)
213parser.add_argument(
214 "-dg-order",
215 "--dg-order",
216 dest="dg_order",
217 type=int,
218 default=dg_order,
219 help="Order of space discretisation for Discontinuous Galerkin.",
220)
221
222parser.add_argument(
223 "--trees",
224 dest="trees",
225 type=int,
226 default=0,
227 help="Number of trees (subpartitions) per rank.",
228)
229parser.add_argument(
230 "--threads", dest="threads", type=int, default=0, help="Number of threads per rank."
231)
232parser.add_argument(
233 "--timeout", dest="timeout", type=int, default=3600, help="MPI timeout in seconds."
234)
235parser.add_argument(
236 "-fpe",
237 "--fp-exception",
238 dest="enable_fpe",
239 action="store_true",
240 help="Enable a floating-point exception handler",
241)
242
243parser.add_argument(
244 "--no-make",
245 dest="no_make",
246 action="store_true",
247 help="Do not compile the code after generation",
248)
249parser.add_argument(
250 "--executable",
251 dest="executable",
252 type=str,
253 help="Specify the executable name. If none is chosen, the name is automatically deduced by given parameters.",
254)
255args = parser.parse_args()
256
257if args.dim not in [2, 3]:
258 print("Error, dimension must be 2 or 3, you supplied {}".format(args.dim))
259 sys.exit(1)
260
261if args.build_mode not in build_modes:
262 print(
263 "Error, mode must be {} or {}, you supplied {}".format(
264 ", ".join(list(build_modes.keys())[:-1]),
265 list(build_modes.keys())[-1],
266 args.build_mode,
267 )
268 )
269 sys.exit(1)
270
271if args.solver not in solvers:
272 print(
273 "Error, solver must be {} or {}, you supplied {}".format(
274 ", ".join(list(solvers.keys())[:-1]), list(solvers.keys())[-1], args.solver
275 )
276 )
277 sys.exit(1)
278
279if args.storage not in storage_type:
280 print(
281 "Error, storage type must be {} or {}, you supplied {}".format(
282 ", ".join(list(storage_type.keys())[:-1]),
283 list(storage_type.keys())[-1],
284 args.storage,
285 )
286 )
287 sys.exit(1)
288
289if args.scenario not in scenarios:
290 print(
291 "Error, scenario must be {} or {}, you supplied {}".format(
292 ", ".join(scenarios[:-1]), scenarios[-1], args.scenario
293 )
294 )
295 sys.exit(1)
296
297if args.executable:
298 executable_name = args.executable
299else:
300 executable_name = (
301 "elastic-wave-"
302 + str(args.dim)
303 + "d-"
304 + str(args.scenario)
305 + "-"
306 + str(args.solver)
307 )
308 if args.use_stateless_pde_terms:
309 executable_name += "-Stateless"
310 executable_name += "-" + str(args.build_mode)
311
312my_project = exahype2.Project(
313 namespace=["applications", "exahype2", "elasticwave"],
314 project_name="ElasticWave" + args.scenario,
315 directory=".",
316 executable=executable_name,
317)
318
319if args.solver in aderdg_rusanov_solvers:
320 solver_name = "ElasticWave_ADERDG"
321elif args.solver in rkdg_rusanov_solvers:
322 solver_name = "ElasticWave_RKDG"
323elif args.solver in fv_rusanov_solvers:
324 solver_name = "ElasticWave_FV"
325
326unknowns = {"v": args.dim, "sigma": 3}
327
328my_solver = None
329if args.solver in fv_rusanov_solvers and "Adaptive" in args.solver:
330 my_solver = solvers[args.solver](
331 name=solver_name,
332 patch_size=args.patch_size,
333 unknowns=unknowns,
334 auxiliary_variables=0,
335 max_volume_h=args.max_cell_size,
336 min_volume_h=args.max_cell_size * 3.0 ** (-args.amr_levels),
337 time_step_relaxation=0.5,
338 pde_terms_without_state=args.use_stateless_pde_terms
339 )
340 # self._patch.generator = peano4.datamodel.PatchToDoubleArrayOnHeap(self._patch,"double")
341 my_solver.switch_to_heap_storage(
342 cell_data_storage=storage_types[args.storage],
343 face_data_storage=storage_types[args.storage],
344 )
345
346 # from exahype2.solvers.fv.rusanov.kernels import create_compute_Riemann_kernel_for_Rusanov
347 # my_solver._fused_compute_kernel_call_cpu = create_compute_Riemann_kernel_for_Rusanov(
348 # my_solver._flux_implementation,
349 # my_solver._ncp_implementation,
350 # my_solver._source_term_implementation,
351 # compute_max_eigenvalue_of_next_time_step=True,
352 # solver_variant = exahype2.solvers.fv.rusanov.SolverVariant.Multicore,
353 # kernel_variant = exahype2.solvers.fv.rusanov.KernelVariant.BatchedAoSHeap
354 # )
355elif args.solver in fv_rusanov_solvers and "Fixed" in args.solver:
356 min_h = args.max_cell_size * 3.0 ** (-args.amr_levels)
357 admissible_time_step_size = min_h / patch_size * 0.5
358 my_solver = solvers[args.solver](
359 name=solver_name,
360 patch_size=args.patch_size,
361 unknowns=unknowns,
362 auxiliary_variables=0,
363 max_volume_h=args.max_cell_size,
364 min_volume_h=min_h,
365 normalised_time_step_size=admissible_time_step_size,
366 pde_terms_without_state=args.use_stateless_pde_terms
367 )
368elif args.solver in rkdg_rusanov_solvers and "Adaptive" in args.solver:
370 FaceProjections,
371 )
372
373 my_solver = solvers[args.solver](
374 name=solver_name,
375 rk_order=args.rk_order,
376 polynomials=exahype2.solvers.GaussLegendreBasis(args.dg_order),
377 face_projections=FaceProjections.Solution,
378 unknowns=unknowns,
379 auxiliary_variables=0,
380 max_cell_h=args.max_cell_size,
381 min_cell_h=args.max_cell_size * 3.0 ** (-args.amr_levels),
382 time_step_relaxation=0.5,
383 pde_terms_without_state=args.use_stateless_pde_terms
384 )
385elif args.solver in rkdg_rusanov_solvers and "Fixed" in args.solver:
387 FaceProjections,
388 )
389
390 min_h = args.max_cell_size * 3.0 ** (-args.amr_levels)
391 admissible_time_step_size = min_h / patch_size * 0.5
392 my_solver = solvers[args.solver](
393 name=solver_name,
394 rk_order=args.rk_order,
395 polynomials=exahype2.solvers.GaussLegendreBasis(args.dg_order),
396 face_projections=FaceProjections.Solution,
397 unknowns=unknowns,
398 auxiliary_variables=0,
399 max_cell_h=args.max_cell_size,
400 min_cell_h=min_h,
401 time_step_size=admissible_time_step_size,
402 pde_terms_without_state=args.use_stateless_pde_terms
403 )
404elif args.solver in aderdg_rusanov_solvers and "Fixed" in args.solver:
405 # stability condition: CFL for Ader-DG
406 # dt <= c*dx / (lambda_max*(2N+1)
407 # cflAder = [ 1.0, 0.33, 0.17, 0.1, 0.069, 0.045, 0.038, 0.03, 0.02, 0.015]
408 # here lambda_max = 1
409 # dt <= 0.1 * 0.03 / 7 \approx 0.0004
410 # dt <= 0.1 * 0.03 / 11 \approx 0.00025
411 min_h = args.max_cell_size * 3.0 ** (-args.amr_levels)
412 admissible_time_step_size = 0.001
413 my_solver = solvers[args.solver](
414 name=solver_name,
415 order=args.dg_order,
416 unknowns=unknowns,
417 auxiliary_variables=0,
418 max_cell_h=args.max_cell_size,
419 min_cell_h=min_h,
420 time_step_size=admissible_time_step_size,
421 )
422 my_solver.add_kernel_optimizations(
423 is_linear=True, use_kernel_generator=True, architecture="noarch"
424 )
425elif args.solver in aderdg_rusanov_solvers and "Adaptive" in args.solver:
426 my_solver = solvers[args.solver](
427 name=solver_name,
428 order=args.dg_order,
429 unknowns=unknowns,
430 auxiliary_variables=0,
431 max_cell_h=args.max_cell_size,
432 min_cell_h=args.max_cell_size * 3.0 ** (-args.amr_levels),
433 time_step_relaxation=0.5,
434 )
435 my_solver.add_kernel_optimizations(
436 is_linear=True, use_kernel_generator=True, architecture="noarch"
437 )
438
439my_solver.set_implementation(
440 initial_conditions=exahype2.solvers.PDETerms.User_Defined_Implementation,
441 boundary_conditions=exahype2.solvers.PDETerms.User_Defined_Implementation,
442 flux="::applications::exahype2::elasticwave::flux(Q, normal, F);",
443 eigenvalues="return ::applications::exahype2::elasticwave::maxEigenvalue();",
444 refinement_criterion=exahype2.solvers.PDETerms.User_Defined_Implementation,
445 source_term="::applications::exahype2::elasticwave::sourceTerm(Q, volumeCentre, t, S);",
446 ncp=exahype2.solvers.PDETerms.None_Implementation,
447)
448
449my_solver.plot_description = "u, v, sxx, syy, sxy"
450
451my_solver.add_user_solver_includes(
452 """
453#include "ElasticWaveKernels.h"
454"""
455)
456
457my_project.add_solver(my_solver)
458
459if args.periodic_boundary_conditions:
460 periodic_BC = [True, True, True]
461else:
462 periodic_BC = [False, False, False]
463
464my_project.set_global_simulation_parameters(
465 dimensions=args.dim,
466 size=[20.0, 20.0, 20.0][0 : args.dim],
467 offset=[-5.0, 0.0, 0.0][0 : args.dim],
468 min_end_time=args.end_time,
469 max_end_time=args.end_time,
470 first_plot_time_stamp=0.0,
471 time_in_between_plots=args.plot_interval,
472 periodic_BC=periodic_BC,
473)
474
475solution_dir = str(args.output_path) + "-" + str(executable_name)
476my_project.set_output_path(solution_dir)
477if not os.path.exists(solution_dir):
478 os.makedirs(solution_dir)
479
480if args.visualise:
481 output_patch_file = (
482 solution_dir + "/solution-" + str(solver_name) + ".peano-patch-file"
483 )
484 if os.path.isfile(output_patch_file):
485 subprocess.call(
486 [
487 "python3",
488 "../../../python/peano4/visualisation/render.py",
489 str(output_patch_file),
490 ]
491 )
492 sys.exit(0)
493
494if args.trees > 0:
495 my_project.set_load_balancing(
496 "toolbox::loadbalancing::strategies::SpreadOutHierarchically",
497 "new ::exahype2::LoadBalancingConfiguration(0.98, 100, "
498 + str(args.trees)
499 + ")",
500 )
501else:
502 my_project.set_load_balancing(
503 "toolbox::loadbalancing::strategies::SpreadOutHierarchically",
504 "new ::exahype2::LoadBalancingConfiguration(0.98)",
505 )
506if args.threads > 0:
507 my_project.set_number_of_threads(args.threads)
508else:
509 my_project.set_number_of_threads(
510 "tarch::multicore::Core::UseDefaultNumberOfThreads"
511 )
512
513my_project.set_timeout(args.timeout)
514
515my_project.set_multicore_orchestration(
516 "::tarch::multicore::orchestration::Hardcoded::createFuseAll(%s, true, true, ::tarch::accelerator::Device::getInstance().getLocalDeviceId())"
517 % str(args.fused_tasks)
518)
519
520my_project.set_Peano4_installation("../../../", mode=build_modes[args.build_mode])
521my_project = my_project.generate_Peano4_project(verbose=True)
522if args.enable_fpe:
523 my_project.set_fenv_handler("FE_DIVBYZERO | FE_INVALID | FE_OVERFLOW")
524
525my_project.constants.define_value("RHO", str(rho))
526my_project.constants.define_value("P_WAVE_SPEED", str(p_wave_speed))
527my_project.constants.define_value("S_WAVE_SPEED", str(s_wave_speed))
528my_project.constants.export_boolean("UseDynamicAMR", args.dynamic_amr)
529
530my_project.output.makefile.add_cpp_file(solver_name + ".cpp")
531
532my_project.build(
533 make=not args.no_make, make_clean_first=True, throw_away_data_after_build=True
534)
535
536print("Executable is ", executable_name)
537print("Clean object files via 'make clean'")
538print("Recompile the generated code via 'make -j'")
539print("Remove all generated code via 'make distclean'")
540print("Regenerate all code by running 'python3 elastic-wave.py' again")
541print(
542 "Generate (and convert) any postprocessing data by calling 'python3 elastic-wave.py --visualise'"
543)
ExaHyPE 2 project.
Definition Project.py:18
The Gauss-Legendre Basis is by construction the only basis which yields diagonal mass matrices.