smooth-convergence-test-with-injected-functions.py

Go to the documentation of this file.

import os, sys
import peano4.api
import exahype2
 
 
project = exahype2.api.Project(
    namespace=["benchmarks", "exahype2", "euler", "smooth_convergence_test"],
    project_name="smooth-convergence-test",
    directory=".",
    executable="smooth-convergence-test",
)
 
solver_name = "EulerFVSmooth"
unknowns = 4
auxiliary_variables = 0
mesh_size = 0.2 / 9.0
patch_size = 9
dimensions = 2
verbosity = peano4.api.LoggerVerbosity.TEXT
# verbosity   = peano4.api.LoggerVerbosity.VISUALISATION
 
thesolver = exahype2.api.Solver(
    solver_name,
    unknowns,
    auxiliary_variables,
)
 
thesolver.add_user_solver_includes("#include <cmath>\n" "#include <algorithm>\n")
 
thesolver.add_pde_term(
    exahype2.api.Flux(
        implementation="""
      const double irho  = 1.0 / Q[0];
      constexpr double gamma = 1.4;
      const double p = (gamma - 1.0) * (Q[3] - 0.5 * irho * (Q[1] * Q[1] + Q[2] * Q[2]));
 
      switch (normal) {
        case 0:
          F[0] = Q[1];
          F[1] = irho * Q[1] * Q[1] + p;
          F[2] = irho * Q[1] * Q[2];
          F[3] = irho * Q[1] * (Q[3] + p);
          break;
        case 1:
          F[0] = Q[2];
          F[1] = irho * Q[2] * Q[1];
          F[2] = irho * Q[2] * Q[2] + p;
          F[3] = irho * Q[2] * (Q[3] + p);
          break;
      }
    """,
        stateless=False,
    )
)
thesolver.add_pde_term(
    exahype2.api.MaxEigenvalue(
        implementation="""
      // Eigenvalues of the Euler flux Jacobian in direction "normal" are
      // u_n +/- c, where u_n is the *velocity* component (momentum/density),
      // not the raw momentum component - several earlier ExaHyPE
      // applications/benchmarks used the raw momentum here by mistake, which
      // silently gives a wrong (and usually too small) CFL time step bound.
      const double irho  = 1.0 / Q[0];
      constexpr double gamma = 1.4;
      const double p = (gamma - 1.0) * (Q[3] - 0.5 * irho * (Q[1] * Q[1] + Q[2] * Q[2]));
      const double c = std::sqrt(gamma * p * irho);
 
      switch (normal) {
        case 0:
          lambdaMax[0] = std::max(std::abs(Q[1] * irho + c), std::abs(Q[1] * irho - c));
          break;
        case 1:
          lambdaMax[0] = std::max(std::abs(Q[2] * irho + c), std::abs(Q[2] * irho - c));
          break;
      }
    """,
        stateless=True,
    )
)
thesolver.add_pde_term(exahype2.api.InitPointWisely(stateless=False,))
 
thesolver.numerical_scheme = exahype2.api.fv.RusanovWithExplicitEuler(thesolver, patch_size)
thesolver.grid_construction = exahype2.api.fv.SetupSimulationOverRegularGrid(
    thesolver, volume_size=0.001
)
thesolver.plotter = exahype2.api.fv.PlotAsVTKFile(
    thesolver, first_plot_time_stamp=0.0, time_in_between_plots=0.01
)
 
#  patch_size,
#  mesh_size,
#  mesh_size,
#  normalised_time_step_size=0.0001,
#  flux = exahype2.solvers.PDETerms.User_Defined_Implementation,
#  eigenvalues = exahype2.solvers.PDETerms.User_Defined_Implementation,
 
 
# oldsolver = exahype2.solvers.fv.musclhancock.GlobalFixedTimeStep(
#  solver_name,
#  patch_size,
#  unknowns, auxiliary_variables,
#  mesh_size,
#  mesh_size,
#  normalised_time_step_size=0.0001,
#  flux = exahype2.solvers.PDETerms.User_Defined_Implementation,
#  eigenvalues = exahype2.solvers.PDETerms.User_Defined_Implementation,
# )
 
project.add_solver(thesolver)
 
# build_mode = peano4.api.output.CompileMode.RELEASE
build_mode = peano4.api.output.CompileMode.DEBUG
project.set_global_simulation_parameters(
    dimensions=dimensions,
    offset=[-0.1, -0.1],
    size=[0.2, 0.2],
    min_end_time=0.01,
    max_end_time=0.01,
    periodic_BC=[False, False, False],
)
 
# project.set_load_balancing("toolbox::loadbalancing::strategies::RecursiveBipartition", "new ::exahype2::LoadBalancingConfiguration()")
project.set_Peano4_installation("../../../../", build_mode)
peano4_project = project.generate_Peano4_project(verbosity=verbosity)
peano4_project.generate()
peano4_project.build(make_clean_first=True)