pml-dynamic-rupture.py

Go to the documentation of this file.

import os, sys, peano4, exahype2
 
from exahype2.solvers.aderdg.ADERDG import Polynomials
 
sys.path.insert(0, os.path.abspath("../pml"))
sys.path.insert(0, os.path.abspath("../dynamic-rupture"))
 
import PMLElastic
import DynamicRuptureElastic
import PMLDynamicRuptureElastic
 
sys.path.insert(0, os.path.abspath("../ExaSeis_core"))
import Scenario
 
available_scenarios = {
    "TPV5":     Scenario.TPV5(),
    "TPV6":     Scenario.TPV6(),
    "TPV16":    Scenario.TPV16(),
    "TPV26":    Scenario.TPV26(),
    "TPV28":    Scenario.TPV28(),
    "TPV29":    Scenario.TPV29(),
    "TPV34":    Scenario.TPV34(),
    "Husavik":  Scenario.Husavik()
}
 
scenario = "TPV28"
scenario_string = scenario.lower()
my_scenario = available_scenarios[scenario]
 
unknowns = {"v": 3, "sigma": 6, "u": 3, "pml": 27}
auxiliary_variables = {
    "rho": 1, "cp": 1, "cs": 1,
    "dmp_pml": 3, "jacobian": 1,
    "metric_derivative": 9,
    "curve_grid": 3,
}
 
offset   = my_scenario.domain_offset
size     = my_scenario.domain_size
end_time = my_scenario.end_time
plot_dt  = 0.0
order = 5
min_level = 3
max_depth = 0
max_h = 1.1 * min(size) / (3.0**min_level)
min_h = max_h / (3.0**max_depth)
 
precision = "fp64"
 
my_scenario.generate_required_files(order)
 
theSolver = exahype2.solvers.aderdg.GlobalAdaptiveTimeStep(
    name="ElasticSolver",
    order=order,
    unknowns=unknowns,
    auxiliary_variables=auxiliary_variables,
    min_cell_h=min_h,
    max_cell_h=max_h,
    time_step_relaxation=0.9
)
 
theSolver.add_kernel_optimisations(
    is_linear=True,
    polynomials=Polynomials.Gauss_Lobatto,
    initialise_patches=True,
    precision=precision
)
 
theSolver.set_implementation(
    initial_conditions    = PMLElastic.initial(my_scenario.initial_conditions()),
    boundary_conditions   = PMLElastic.boundary(),
    max_eigenvalue        = PMLElastic.eigenvalue(),
    flux                  = PMLElastic.flux(),
    ncp                   = PMLElastic.ncp(),
    material_parameters   = PMLElastic.multiplyMaterialParameterMatrix(),
    riemann_solver        = DynamicRuptureElastic.riemann_solver(),
    source_term           = PMLElastic.algebraic_source()+PMLDynamicRuptureElastic.algebraic_source()
)
 
theSolver._abstract_solver_user_definitions  += DynamicRuptureElastic.abstractUserDefinitions()
theSolver._abstract_solver_user_declarations += PMLDynamicRuptureElastic.abstractDeclarations()
theSolver._start_grid_initialisation_implementation += DynamicRuptureElastic.init_grid_step_implementation(scenario_string)
theSolver._destructor_implementation += PMLElastic.abstractDestructor()
theSolver.add_user_solver_includes(PMLDynamicRuptureElastic.userIncludes())
 
project = exahype2.Project(["applications", "exahype2", "exaseis"], "pml", executable=scenario)
 
project.add_solver(theSolver)
 
filename = scenario+"_l_"+str(min_level)+"_o_"+str(order)+"_pr_"+str(precision)
 
project.add_action_set_to_timestepping(
    exahype2.tracer.CurviTracer(
    coordinates= my_scenario.tracer_coordinates,
    fault_coordinates=my_scenario.fault_coordinates,
    solver=theSolver,
    filename=filename,
    data_delta_between_two_snapshots=1e16, time_delta_between_two_snapshots=0.01,
    output_precision=10, clear_database_after_flush=False,
    unknowns_to_store=12
))
 
if plot_dt > 0.0:
  project.set_output_path("solutions"+filename)
 
project.set_global_simulation_parameters(
    dimensions=3,
    offset=offset,
    size=size,
    min_end_time=end_time,
    max_end_time=end_time,
    first_plot_time_stamp=0.0,
    time_in_between_plots=plot_dt,
    periodic_BC=[False, False, False],
)
 
project.set_load_balancer("new ::exahype2::LoadBalancingConfiguration()")
project.set_build_mode(peano4.output.CompileMode.Release)
peano4_project = project.generate_Peano4_project(verbose=False)
peano4_project.build(make_clean_first=True)