7 import peano4.toolbox.particles
9 sys.path.insert(0, os.path.abspath(
"../../../../applications/exahype2/ccz4"))
11 from Probe_file_gene
import tracer_seeds_generate
12 from ComputeFirstDerivatives
import ComputeFirstDerivativesFD4RK
15 "release": peano4.output.CompileMode.Release,
16 "trace": peano4.output.CompileMode.Trace,
17 "assert": peano4.output.CompileMode.Asserts,
18 "stats": peano4.output.CompileMode.Stats,
19 "debug": peano4.output.CompileMode.Debug,
66 if __name__ ==
"__main__":
67 parser = argparse.ArgumentParser(description=
"ExaHyPE 2 - CCZ4 script")
74 help=
"upper limit for refinement. Refers to volume size, i.e. not to patch size",
82 help=
"lower limit for refinement (set to 0 to make it equal to max_h - default). Refers to volume size, i.e. not to patch size",
90 help=
"Patch size, i.e. number of volumes per patch per direction",
95 dest=
"plot_step_size",
98 help=
"Plot step size (0 to switch it off)",
101 "-m",
"--mode", dest=
"mode", default=
"release", help=
"|".join(modes.keys())
108 help=
"Run with accelerator support",
113 dest=
"implementation",
114 default=
"fd4-rk1-adaptive",
115 choices=[
"fd4-rk1-adaptive",
"fd4-rk1-adaptive-enclave"],
116 help=
"Pick solver type",
120 "--no-periodic-boundary-conditions",
122 action=
"store_false",
124 help=
"switch on or off the periodic BC",
128 "--sommerfeld-boundary-conditions",
129 dest=
"sommerfeld_bc",
132 help=
"switch on or off the Sommerfeld radiative BC",
138 choices=[
"none",
"poisson"],
140 help=
"switch on or off the AMR boundary marker",
148 help=
"End (terminal) time",
153 dest=
"plot_start_time",
156 help=
"start time for plot",
162 default=
"single-puncture",
174 "-cfl",
"--CFL-ratio", dest=
"cfl", type=float, default=0.1, help=
"Set CFL ratio"
182 help=
"Add tracers and specify the seeds. 0-switch off, 1-static point tracer, 2-moving point tracer",
190 help=
"name of output tracer file (temporary)",
198 help=
"name of output executable file",
202 "--output-directory",
206 help=
"specify the output directory, include the patch file and tracer file",
211 dest=
"interpolation",
215 "linear-constant-extrap",
216 "linear-linear-extrap",
217 "linear-con-extrap-lin-normal-interp",
218 "linear-lin-extrap-lin-normal-interp",
223 default=
"linear-lin-extrap-lin-normal-interp",
224 help=
"interpolation scheme for AMR",
230 choices=[
"average",
"inject",
"matrix",
"second_order"],
232 help=
"restriction scheme for AMR",
240 help=
"enable double communication per timestep, used in the soccz4 formulation.",
243 for k, v
in floatparams.items():
246 dest=
"CCZ4{}".format(k),
249 help=
"default: %(default)s",
251 for k, v
in intparams.items():
255 dest=
"CCZ4{}".format(k),
258 help=
"default: %(default)s, choose refinement criterion, 0-no refinement, 1-radius based, 2-SBH phi gradient based, 3-BBH phi gradient based. Notice: 2 and 3 only work with -ext Full",
263 dest=
"CCZ4{}".format(k),
266 help=
"default: %(default)s",
269 args = parser.parse_args()
273 print(args.implementation)
275 args.implementation ==
"fd4-rk1-adaptive"
276 or args.implementation ==
"fd4-rk4-adaptive"
278 SuperClass = exahype2.solvers.rkfd.fd4.GlobalAdaptiveTimeStep
279 if args.implementation ==
"fd4-rk1-adaptive-enclave":
280 SuperClass = exahype2.solvers.rkfd.fd4.GlobalAdaptiveTimeStepWithEnclaveTasking
282 if args.implementation ==
"fd4-rk1-fixed" or args.implementation ==
"fd4-rk4-fixed":
283 SuperClass = exahype2.solvers.rkfd.fd4.GlobalFixedTimeStep
287 self, name, patch_size, min_volume_h, max_volume_h, cfl, domain_r, KOSig
310 number_of_unknowns = sum(unknowns.values())
313 #include "../CCZ4Kernels.h"
315 if SuperClass == exahype2.solvers.rkfd.fd4.GlobalAdaptiveTimeStep:
317 args.implementation ==
"fd4-rk1-adaptive"
318 or args.implementation ==
"fd4-rk1-adaptive-enclave"
326 patch_size=patch_size,
328 unknowns=number_of_unknowns,
329 auxiliary_variables=0,
330 min_meshcell_h=min_volume_h,
331 max_meshcell_h=max_volume_h,
332 time_step_relaxation=cfl,
334 reconstruction_with_rk=args.SO_flag,
338 == exahype2.solvers.rkfd.fd4.GlobalAdaptiveTimeStepWithEnclaveTasking
341 args.implementation ==
"fd4-rk1-adaptive"
342 or args.implementation ==
"fd4-rk1-adaptive-enclave"
350 patch_size=patch_size,
352 unknowns=number_of_unknowns,
353 auxiliary_variables=0,
354 min_meshcell_h=min_volume_h,
355 max_meshcell_h=max_volume_h,
356 time_step_relaxation=cfl,
358 pde_terms_without_state=
True,
361 SuperClass == exahype2.solvers.rkfd.fd4.GlobalFixedTimeStep
363 == exahype2.solvers.rkfd.fd4.GlobalFixedTimeStepWithEnclaveTasking
365 if args.implementation ==
"fd4-rk1-fixed":
372 patch_size=patch_size,
374 unknowns=number_of_unknowns,
375 auxiliary_variables=0,
376 min_meshcell_h=min_volume_h,
377 max_meshcell_h=max_volume_h,
378 normalised_time_step_size=0.01,
385 patch_size=patch_size,
386 unknowns=number_of_unknowns,
387 auxiliary_variables=0,
388 min_volume_h=min_volume_h,
389 max_volume_h=max_volume_h,
390 time_step_relaxation=cfl
398 self.set_implementation(
399 boundary_conditions=exahype2.solvers.PDETerms.User_Defined_Implementation,
400 ncp=exahype2.solvers.PDETerms.User_Defined_Implementation,
401 flux=exahype2.solvers.PDETerms.None_Implementation,
402 source_term=exahype2.solvers.PDETerms.User_Defined_Implementation,
403 refinement_criterion=exahype2.solvers.PDETerms.User_Defined_Implementation,
404 eigenvalues=exahype2.solvers.PDETerms.User_Defined_Implementation,
409 == exahype2.solvers.rkfd.fd4.GlobalAdaptiveTimeStepWithEnclaveTasking
412 self._flux_implementation,
413 self._ncp_implementation,
414 self._source_term_implementation,
415 compute_max_eigenvalue_of_next_time_step=
True,
416 solver_variant=exahype2.solvers.rkfd.kernels.SolverVariant.Multicore,
417 kernel_variant=exahype2.solvers.rkfd.kernels.KernelVariant.BatchedAoSHeap,
418 KOSigma=self._KO_Sigma,
426 SuperClass == exahype2.solvers.rkfd.fd4.GlobalAdaptiveTimeStep
427 or SuperClass == exahype2.solvers.rkfd.fd4.GlobalFixedTimeStep
429 == exahype2.solvers.rkfd.fd4.GlobalAdaptiveTimeStepWithEnclaveTasking
434 constexpr int itmax = {{NUMBER_OF_GRID_CELLS_PER_PATCH_PER_AXIS}} * {{NUMBER_OF_GRID_CELLS_PER_PATCH_PER_AXIS}};
438 constexpr int itmax = {{NUMBER_OF_GRID_CELLS_PER_PATCH_PER_AXIS}} * {{NUMBER_OF_GRID_CELLS_PER_PATCH_PER_AXIS}} * {{NUMBER_OF_GRID_CELLS_PER_PATCH_PER_AXIS}};
448 constexpr int itmax = {{NUMBER_OF_VOLUMES_PER_AXIS}} * {{NUMBER_OF_VOLUMES_PER_AXIS}};
452 constexpr int itmax = {{NUMBER_OF_VOLUMES_PER_AXIS}} * {{NUMBER_OF_VOLUMES_PER_AXIS}} * {{NUMBER_OF_VOLUMES_PER_AXIS}};
458 for (int i=0;i<itmax;i++)
460 applications::exahype2::ccz4::enforceCCZ4constraints( newQ+index );
461 index += {{NUMBER_OF_UNKNOWNS}} + {{NUMBER_OF_AUXILIARY_VARIABLES}};
467 SuperClass.create_action_sets(self)
468 self._action_set_couple_resolution_transitions_and_handle_dynamic_mesh_refinement.additional_includes +=
"""
469 #include "../CCZ4Kernels.h"
474 We take this routine to add a few additional include statements.
477 SuperClass.get_user_action_set_includes(self) + self.
_my_user_includes_my_user_includes
486 if args.exe_name !=
"":
489 if not args.tra_name ==
"de":
490 exe +=
"_" + args.tra_name
491 project = exahype2.Project(
492 [
"benchmarks",
"exahype2",
"ccz4"],
"ccz4", executable=exe
504 == exahype2.solvers.aderdg.NonFusedGenericRusanovFixedTimeStepSize
509 time_step_size = 0.001
510 except Exception
as e:
517 if SuperClass == exahype2.solvers.rkdg.rusanov.GlobalAdaptiveTimeStep:
519 except Exception
as e:
521 msg =
"Warning: RKDG not supported on this machine"
523 userinfo.append((msg, e))
526 solver_name =
"ADERDG" + solver_name
528 solver_name =
"RKDG" + solver_name
530 solver_name = solver_name
537 print(
"No minimal mesh size chosen. Set it to max mesh size (regular grid)")
541 my_solver = exahype2.solvers.aderdg.NonFusedGenericRusanovFixedTimeStepSize(
546 exahype2.solvers.aderdg.Polynomials.Gauss_Legendre,
551 ncp=exahype2.solvers.PDETerms.User_Defined_Implementation,
552 sources=exahype2.solvers.PDETerms.User_Defined_Implementation,
564 userinfo.append((
"CFL ratio set as " +
str(args.cfl),
None))
566 userinfo.append((
"The solver is " + args.implementation,
None))
572 if args.interpolation ==
"order-2":
573 my_solver.overlap = 2
575 if args.interpolation ==
"constant":
576 my_solver.interpolation =
"piecewise_constant"
577 print(
"Interpolation rule: piecewise_constant")
578 if args.interpolation ==
"linear-constant-extrap":
579 my_solver.interpolation =
"linear_with_constant_extrapolation"
580 print(
"Interpolation rule: linear constant extrapolation")
581 if args.interpolation ==
"linear-linear-extrap":
582 my_solver.interpolation =
"linear_with_linear_extrapolation"
583 print(
"Interpolation rule: linear extrapolation")
584 if args.interpolation ==
"linear-con-extrap-lin-normal-interp":
585 my_solver.interpolation = (
586 "linear_with_constant_extrapolation_and_linear_normal_interpolation"
589 "Interpolation rule: linear+constant extrapolation and linear normal interpolation"
592 if args.interpolation ==
"order-2":
593 my_solver.interpolation =
"linear"
595 tem_interp = [
"TP_constant",
"TP_linear_with_linear_extrap_normal_interp"]
596 tem_restrict = [
"TP_inject_normal_extrap",
"TP_average_normal_extrap"]
597 if "fd4-rk" in args.implementation:
598 if args.interpolation ==
"matrix":
599 exahype2.solvers.rkfd.fd4.switch_to_FD4_matrix_interpolation(my_solver)
600 elif args.interpolation ==
"second_order":
601 exahype2.solvers.rkfd.fd4.switch_to_FD4_second_order_interpolation(
604 elif args.interpolation ==
"third_order":
605 exahype2.solvers.rkfd.fd4.switch_to_FD4_third_order_interpolation(my_solver)
607 exahype2.solvers.rkfd.fd4.switch_to_FD4_tensor_product_interpolation(
608 my_solver, tem_interp[1]
610 userinfo.append((
"FD4 Interpolation: " + tem_interp[1],
None))
612 if args.restriction ==
"matrix":
613 exahype2.solvers.rkfd.fd4.switch_to_FD4_matrix_restriction(my_solver)
614 elif args.restriction ==
"second_order":
615 exahype2.solvers.rkfd.fd4.switch_to_FD4_second_order_restriction(my_solver)
619 exahype2.solvers.rkfd.fd4.switch_to_FD4_tensor_product_restriction(
620 my_solver, tem_restrict[1]
622 userinfo.append((
"FD4 Restriction: " + tem_restrict[1],
None))
629 for k, v
in intparams.items():
630 intparams.update({k: eval(
"args.CCZ4{}".format(k))})
631 for k, v
in floatparams.items():
632 floatparams.update({k: eval(
"args.CCZ4{}".format(k))})
634 if args.SO_flag ==
True:
635 intparams.update({
"SO": 1})
637 if args.scenario ==
"two-punctures":
638 msg =
"Periodic BC deactivated because you pick Puncture scenario\nInitialize binary black holes"
640 periodic_boundary_conditions = [
False,
False,
False]
645 userinfo.append((msg,
None))
646 elif args.scenario ==
"single-puncture":
647 msg =
"Periodic BC deactivated because you pick Puncture scenario\nInitialize single black hole"
649 periodic_boundary_conditions = [
False,
False,
False]
650 intparams.update({
"swi": 0})
651 userinfo.append((msg,
None))
652 elif args.periodic_bc ==
True:
653 msg =
"Periodic BC set"
655 periodic_boundary_conditions = [
True,
True,
True]
656 userinfo.append((msg,
None))
658 msg =
"WARNING: Periodic BC deactivated by hand"
660 periodic_boundary_conditions = [
False,
False,
False]
661 userinfo.append((msg,
None))
665 if args.scenario ==
"gauge":
666 solverconstants +=
"static constexpr int Scenario=0; /* Gauge wave */ \n "
667 userinfo.append((
"picking gauge wave scenario",
None))
668 floatparams.update({
"sk": 0.0})
669 floatparams.update({
"bs": 0.0})
670 intparams.update({
"LapseType": 0})
671 elif (args.scenario ==
"two-punctures")
or (args.scenario ==
"single-puncture"):
672 solverconstants +=
"static constexpr int Scenario=2; /* Two-puncture */ \n"
673 userinfo.append((
"picking black hole scenario",
None))
675 raise Exception(
"Scenario " + args.scenario +
" is now unknown")
677 for k, v
in floatparams.items():
678 solverconstants +=
"static constexpr double {} = {};\n".format(
679 "CCZ4{}".format(k), v
681 for k, v
in intparams.items():
682 solverconstants +=
"static constexpr int {} = {};\n".format(
683 "CCZ4{}".format(k), v
685 my_solver.add_solver_constants(solverconstants)
687 project.add_solver(my_solver)
689 build_mode = modes[args.mode]
696 floatparams.update({
"domain_r": args.CCZ4domain_r})
697 dr = floatparams[
"domain_r"]
698 offset = [-dr, -dr, -dr]
699 domain_size = [2 * dr, 2 * dr, 2 * dr]
700 msg =
"domain set as " +
str(offset) +
" and " +
str(domain_size)
702 userinfo.append((msg,
None))
704 project.set_global_simulation_parameters(
709 args.plot_start_time,
711 periodic_boundary_conditions,
717 +
str(args.plot_start_time)
718 +
", plot step size: "
719 +
str(args.plot_step_size),
723 userinfo.append((
"Terminal time: " +
str(args.end_time),
None))
725 project.set_Peano4_installation(
"../../../..", build_mode)
731 if not args.path ==
"./":
735 project.set_output_path(path)
736 probe_point = [-12, -12, 0.0]
737 project.add_plot_filter(probe_point, [24.0, 24.0, 0.01], 1)
739 project.set_load_balancer(
"new ::exahype2::LoadBalancingConfiguration")
744 peano4_project = project.generate_Peano4_project(verbose=
True)
745 peano4_project.output.makefile.add_header_search_path(
746 "../../../../applications/exahype2/ccz4/"
750 args.scenario ==
"gauge"
751 or args.scenario ==
"linear"
752 or args.scenario ==
"dia_gauge"
753 or args.scenario ==
"flat"
756 elif args.scenario ==
"two-punctures" or args.scenario ==
"single-puncture":
757 peano4_project.output.makefile.add_linker_flag(
"-lm -lgsl -lgslcblas")
758 peano4_project.output.makefile.add_cpp_file(
759 "../../../../applications/exahype2/ccz4/libtwopunctures/TP_Utilities.cpp"
761 peano4_project.output.makefile.add_cpp_file(
762 "../../../../applications/exahype2/ccz4/libtwopunctures/TP_Parameters.cpp"
764 peano4_project.output.makefile.add_cpp_file(
765 "../../../../applications/exahype2/ccz4/libtwopunctures/TP_Logging.cpp"
767 peano4_project.output.makefile.add_cpp_file(
768 "../../../../applications/exahype2/ccz4/libtwopunctures/TwoPunctures.cpp"
770 peano4_project.output.makefile.add_cpp_file(
771 "../../../../applications/exahype2/ccz4/libtwopunctures/CoordTransf.cpp"
773 peano4_project.output.makefile.add_cpp_file(
774 "../../../../applications/exahype2/ccz4/libtwopunctures/Equations.cpp"
776 peano4_project.output.makefile.add_cpp_file(
777 "../../../../applications/exahype2/ccz4/libtwopunctures/FuncAndJacobian.cpp"
779 peano4_project.output.makefile.add_cpp_file(
780 "../../../../applications/exahype2/ccz4/libtwopunctures/Newton.cpp"
782 peano4_project.output.makefile.add_CXX_flag(
"-DIncludeTwoPunctures")
784 raise Exception(
"Scenario " + args.scenario +
" is now unknown")
786 peano4_project.output.makefile.add_CXX_flag(
"-DCCZ4EINSTEIN")
788 peano4_project.output.makefile.add_h_file(
"CCZ4.h")
789 peano4_project.output.makefile.add_cpp_file(
"CCZ4.cpp")
790 peano4_project.output.makefile.add_cpp_file(
791 "../../../../applications/exahype2/ccz4/InitialValues.cpp"
793 peano4_project.output.makefile.add_cpp_file(
794 "../../../../applications/exahype2/ccz4/CCZ4Kernels.cpp"
797 peano4_project.generate(throw_away_data_after_generation=
False)
798 peano4_project.build(make_clean_first=
True)
801 userinfo.append((
"the executable file name: " + exe,
None))
802 userinfo.append((
"output directory: " + path,
None))
803 print(
"=========================================================")
804 if not args.add_tracer == 0:
805 userinfo.append((
"tracer output file: " + args.path +
"Tracer-test",
None))
806 if len(userinfo) > 0:
807 print(
"The building information:")
808 for msg, exception
in userinfo:
809 if exception
is None:
812 print(msg,
"Exception: {}".format(exception))
815 print(
"=========================================================")