5from abc
import abstractmethod
11 Abstract base class for any CCZ4 solver
13 Each CCZ4 solver inherits from this abstract base class which really only
14 defines some generic stuff such as the unknowns and includes that every
15 single solver will use.
17 The solver should, more or less, work out of the box, but you have to do
18 three things if you use a subclass:
20 1. If you use a CCZ4 solver, you will still have to add all the libraries to
21 your Peano project such that the Makefile picks them up. For this, the
22 solver offers an add_makefile_parameters().
24 2. You have to set the initial conditions via
26 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
27 my_solver.set_implementation(initial_conditions = " " "
28 for (int i=0; i<NumberOfUnknowns+NumberOfAuxiliaryVariables; i++) Q[i] = 0.0;
29 ::applications::exahype2::ccz4::gaugeWave(Q, volumeCentre, 0);
31 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
33 At this point, different CCZ4 solver variants might require different
34 syntax. The term volumeCentre for example above is only defined in a
37 3. Finally, you have to add domain-specific constants to the project.
38 For this, call add_all_solver_constants(). See the comment below.
40 Further to that, you might want to have to set boundary conditions. By
41 default, we do not set any boundary conditions. This works fine if
42 periodic boundary conditions are used. But once you switch off periodic
43 boundary conditions, you have to tell the solver how to treat the boundary.
44 This is typically done via set_implementation(), too.
46 ## More complex scenarios
48 Setting particular implementations via set_implementation() is not always
49 convenient or possible. You might want to add new functions to your classes,
50 do something in the solver constructor, and so forth. If so, feel free to
51 modify the file MySolverName.cpp which the tool generates. In this context,
52 you might want to pass in
54 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
55 my_solver.set_implementation(initial_conditions = exahype2.solvers.PDETerms.User_Defined_Implementation,
56 refinement_criterion = exahype2.solvers.PDETerms.User_Defined_Implementation,
57 boundary_conditions=exahype2.solvers.PDETerms.User_Defined_Implementation
60 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
62 which ensures that you get the right hook-in methods generated when you
63 invoke the Python script for the first time. These methods will contain
64 todo comments for you. Subsequent runs of the Python API should not
65 overwrite the solver implementation.
69 Each CCZ4 solver requires a minimal set of constants. These are represented
70 by integer_constants and double_constants. Please augment these dictionaries.
71 Eventually, you have to submit all the constants via add_all_solver_constants().
77 Dictionary which specifies the unknown names plus their cardinality
79 Has to be class attribute, as we need it in the constructor, i.e. before the
80 abstract object is created.
83 _FO_formulation_unknowns = {
106 Primary unknowns of the CCZ4 formulation which are there in the initial
107 formulation. All the other variables are auxiliary variables, i.e. ones
108 introduced to return to a first-order formulation. Unfortunately, the
109 ordering in _FO_formulation_unknows is motivated by the original papers
110 and not by the fact which quantities are original ones and which one are
111 helper or auxiliary variables.
114 _SO_formulation_unknowns = {
126 Default_Time_Step_Size_Relaxation = 0.1
133 Initialise the two dictionaries with default values (which work).
167 Add the headers for the compute kernels and initial condition implementations
169 Usually called by the subclass constructor.
172 self.add_user_action_set_includes(
174#include "CCZ4Kernels.h"
175#include "SecondOrderAuxiliaryVariablesReconstruction.h"
178 self.add_user_solver_includes(
180#include "CCZ4Kernels.h"
181#include "InitialValues.h"
182#include "SecondOrderAuxiliaryVariablesReconstruction.h"
190 Add domain-specific constants
192 I need a couple of constants. I could either replace them directly
193 within the Python snippets below, but I prefer here to go a different
194 way and to export them as proper C++ constants.
196 There are two ways to inject solver constants into Peano: We can either
197 add them to the Makefile as global const expressions, or we can add
198 them to the ExaHyPE2 solver. The latter is the route we go down here,
199 as these constants logically belong to the solver and not to the project.
201 This operation uses the parent class' add_solver_constants(). You still
202 can use this operation to add further parameters. Or you can, as a user,
203 always add new entries to integer_constants or double_constants and then
204 call this routine rather than adding individual constants one by one.
208 self.add_solver_constants(
209 "static constexpr int {} = {};".format(key, value)
212 self.add_solver_constants(
213 "static constexpr double {} = {};".format(key, value)
219 Add include path and minimal required cpp files to makefile
221 If you have multiple CCZ4 solvers, i.e. different solvers of CCZ4 or multiple
222 instances of the CCZ4 type, please call this operation only once on one of
223 your solvers. At the moment, I add hte following cpp files to the setup:
227 - SecondOrderAuxiliaryVariablesReconstruction.cpp
229 You can always add further files via
230 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
231 peano4_project.output.makefile.add_cpp_file( "mypath/myfile.cpp" )
232 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
235 if path_of_ccz4_application[-1] !=
"/":
236 path_of_ccz4_application +=
"/"
238 peano4_project.output.makefile.add_cpp_file(
239 path_of_ccz4_application +
"InitialValues.cpp"
241 peano4_project.output.makefile.add_cpp_file(
242 path_of_ccz4_application +
"CCZ4Kernels.cpp"
244 peano4_project.output.makefile.add_cpp_file(
245 path_of_ccz4_application +
"SecondOrderAuxiliaryVariablesReconstruction.cpp"
247 peano4_project.output.makefile.add_header_search_path(path_of_ccz4_application)
255 number_of_entries_between_two_db_flushes,
256 data_delta_between_two_snapsots,
257 time_delta_between_two_snapsots,
258 clear_database_after_flush,
259 tracer_unknowns=None,
263 Add tracer to project
265 Consult exahype2.tracer.DumpTracerIntoDatabase for an explanation of
266 some of the arguments. Most of them are simply piped through to this
269 The tracer is given a name and initial coordinates (list of three-tuples).
270 We need to know the underlying project as well, as we have to add the
271 tracing to the time stepping and the database update to the plotting.
272 ~~~~~~~~~~~~~~~~~~~~~~~
273 project.add_action_set_to_timestepping(my_interpolation)
274 project.add_action_set_to_timestepping(exahype2.tracer.DumpTracerIntoDatabase(
275 particle_set=tracer_particles,
277 filename=name + "-" + self._name,
278 number_of_entries_between_two_db_flushes=number_of_entries_between_two_db_flushes,
280 data_delta_between_two_snapsots = data_delta_between_two_snapsots,
281 time_delta_between_two_snapsots = time_delta_between_two_snapsots,
282 clear_database_after_flush = True,
284 ~~~~~~~~~~~~~~~~~~~~~~~
287 assert "should not be called"
296 CCZ4 solver using finite volumes and global adaptive time stepping incl enclave tasking
298 Please consult CCZ4Solver_FV_GlobalAdaptiveTimeStepWithEnclaveTasking.
303 self, name, patch_size, min_volume_h, max_volume_h, pde_terms_without_state
305 AbstractCCZ4Solver.__init__(self)
306 exahype2.solvers.fv.rusanov.GlobalAdaptiveTimeStep.__init__(
309 patch_size=patch_size,
311 auxiliary_variables=0,
312 min_volume_h=min_volume_h,
313 max_volume_h=max_volume_h,
314 time_step_relaxation=AbstractCCZ4Solver.Default_Time_Step_Size_Relaxation,
319 boundary_conditions=exahype2.solvers.PDETerms.Empty_Implementation,
321 flux=exahype2.solvers.PDETerms.None_Implementation,
323 refinement_criterion=exahype2.solvers.PDETerms.Empty_Implementation,
334 number_of_entries_between_two_db_flushes,
335 data_delta_between_two_snapsots,
336 time_delta_between_two_snapsots,
337 clear_database_after_flush,
342 Add tracer to project
344 Consult exahype2.tracer.DumpTracerIntoDatabase for an explanation of
345 some of the arguments. Most of them are simply piped through to this
348 project: exahype2.Project
356 number_of_entries_between_two_db_flushes,
357 data_delta_between_two_snapsots,
358 time_delta_between_two_snapsots,
359 clear_database_after_flush,
370 CCZ4 solver using finite volumes and global adaptive time stepping incl enclave tasking
372 The constructor of this classs is straightforward and realises the standard
373 steps of any numerical implementation of the CCZ4 scheme:
375 1. Init the actual numerical scheme. This happens through the constructor
378 2. Add the header files that we need, i.e. those files which contain the
379 actual CCZ4 implementation.
381 3. Add some constants that any CCZ4 C++ code requires.
383 4. Set the actual implementation, i.e. link the generic PDE terms to the
384 CCZ4-specific function calls.
386 5. Add the CCZ4-specific postprocessing.
396 pde_terms_without_state,
400 Construct solver with enclave tasking and adaptive time stepping
403 AbstractCCZ4Solver.__init__(self)
404 exahype2.solvers.fv.rusanov.GlobalAdaptiveTimeStepWithEnclaveTasking.__init__(
407 patch_size=patch_size,
409 auxiliary_variables=0,
410 min_volume_h=min_volume_h,
411 max_volume_h=max_volume_h,
412 time_step_relaxation=AbstractCCZ4Solver.Default_Time_Step_Size_Relaxation,
413 pde_terms_without_state=pde_terms_without_state,
418 boundary_conditions=exahype2.solvers.PDETerms.Empty_Implementation,
420 flux=exahype2.solvers.PDETerms.None_Implementation,
422 refinement_criterion=exahype2.solvers.PDETerms.Empty_Implementation,
433 number_of_entries_between_two_db_flushes,
434 data_delta_between_two_snapsots,
435 time_delta_between_two_snapsots,
436 clear_database_after_flush,
441 Add tracer to project
443 Consult exahype2.tracer.DumpTracerIntoDatabase for an explanation of
444 some of the arguments. Most of them are simply piped through to this
448 project: exahype2.Project
456 number_of_entries_between_two_db_flushes,
457 data_delta_between_two_snapsots,
458 time_delta_between_two_snapsots,
459 clear_database_after_flush,
466 double deltaQSerialised[NumberOfUnknowns*3];
467 for (int i=0; i<NumberOfUnknowns; i++) {
468 deltaQSerialised[i+0*NumberOfUnknowns] = 0.0;
469 deltaQSerialised[i+1*NumberOfUnknowns] = 0.0;
470 deltaQSerialised[i+2*NumberOfUnknowns] = 0.0;
472 deltaQSerialised[i+normal*NumberOfUnknowns] = deltaQ[i];
474 ::applications::exahype2::ccz4::ncp(BTimesDeltaQ, Q, deltaQSerialised, normal%Dimensions, CCZ4LapseType, CCZ4ds, CCZ4c, CCZ4e, CCZ4f, CCZ4bs, CCZ4sk, CCZ4xi, CCZ4mu, CCZ4SO);
480 double deltaQSerialised[NumberOfUnknowns*3];
481 for (int i=0; i<NumberOfUnknowns; i++) {
482 deltaQSerialised[i+0*NumberOfUnknowns] = 0.0;
483 deltaQSerialised[i+1*NumberOfUnknowns] = 0.0;
484 deltaQSerialised[i+2*NumberOfUnknowns] = 0.0;
486 deltaQSerialised[i+normal*NumberOfUnknowns] = deltaQ[i];
488 ::applications::exahype2::ccz4::ncp(BTimesDeltaQ, Q, deltaQSerialised, normal%Dimensions, CCZ4LapseType, CCZ4ds, CCZ4c, CCZ4e, CCZ4f, CCZ4bs, CCZ4sk, CCZ4xi, CCZ4mu, CCZ4SO);
494 tarch::memset(S, 0.0, NumberOfUnknowns*sizeof(double));
495 ::applications::exahype2::ccz4::source(S,Q, CCZ4LapseType, CCZ4ds, CCZ4c, CCZ4e, CCZ4f, CCZ4bs, CCZ4sk, CCZ4xi, CCZ4itau, CCZ4eta, CCZ4k1, CCZ4k2, CCZ4k3, CCZ4SO);
501 tarch::memset(S, 0.0, NumberOfUnknowns*sizeof(double));
502 ::applications::exahype2::ccz4::source(S,Q, CCZ4LapseType, CCZ4ds, CCZ4c, CCZ4e, CCZ4f, CCZ4bs, CCZ4sk, CCZ4xi, CCZ4itau, CCZ4eta, CCZ4k1, CCZ4k2, CCZ4k3, CCZ4SO);
508 return ::applications::exahype2::ccz4::maxEigenvalue(Q, normal%Dimensions, CCZ4e, CCZ4ds, CCZ4GLMc, CCZ4GLMd );
514 return ::applications::exahype2::ccz4::maxEigenvalue(Q, normal%Dimensions, CCZ4e, CCZ4ds, CCZ4GLMc, CCZ4GLMd );
521 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}};
523 for (int i=0;i<itmax;i++)
525 applications::exahype2::ccz4::enforceCCZ4constraints( newQ+index );
526 index += {{NUMBER_OF_UNKNOWNS}} + {{NUMBER_OF_AUXILIARY_VARIABLES}};
535 constexpr int itmax = {{NUMBER_OF_VOLUMES_PER_AXIS}} * {{NUMBER_OF_VOLUMES_PER_AXIS}} * {{NUMBER_OF_VOLUMES_PER_AXIS}};
537 for (int i=0;i<itmax;i++)
539 applications::exahype2::ccz4::enforceCCZ4constraints( newQ+index );
540 index += {{NUMBER_OF_UNKNOWNS}} + {{NUMBER_OF_AUXILIARY_VARIABLES}};
551 number_of_entries_between_two_db_flushes,
552 data_delta_between_two_snapsots,
553 time_delta_between_two_snapsots,
554 clear_database_after_flush,
559 Add tracer to project
561 Consult exahype2.tracer.DumpTracerIntoDatabase for an explanation of
562 some of the arguments. Most of them are simply piped through to this
565 I realise this as a separate routine, as we need it for all FV flavours
568 number_of_attributes = (
569 (solver.unknowns + solver.auxiliary_variables)
570 if tracer_unknowns ==
None
571 else len(tracer_unknowns)
573 tracer_particles = project.add_tracer(
574 name=name, attribute_count=number_of_attributes
577 particle_set=tracer_particles, coordinates=coordinates
579 init_action_set.descend_invocation_order = 0
580 project.add_action_set_to_initialisation(init_action_set)
582 project_on_tracer_properties_kernel =
""
583 if tracer_unknowns ==
None:
584 project_on_tracer_properties_kernel = (
585 "::exahype2::fv::projectAllValuesOntoParticle_piecewiseLinear"
588 elif len(tracer_unknowns) == 1:
589 project_on_tracer_properties_kernel = (
590 "::exahype2::fv::projectValueOntoParticle_piecewiseLinear<{},{}>".format(
591 i, tracer_unknowns.index(i)
595 project_on_tracer_properties_kernel = (
596 "::exahype2::fv::projectValuesOntoParticle_piecewiseLinear<{}>".format(
606 project_on_tracer_properties_kernel=project_on_tracer_properties_kernel,
608 tracing_action_set.descend_invocation_order = (
609 solver._action_set_update_cell.descend_invocation_order + 1
611 project.add_action_set_to_timestepping(tracing_action_set)
612 project.add_action_set_to_initialisation(tracing_action_set)
615 particle_set=tracer_particles,
617 filename=name +
"-" + solver._name,
618 number_of_entries_between_two_db_flushes=number_of_entries_between_two_db_flushes,
620 data_delta_between_two_snapsots=data_delta_between_two_snapsots,
621 time_delta_between_two_snapsots=time_delta_between_two_snapsots,
622 clear_database_after_flush=clear_database_after_flush,
624 dump_into_database_action_set.descend_invocation_order = (
625 solver._action_set_update_cell.descend_invocation_order + 2
627 project.add_action_set_to_timestepping(dump_into_database_action_set)
635 number_of_entries_between_two_db_flushes,
636 data_delta_between_two_snapsots,
637 time_delta_between_two_snapsots,
638 clear_database_after_flush,
643 I realise this as a separate routine, as we need it for all FD4 flavours
645 This is a wrapper around all the tracer handling. It adds the tracer to the
646 exahype2.Project, but it also instantiates the solution to tracer mapping
647 as well as the database bookkeeping.
649 @param tracer_unknowns: Integer
650 You can set this variable to None. In this case, all variables are
654 number_of_attributes = (
655 (solver.unknowns + solver.auxiliary_variables)
656 if tracer_unknowns ==
None
657 else len(tracer_unknowns)
659 tracer_particles = project.add_tracer(
660 name=name, attribute_count=number_of_attributes
662 project.add_action_set_to_initialisation(
664 particle_set=tracer_particles, coordinates=coordinates
667 project_on_tracer_properties_kernel =
""
668 if tracer_unknowns ==
None:
669 project_on_tracer_properties_kernel = (
670 "::exahype2::fv::projectAllValuesOntoParticle_piecewiseLinear"
672 elif len(tracer_unknowns) == 1:
673 project_on_tracer_properties_kernel = (
674 "::exahype2::fv::projectValueOntoParticle_piecewiseLinear<{},{}>".format(
675 i, tracer_unknowns.index(i)
679 project_on_tracer_properties_kernel = (
680 "::exahype2::fv::projectValuesOntoParticle_piecewiseLinear<{}>".format(
690 project_on_tracer_properties_kernel=project_on_tracer_properties_kernel,
692 tracing_action_set.descend_invocation_order = (
693 solver._action_set_compute_final_linear_combination.descend_invocation_order + 1
695 project.add_action_set_to_timestepping(tracing_action_set)
696 project.add_action_set_to_initialisation(tracing_action_set)
699 particle_set=tracer_particles,
701 filename=name +
"-" + solver._name,
702 number_of_entries_between_two_db_flushes=number_of_entries_between_two_db_flushes,
704 data_delta_between_two_snapsots=data_delta_between_two_snapsots,
705 time_delta_between_two_snapsots=time_delta_between_two_snapsots,
706 clear_database_after_flush=clear_database_after_flush,
708 dump_into_database_action_set.descend_invocation_order = (
709 solver._action_set_compute_final_linear_combination.descend_invocation_order + 2
711 project.add_action_set_to_timestepping(dump_into_database_action_set)
720 CCZ4 solver using fourth-order finite differences and global adaptive time stepping incl enclave tasking
722 The constructor of this classs is straightforward and realises the standard
723 steps of any numerical implementation of the CCZ4 scheme:
725 1. Init the actual numerical scheme. This happens through the constructor
728 2. Add the header files that we need, i.e. those files which contain the
729 actual CCZ4 implementation.
731 3. Add some constants that any CCZ4 C++ code requires.
733 4. Set the actual implementation, i.e. link the generic PDE terms to the
734 CCZ4-specific function calls.
736 5. Add the CCZ4-specific postprocessing.
738 6. Switch to higher-order interpolation and restriction.
749 pde_terms_without_state,
756 Calibrate the default time step size calibration with 1/16 to take into
757 account that we have a higher-order numerical scheme.
760 AbstractCCZ4Solver.__init__(self)
762 AbstractCCZ4Solver.enable_second_order(self)
763 exahype2.solvers.rkfd.fd4.GlobalAdaptiveTimeStepWithEnclaveTasking.__init__(
766 patch_size=patch_size,
769 auxiliary_variables=0,
770 min_meshcell_h=min_meshcell_h,
771 max_meshcell_h=max_meshcell_h,
773 pde_terms_without_state=pde_terms_without_state,
779 boundary_conditions=exahype2.solvers.PDETerms.Empty_Implementation,
781 flux=exahype2.solvers.PDETerms.None_Implementation,
783 refinement_criterion=exahype2.solvers.PDETerms.Empty_Implementation,
790 # Use second order interpolation and restriction
791 exahype2.solvers.rkfd.fd4.switch_to_FD4_second_order_interpolation(
794 exahype2.solvers.rkfd.fd4.switch_to_FD4_second_order_restriction(
800 exahype2.solvers.rkfd.fd4.switch_to_FD4_third_order_interpolation(self)
801 exahype2.solvers.rkfd.fd4.switch_to_FD4_third_order_restriction(self)
804 # Use matrix interpolation and restriction
805 exahype2.solvers.rkfd.fd4.switch_to_FD4_matrix_interpolation(
808 exahype2.solvers.rkfd.fd4.switch_to_FD4_matrix_restriction(
814 # Use tensor product interpolation and restriction
815 exahype2.solvers.rkfd.fd4.switch_to_FD4_tensor_product_interpolation(
817 "TP_linear_with_linear_extrap_normal_interp"
819 exahype2.solvers.rkfd.fd4.switch_to_FD4_tensor_product_restriction(
821 "TP_average_normal_extrap"
830 number_of_entries_between_two_db_flushes,
831 data_delta_between_two_snapsots,
832 time_delta_between_two_snapsots,
833 clear_database_after_flush,
838 Add tracer to project
840 This is a delegate to add_tracer_to_FD4_solver() which passes the
841 object in as first argument.
843 Consult exahype2.tracer.DumpTracerIntoDatabase for an explanation of
844 some of the arguments. Most of them are simply piped through to this
847 @param project: exahype2.Project
849 @param tracer_unknowns: Integer
850 You can set this variable to None. In this case, all variables are
859 number_of_entries_between_two_db_flushes,
860 data_delta_between_two_snapsots,
861 time_delta_between_two_snapsots,
862 clear_database_after_flush,
872 CCZ4 solver using fourth-order finite differences and global adaptive time stepping without enclave tasking
874 Consult CCZ4Solver_FD4_GlobalAdaptiveTimeStepWithEnclaveTasking please.
891 Calibrate the default time step size calibration with 1/16 to take into
892 account that we have a higher-order numerical scheme.
895 AbstractCCZ4Solver.__init__(self)
897 AbstractCCZ4Solver.enable_second_order(self)
898 exahype2.solvers.rkfd.fd4.GlobalAdaptiveTimeStep.__init__(
901 patch_size=patch_size,
904 auxiliary_variables=0,
905 min_meshcell_h=min_meshcell_h,
906 max_meshcell_h=max_meshcell_h,
913 boundary_conditions=exahype2.solvers.PDETerms.Empty_Implementation,
915 flux=exahype2.solvers.PDETerms.None_Implementation,
917 refinement_criterion=exahype2.solvers.PDETerms.Empty_Implementation,
924 # Use second order interpolation and restriction
925 exahype2.solvers.rkfd.fd4.switch_to_FD4_second_order_interpolation(
928 exahype2.solvers.rkfd.fd4.switch_to_FD4_second_order_restriction(
934 exahype2.solvers.rkfd.fd4.switch_to_FD4_third_order_interpolation(self)
935 exahype2.solvers.rkfd.fd4.switch_to_FD4_third_order_restriction(self)
938 # Use matrix interpolation and restriction
939 exahype2.solvers.rkfd.fd4.switch_to_FD4_matrix_interpolation(
942 exahype2.solvers.rkfd.fd4.switch_to_FD4_matrix_restriction(
948 # Use tensor product interpolation and restriction
949 exahype2.solvers.rkfd.fd4.switch_to_FD4_tensor_product_interpolation(
951 "TP_linear_with_linear_extrap_normal_interp"
953 exahype2.solvers.rkfd.fd4.switch_to_FD4_tensor_product_restriction(
955 "TP_average_normal_extrap"
964 number_of_entries_between_two_db_flushes,
965 data_delta_between_two_snapsots,
966 time_delta_between_two_snapsots,
967 clear_database_after_flush,
972 Add tracer to project
974 Consult exahype2.tracer.DumpTracerIntoDatabase for an explanation of
975 some of the arguments. Most of them are simply piped through to this
978 project: exahype2.Project
986 number_of_entries_between_two_db_flushes,
987 data_delta_between_two_snapsots,
988 time_delta_between_two_snapsots,
989 clear_database_after_flush=clear_database_after_flush,
990 tracer_unknowns=tracer_unknowns,
1000 Variation of classic FD4 which relies on second order PDE formulation
1002 The traditional ExaHyPE CCZ4 formulation is the first order formulation
1003 introduced by Dumbser et al. In this formulation, the second order terms
1004 in CCZ4 are substituted with helper variables which represent first order
1005 derivatives. While formally straightforward, keeping the whole system
1006 consistent and stricly hyperbolic is a different challenge.
1008 In this revised version, we have to evolve the primary quantities of CCZ4
1009 and also the helper variables, which blows the overall system up to 59
1010 equations in its simplest form. The work by Dumbser and others suggest that
1011 this is a consistent and stable approach, but limited work is actually
1012 published on proper physical simulations. We therefore also implemented a
1013 second order PDE version within ExaHyPE.
1015 This second order variant is not really second order from the start.
1016 Instead, we use the first order formulation, and we reconstruct the helper
1017 term via finite differences prior to the compute kernel application. That is,
1018 the compute kernels see variables representing first order derivatives, and
1019 they also evolve these guys. Afterwards, we throw away the evolved quantities
1020 and reconstruct them from the primary unknowns prior to the next time step.
1022 This might not be super efficient (it would be faster to stick to the
1023 second order formulation right from the start), but it allows us to reuse
1024 the compute kernels written for the first order PDE formulation.
1028 We have now a smaller number of real unknowns, i.e. only those guys who
1029 belong to the "original" second-order formulation. The remaining quantities
1030 compared to a first-order formulation are technically material or auxiliary
1031 quantities. We model them as such, which allows ExaHyPE`s data management
1032 to deal more efficiently with them.
1035 reconstruction_type: "4thOrder", "centralDifferences", "leftDifference", "rightDifference"
1046 reconstruction_type,
1052 Calibrate the default time step size calibration with 1/16 to take into
1053 account that we have a higher-order numerical scheme.
1056 AbstractCCZ4Solver.__init__(self)
1057 exahype2.solvers.rkfd.fd4.GlobalAdaptiveTimeStepWithEnclaveTasking.__init__(
1060 patch_size=patch_size,
1065 min_meshcell_h=min_meshcell_h,
1066 max_meshcell_h=max_meshcell_h,
1067 time_step_relaxation=AbstractCCZ4Solver.Default_Time_Step_Size_Relaxation
1074 boundary_conditions=exahype2.solvers.PDETerms.Empty_Implementation,
1076 double deltaQSerialised[NumberOfUnknowns*3];
1077 for (int i=0; i<NumberOfUnknowns; i++) {
1078 deltaQSerialised[i+0*NumberOfUnknowns] = 0.0;
1079 deltaQSerialised[i+1*NumberOfUnknowns] = 0.0;
1080 deltaQSerialised[i+2*NumberOfUnknowns] = 0.0;
1082 deltaQSerialised[i+normal*NumberOfUnknowns] = deltaQ[i];
1084 ::applications::exahype2::ccz4::ncpSecondOrderFormulation(BTimesDeltaQ, Q, deltaQSerialised, normal%Dimensions, CCZ4LapseType, CCZ4ds, CCZ4c, CCZ4e, CCZ4f, CCZ4bs, CCZ4sk, CCZ4xi, CCZ4mu, CCZ4SO);
1086 flux=exahype2.solvers.PDETerms.None_Implementation,
1088 tarch::memset(S, 0.0, NumberOfUnknowns*sizeof(double));
1089 ::applications::exahype2::ccz4::sourceSecondOrderFormulation(S,Q, CCZ4LapseType, CCZ4ds, CCZ4c, CCZ4e, CCZ4f, CCZ4bs, CCZ4sk, CCZ4xi, CCZ4itau, CCZ4eta, CCZ4k1, CCZ4k2, CCZ4k3);
1091 refinement_criterion=exahype2.solvers.PDETerms.Empty_Implementation,
1094 return ::applications::exahype2::ccz4::maxEigenvalueSecondOrderFormulation(Q, normal%Dimensions, CCZ4e, CCZ4ds, CCZ4GLMc, CCZ4GLMd );
1100 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}};
1102 for (int i=0;i<itmax;i++)
1104 applications::exahype2::ccz4::enforceCCZ4constraintsSecondOrderFormulation( newQ+index );
1105 index += {{NUMBER_OF_UNKNOWNS}} + {{NUMBER_OF_AUXILIARY_VARIABLES}};
1111 # Use second order interpolation and restriction
1112 exahype2.solvers.rkfd.fd4.switch_to_FD4_second_order_interpolation(
1115 exahype2.solvers.rkfd.fd4.switch_to_FD4_second_order_restriction(
1121 exahype2.solvers.rkfd.fd4.switch_to_FD4_third_order_interpolation(self)
1122 exahype2.solvers.rkfd.fd4.switch_to_FD4_third_order_restriction(self)
1125 # Use matrix interpolation and restriction
1126 exahype2.solvers.rkfd.fd4.switch_to_FD4_matrix_interpolation(
1129 exahype2.solvers.rkfd.fd4.switch_to_FD4_matrix_restriction(
1135 # Use tensor product interpolation and restriction
1136 exahype2.solvers.rkfd.fd4.switch_to_FD4_tensor_product_interpolation(
1138 "TP_linear_with_linear_extrap_normal_interp"
1140 exahype2.solvers.rkfd.fd4.switch_to_FD4_tensor_product_restriction(
1142 "TP_average_normal_extrap"
1148::exahype2::CellData reconstructedPatchData(
1154 nullptr // targetPatch
1156::applications::exahype2::ccz4::recomputeAuxiliaryVariablesFD4_"""
1157 + reconstruction_type
1159 reconstructedPatchData,
1160 {{NUMBER_OF_GRID_CELLS_PER_PATCH_PER_AXIS}},
1162 {{NUMBER_OF_UNKNOWNS}},
1163 {{NUMBER_OF_AUXILIARY_VARIABLES}}
1173 number_of_entries_between_two_db_flushes,
1174 data_delta_between_two_snapsots,
1175 time_delta_between_two_snapsots,
1176 clear_database_after_flush,
1181 Add tracer to project
1183 Consult exahype2.tracer.DumpTracerIntoDatabase for an explanation of
1184 some of the arguments. Most of them are simply piped through to this
1187 project: exahype2.Project
1190 number_of_attributes = (
1192 if tracer_unknowns ==
None
1193 else len(tracer_unknowns)
1195 tracer_particles = project.add_tracer(
1196 name=name, attribute_count=number_of_attributes
1199 particle_set=tracer_particles, coordinates=coordinates
1201 init_action_set.descend_invocation_order = 0
1202 project.add_action_set_to_initialisation(init_action_set)
1204 project_on_tracer_properties_kernel =
""
1205 if tracer_unknowns ==
None:
1206 project_on_tracer_properties_kernel = (
1207 "::exahype2::fv::projectAllValuesOntoParticle_piecewiseLinear"
1209 elif len(tracer_unknowns) == 1:
1210 project_on_tracer_properties_kernel =
"::exahype2::fv::projectValueOntoParticle_piecewiseLinear<{},{}>".format(
1211 i, tracer_unknowns.index(i)
1214 project_on_tracer_properties_kernel = (
1215 "::exahype2::fv::projectValuesOntoParticle_piecewiseLinear<{}>".format(
1225 project_on_tracer_properties_kernel=project_on_tracer_properties_kernel,
1227 tracing_action_set.descend_invocation_order = (
1231 project.add_action_set_to_timestepping(tracing_action_set)
1232 project.add_action_set_to_initialisation(tracing_action_set)
1235 particle_set=tracer_particles,
1238 number_of_entries_between_two_db_flushes=number_of_entries_between_two_db_flushes,
1239 output_precision=10,
1240 data_delta_between_two_snapsots=data_delta_between_two_snapsots,
1241 time_delta_between_two_snapsots=time_delta_between_two_snapsots,
1242 clear_database_after_flush=
True,
1244 dump_into_database_action_set.descend_invocation_order = (
1248 project.add_action_set_to_timestepping(dump_into_database_action_set)
1253 double dQdxSerialised[NumberOfUnknowns*3];
1254 for (int i=0; i<NumberOfUnknowns; i++) {
1255 dQdxSerialised[i+0*NumberOfUnknowns] = 0.0;
1256 dQdxSerialised[i+1*NumberOfUnknowns] = 0.0;
1257 dQdxSerialised[i+2*NumberOfUnknowns] = 0.0;
1259 dQdxSerialised[i+normal*NumberOfUnknowns] = deltaQ[i];
1261 ::applications::exahype2::ccz4::ncp(BTimesDeltaQ, Q, dQdxSerialised, normal%Dimensions, CCZ4LapseType, CCZ4ds, CCZ4c, CCZ4e, CCZ4f, CCZ4bs, CCZ4sk, CCZ4xi, CCZ4mu, CCZ4SO);
1267 tarch::memset(S, 0.0, NumberOfUnknowns*sizeof(double));
1268 ::applications::exahype2::ccz4::source(S,Q, CCZ4LapseType, CCZ4ds, CCZ4c, CCZ4e, CCZ4f, CCZ4bs, CCZ4sk, CCZ4xi, CCZ4itau, CCZ4eta, CCZ4k1, CCZ4k2, CCZ4k3, CCZ4SO);
1274 return ::applications::exahype2::ccz4::maxEigenvalue(Q, normal%Dimensions, CCZ4e, CCZ4ds, CCZ4GLMc, CCZ4GLMd );
1281 constexpr int itmax = ({{DG_ORDER}}+1) * ({{DG_ORDER}}+1) * ({{DG_ORDER}}+1);
1283 for (int i=0;i<itmax;i++)
1285 applications::exahype2::ccz4::enforceCCZ4constraints( newQ+index );
1286 index += {{NUMBER_OF_UNKNOWNS}} + {{NUMBER_OF_AUXILIARY_VARIABLES}};
1297 number_of_entries_between_two_db_flushes,
1298 data_delta_between_two_snapsots,
1299 time_delta_between_two_snapsots,
1300 clear_database_after_flush,
1303 number_of_attributes = (
1304 (self.unknowns + self.auxiliary_variables)
1305 if tracer_unknowns ==
None
1306 else len(tracer_unknowns)
1308 tracer_particles = project.add_tracer(
1309 name=name, attribute_count=number_of_attributes
1312 particle_set=tracer_particles, coordinates=coordinates
1314 init_action_set.descend_invocation_order = 0
1315 project.add_action_set_to_initialisation(init_action_set)
1317 assert tracer_unknowns ==
None
1322 project_on_tracer_properties_kernel=
"::exahype2::dg::projectAllValuesOntoParticle",
1324 tracing_action_set.descend_invocation_order = (
1325 self._action_set_compute_final_linear_combination_and_project_solution_onto_faces.descend_invocation_order
1328 project.add_action_set_to_timestepping(tracing_action_set)
1329 project.add_action_set_to_initialisation(tracing_action_set)
1332 particle_set=tracer_particles,
1334 filename=name +
"-" + self._name,
1335 number_of_entries_between_two_db_flushes=number_of_entries_between_two_db_flushes,
1336 output_precision=10,
1337 data_delta_between_two_snapsots=data_delta_between_two_snapsots,
1338 time_delta_between_two_snapsots=time_delta_between_two_snapsots,
1339 clear_database_after_flush=clear_database_after_flush,
1341 dump_into_database_action_set.descend_invocation_order = (
1342 self._action_set_compute_final_linear_combination_and_project_solution_onto_faces.descend_invocation_order
1345 project.add_action_set_to_timestepping(dump_into_database_action_set)
1354 CCZ4 solver using Runge-Kutta Discontinuous Galerkin and global adaptive time stepping incl enclave tasking
1356 The constructor of this classs is straightforward and realises the standard
1357 steps of any numerical implementation of the CCZ4 scheme:
1359 1. Init the actual numerical scheme. This happens through the constructor
1362 2. Add the header files that we need, i.e. those files which contain the
1363 actual CCZ4 implementation.
1365 3. Add some constants that any CCZ4 C++ code requires.
1367 4. Set the actual implementation, i.e. link the generic PDE terms to the
1368 CCZ4-specific function calls.
1370 5. Add the CCZ4-specific postprocessing.
1372 6. Switch to higher-order interpolation and restriction.
1383 pde_terms_without_state,
1387 Construct solver with enclave tasking
1390 AbstractCCZ4Solver.__init__(self)
1391 exahype2.solvers.rkdg.rusanov.GlobalAdaptiveTimeStepWithEnclaveTasking.__init__(
1395 polynomials=polynomials,
1397 auxiliary_variables=0,
1398 min_cell_h=min_cell_h,
1399 max_cell_h=max_cell_h,
1400 time_step_relaxation=AbstractCCZ4Solver.Default_Time_Step_Size_Relaxation,
1401 pde_terms_without_state=pde_terms_without_state,
1407 boundary_conditions=exahype2.solvers.PDETerms.Empty_Implementation,
1409 flux=exahype2.solvers.PDETerms.None_Implementation,
1411 refinement_criterion=exahype2.solvers.PDETerms.Empty_Implementation,
1424 number_of_entries_between_two_db_flushes,
1425 data_delta_between_two_snapsots,
1426 time_delta_between_two_snapsots,
1427 clear_database_after_flush,
1432 Add tracer to project
1434 Consult exahype2.tracer.DumpTracerIntoDatabase for an explanation of
1435 some of the arguments. Most of them are simply piped through to this
1438 At this point, we have not yet created the Peano 4 project. Therefore, we
1439 have not yet befilled the time stepping action set.
1441 project: exahype2.Project
1449 number_of_entries_between_two_db_flushes,
1450 data_delta_between_two_snapsots,
1451 time_delta_between_two_snapsots,
1452 clear_database_after_flush,
1463 CCZ4 solver using Runge-Kutta Discontinuous Galerkin and global adaptive time stepping incl enclave tasking
1465 The constructor of this classs is straightforward and realises the standard
1466 steps of any numerical implementation of the CCZ4 scheme:
1468 1. Init the actual numerical scheme. This happens through the constructor
1471 2. Add the header files that we need, i.e. those files which contain the
1472 actual CCZ4 implementation.
1474 3. Add some constants that any CCZ4 C++ code requires.
1476 4. Set the actual implementation, i.e. link the generic PDE terms to the
1477 CCZ4-specific function calls.
1479 5. Add the CCZ4-specific postprocessing.
1481 6. Switch to higher-order interpolation and restriction.
1492 pde_terms_without_state,
1496 Construct solver with enclave tasking
1499 AbstractCCZ4Solver.__init__(self)
1500 exahype2.solvers.rkdg.rusanov.GlobalAdaptiveTimeStep.__init__(
1504 polynomials=polynomials,
1506 auxiliary_variables=0,
1507 min_cell_h=min_cell_h,
1508 max_cell_h=max_cell_h,
1509 time_step_relaxation=AbstractCCZ4Solver.Default_Time_Step_Size_Relaxation,
1510 pde_terms_without_state=pde_terms_without_state,
1516 boundary_conditions=exahype2.solvers.PDETerms.Empty_Implementation,
1518 flux=exahype2.solvers.PDETerms.None_Implementation,
1520 refinement_criterion=exahype2.solvers.PDETerms.Empty_Implementation,
1533 number_of_entries_between_two_db_flushes,
1534 data_delta_between_two_snapsots,
1535 time_delta_between_two_snapsots,
1536 clear_database_after_flush,
1541 Add tracer to project
1543 Consult exahype2.tracer.DumpTracerIntoDatabase for an explanation of
1544 some of the arguments. Most of them are simply piped through to this
1547 At this point, we have not yet created the Peano 4 project. Therefore, we
1548 have not yet befilled the time stepping action set.
1550 project: exahype2.Project
1558 number_of_entries_between_two_db_flushes,
1559 data_delta_between_two_snapsots,
1560 time_delta_between_two_snapsots,
1561 clear_database_after_flush,
Abstract base class for any CCZ4 solver.
dict _SO_formulation_unknowns
__init__(self)
Constructor.
float Default_Time_Step_Size_Relaxation
_add_standard_includes(self)
Add the headers for the compute kernels and initial condition implementations.
add_tracer(self, name, coordinates, project, number_of_entries_between_two_db_flushes, data_delta_between_two_snapsots, time_delta_between_two_snapsots, clear_database_after_flush, tracer_unknowns=None)
Add tracer to project.
add_all_solver_constants(self)
Add domain-specific constants.
enable_second_order(self)
Default_Time_Step_Size_Relaxation
dict _FO_formulation_unknowns
add_makefile_parameters(self, peano4_project, path_of_ccz4_application)
Add include path and minimal required cpp files to makefile.
CCZ4 solver using fourth-order finite differences and global adaptive time stepping incl enclave task...
add_tracer(self, name, coordinates, project, number_of_entries_between_two_db_flushes, data_delta_between_two_snapsots, time_delta_between_two_snapsots, clear_database_after_flush, tracer_unknowns)
Add tracer to project.
__init__(self, name, patch_size, rk_order, min_meshcell_h, max_meshcell_h, pde_terms_without_state, second_order=False)
Constructor.
CCZ4 solver using fourth-order finite differences and global adaptive time stepping without enclave t...
add_tracer(self, name, coordinates, project, number_of_entries_between_two_db_flushes, data_delta_between_two_snapsots, time_delta_between_two_snapsots, clear_database_after_flush, tracer_unknowns)
Add tracer to project.
__init__(self, name, patch_size, rk_order, min_meshcell_h, max_meshcell_h, second_order=False)
Constructor.
CCZ4 solver using finite volumes and global adaptive time stepping incl enclave tasking.
add_tracer(self, name, coordinates, project, number_of_entries_between_two_db_flushes, data_delta_between_two_snapsots, time_delta_between_two_snapsots, clear_database_after_flush, tracer_unknowns)
Add tracer to project.
__init__(self, name, patch_size, min_volume_h, max_volume_h, pde_terms_without_state)
Construct solver with enclave tasking and adaptive time stepping.
CCZ4 solver using finite volumes and global adaptive time stepping incl enclave tasking.
add_tracer(self, name, coordinates, project, number_of_entries_between_two_db_flushes, data_delta_between_two_snapsots, time_delta_between_two_snapsots, clear_database_after_flush, tracer_unknowns)
Add tracer to project.
__init__(self, name, patch_size, min_volume_h, max_volume_h, pde_terms_without_state)
Constructor.
CCZ4 solver using Runge-Kutta Discontinuous Galerkin and global adaptive time stepping incl enclave t...
__init__(self, name, rk_order, polynomials, min_cell_h, max_cell_h, pde_terms_without_state)
Construct solver with enclave tasking.
add_tracer(self, name, coordinates, project, number_of_entries_between_two_db_flushes, data_delta_between_two_snapsots, time_delta_between_two_snapsots, clear_database_after_flush, tracer_unknowns)
Add tracer to project.
CCZ4 solver using Runge-Kutta Discontinuous Galerkin and global adaptive time stepping incl enclave t...
__init__(self, name, rk_order, polynomials, min_cell_h, max_cell_h, pde_terms_without_state)
Construct solver with enclave tasking.
add_tracer(self, name, coordinates, project, number_of_entries_between_two_db_flushes, data_delta_between_two_snapsots, time_delta_between_two_snapsots, clear_database_after_flush, tracer_unknowns)
Add tracer to project.
set_implementation(self, boundary_conditions, refinement_criterion, initial_conditions, memory_location, use_split_loop, additional_action_set_includes, additional_user_includes)
If you pass in User_Defined, then the generator will create C++ stubs that you have to befill manuall...
postprocess_updated_patch(self)
postprocess_updated_patch(self, kernel)
Define a postprocessing routine over the data.
set_implementation(self, boundary_conditions, refinement_criterion, initial_conditions, memory_location, use_split_loop, additional_action_set_includes, additional_user_includes)
If you pass in User_Defined, then the generator will create C++ stubs that you have to befill manuall...
set_implementation(self, flux=None, ncp=None, eigenvalues=None, boundary_conditions=None, refinement_criterion=None, initial_conditions=None, source_term=None, memory_location=None, use_split_loop=False, additional_action_set_includes="", additional_user_includes="")
If you pass in User_Defined, then the generator will create C++ stubs that you have to befill manuall...
set_implementation(self, flux=None, ncp=None, eigenvalues=None, boundary_conditions=None, refinement_criterion=None, initial_conditions=None, source_term=None, memory_location=None, use_split_loop=False, additional_action_set_includes="", additional_user_includes="")
If you pass in User_Defined, then the generator will create C++ stubs that you have to befill manuall...
postprocess_updated_cell_after_final_linear_combination(self)
postprocess_updated_cell_after_final_linear_combination(self, kernel)
Define a postprocessing routine over the data.
set_implementation(self, flux=None, ncp=None, eigenvalues=None, boundary_conditions=None, refinement_criterion=None, initial_conditions=None, source_term=None, point_source=None, additional_action_set_includes="", additional_user_includes="")
If you pass in User_Defined, then the generator will create C++ stubs that you have to befill manuall...
RKDG solver with global adaptive time step.
set_implementation(self, flux=None, ncp=None, eigenvalues=None, boundary_conditions=None, refinement_criterion=None, initial_conditions=None, source_term=None, point_source=None, additional_action_set_includes="", additional_user_includes="")
If you pass in User_Defined, then the generator will create C++ stubs that you have to befill manuall...
RKDG solver with Rusanov Riemann solver employing global adaptive time stepping.
set_implementation(self, flux=None, ncp=None, eigenvalues=None, boundary_conditions=None, refinement_criterion=None, initial_conditions=None, source_term=None, point_source=None, additional_action_set_includes="", additional_user_includes="")
If you pass in User_Defined, then the generator will create C++ stubs that you have to befill manuall...
postprocess_updated_patch(self)
auxiliary_variables(self)
postprocess_updated_patch(self, kernel)
Define a postprocessing routine over the data.
auxiliary_variables(self, value)
preprocess_reconstructed_patch(self, kernel)
Please consult exahype2.solvers.fv.FV.preprocess_reconstructed_patch() for a documentation on this ro...
_action_set_compute_final_linear_combination
preprocess_reconstructed_patch(self)
set_implementation(self, flux, ncp, source_term, eigenvalues, boundary_conditions, refinement_criterion, initial_conditions, memory_location, additional_action_set_includes, additional_user_includes)
If you pass in User_Defined, then the generator will create C++ stubs that you have to befill manuall...
set_implementation(self, flux=None, ncp=None, source_term=None, eigenvalues=None, boundary_conditions=None, refinement_criterion=None, initial_conditions=None, memory_location=None, additional_action_set_includes="", additional_user_includes="", KOSigma=None)
If you pass in User_Defined, then the generator will create C++ stubs that you have to befill manuall...
set_implementation(self, flux=None, ncp=None, source_term=None, eigenvalues=None, boundary_conditions=None, refinement_criterion=None, initial_conditions=None, memory_location=None, additional_action_set_includes="", additional_user_includes="", KOSigma=None)
If you pass in User_Defined, then the generator will create C++ stubs that you have to befill manuall...
Particle tracing over the DG solver.
Dump the tracer data into a csv database.
Particle tracing over the Finite Volumes solver.
Basically superclass, though we add these numbers.
add_tracer_to_FV_solver(name, coordinates, project, solver, number_of_entries_between_two_db_flushes, data_delta_between_two_snapsots, time_delta_between_two_snapsots, clear_database_after_flush, tracer_unknowns)
Add tracer to project.
construct_FV_source_term()
add_tracer_to_DG_solver(name, coordinates, project, self, number_of_entries_between_two_db_flushes, data_delta_between_two_snapsots, time_delta_between_two_snapsots, clear_database_after_flush, tracer_unknowns)
construct_FV_eigenvalues()
add_tracer_to_FD4_solver(name, coordinates, project, solver, number_of_entries_between_two_db_flushes, data_delta_between_two_snapsots, time_delta_between_two_snapsots, clear_database_after_flush, tracer_unknowns)
I realise this as a separate routine, as we need it for all FD4 flavours.
construct_FD4_postprocessing_kernel()
construct_FD4_source_term()
construct_DG_eigenvalues()
construct_FV_postprocessing_kernel()
construct_FD4_eigenvalues()
construct_DG_postprocessing_kernel()
construct_DG_source_term()