static-limiting-swe.py

Go to the documentation of this file.

# This file is part of the ExaHyPE2 project. For conditions of distribution and
# use, please see the copyright notice at www.peano-framework.org
import peano4
import exahype2
 
project = exahype2.Project(
    ["tests", "exahype2", "aderdg"], ".", executable="ADERDG-Limiting"
)
 
parser = exahype2.ArgumentParser("ExaHyPE 2 - ADER-DG Limiting Testing Script")
parser.set_defaults(
    min_depth=3,
    degrees_of_freedom=6,
)
args = parser.parse_args()
 
dimensions = 2
order = args.degrees_of_freedom - 1
end_time = 1.0
plot_interval = 0.0  # 0.1
 
size = [2.0, 2.0, 2.0][0:dimensions]
offset = [-1.0, -1.0, -1.0][0:dimensions]
max_h = 1.1 * min(size) / (3.0**args.min_depth)
min_h = max_h / (3.0**args.amr_levels)
 
initial_conditions = """
  Q[0] = 4.0;
  Q[1] = 0.0;
  Q[2] = 0.0;
  Q[3] = 2.0 - tarch::la::norm2(x);
"""
 
boundary_conditions = """
  Qoutside[0] = Qinside[0];
  Qoutside[1] = Qinside[1];
  Qoutside[2] = Qinside[2];
  Qoutside[3] = Qinside[3];
"""
 
flux = """
  double ih = 1.0 / Q[0];
 
  F[0] = Q[1 + normal];
  F[1] = Q[1 + normal] * Q[1] * ih;
  F[2] = Q[1 + normal] * Q[2] * ih;
  F[3] = 0.0;
"""
 
ncp = """
  constexpr double grav = 9.81;
 
  BTimesDeltaQ[0] = 0.0;
  switch (normal) {
  case 0:
    BTimesDeltaQ[1] = grav * Q[0] * (deltaQ[0] + deltaQ[3]);
    BTimesDeltaQ[2] = 0.0;
    break;
  case 1:
    BTimesDeltaQ[1] = 0.0;
    BTimesDeltaQ[2] = grav * Q[0] * (deltaQ[0] + deltaQ[3]);
    break;
  }
  BTimesDeltaQ[3] = 0.0;
"""
 
max_eval = """
  constexpr double grav = 9.81;
 
  double u = Q[1 + normal] / Q[0];
  double c = std::sqrt(grav * Q[0]);
 
  return std::max(std::abs(u + c), std::abs(u - c));
"""
 
limiting_criterion = """
bool isBoundaryCell = (
  x[0] + 1.00 <  0.5 * h[0] or
  x[1] + 1.00 <  0.5 * h[1] or
  x[0] - 1.00 > -0.5 * h[0] or
  x[1] - 1.00 > -0.5 * h[1]
);
return !isBoundaryCell;
"""
 
aderdg_solver = exahype2.solvers.aderdg.GlobalAdaptiveTimeStep(
    name="ADERDGSolver",
    order=order,
    unknowns=4,
    auxiliary_variables=0,
    min_cell_h=min_h,
    max_cell_h=max_h,
    time_step_relaxation=0.9,
)
aderdg_solver.set_implementation(
    flux=flux,
    ncp=ncp,
    max_eigenvalue=max_eval,
    initial_conditions=initial_conditions,
    boundary_conditions=boundary_conditions,
)
aderdg_solver.add_kernel_optimisations(
    polynomials=exahype2.solvers.aderdg.Polynomials.Gauss_Lobatto
)
project.add_solver(aderdg_solver)
 
fv_solver = exahype2.solvers.fv.godunov.GlobalAdaptiveTimeStep(
    name="FVSolver",
    patch_size=2 * order + 1,
    unknowns=4,
    auxiliary_variables=0,
    min_volume_h=min_h,
    max_volume_h=max_h,
    time_step_relaxation=0.9,
 
)
fv_solver.set_implementation(
    flux=flux,
    ncp=ncp,
    max_eigenvalue=max_eval,
    initial_conditions=initial_conditions,
    boundary_conditions=boundary_conditions,
)
project.add_solver(fv_solver)
 
limiter = exahype2.solvers.limiting.StaticLimiting(
    name="LimitingSolver",
    regular_solver=aderdg_solver,
    limiting_solver=fv_solver,
    physical_admissibility_criterion=limiting_criterion,
)
project.add_solver(limiter)
 
project.set_output_path("solutions")
 
project.set_global_simulation_parameters(
    dimensions=dimensions,
    offset=offset[0:dimensions],
    size=size[0:dimensions],
    min_end_time=end_time,
    max_end_time=end_time,
    first_plot_time_stamp=0.0,
    time_in_between_plots=plot_interval,
    periodic_BC=[False, False, False][0:dimensions],
)
 
project.set_load_balancer("new ::exahype2::LoadBalancingConfiguration")
project.set_Peano4_installation("../../../", mode=peano4.output.string_to_mode(args.build_mode))
project = project.generate_Peano4_project(verbose=False)
project.set_fenv_handler(False)
project.build()