rusanov-with-explicit-euler-kernel-microbenchmark.py

Go to the documentation of this file.

import os, sys
import argparse
import peano4.api
import exahype2
from exahype2.api.generators.KernelMicrobenchmarkMain import KernelMicrobenchmarkMain
 
parser = argparse.ArgumentParser(
    description="ExaHyPE 2 - RusanovWithExplicitEuler Microbenchmarking Script"
)
parser.add_argument(
    "-ps",
    "--patch-size",
    default=16,
    type=int,
    help="Number of finite volumes per axis (dimension) per patch.",
)
parser.add_argument(
    "-p",
    "--num-patches",
    nargs="+",
    default=[32, 64, 128, 256, 512, 1024, 2048],
    type=int,
    help="Number of patches to study",
)
parser.add_argument(
    "-samples",
    "--samples",
    default=10,
    type=int,
    help="Number of samples per measurement",
)
parser.add_argument(
    "-a",
    "--accuracy",
    default=0.0,
    type=float,
    help="Floating point accuracy to which the different kernel variants have to match (absolute). Pass in 0 to disable correctness check. Pass in values < 0 to use machine epsilon (default).",
)
parser.add_argument(
    "-cpu",
    "--cpu",
    action="store_true",
    help="Assess host kernels",
)
parser.add_argument(
    "-gpu",
    "--gpu",
    action="store_true",
    help="Assess device kernels",
)
args = parser.parse_args()
 
project = exahype2.api.Project(
    namespace=["exahype2", "fv", "rusanov", "microbenchmark"],
    project_name="kernel-microbenchmark",
    directory=".",
    executable="kernel-microbenchmark",
)
solver_name = "EulerFVSmooth"
 
addition_includes = f"""
#include "kernels/{solver_name}_rusanov.h"
"""
 
function_calls = ""
if args.cpu == True or args.gpu == False:
    function_calls += f"""
    assessKernel(wrapHostKernel<
        {solver_name}_timeStepWithRusanovStateless<
            Abstract{solver_name},
            Abstract{solver_name}::PatchSize,
            Abstract{solver_name}::HaloSize,
            Abstract{solver_name}::NumberOfUnknowns,
            Abstract{solver_name}::NumberOfAuxiliaryVariables
            >
        >,
        "host, AoS, serial, cpp, stateless",
        launchingThreads,
        ::tarch::accelerator::Device::HostDevice,
        localPatches
    );
    """
    function_calls += f"""
    assessKernel(wrapHostKernel<
        {solver_name}_timeStepWithRusanov<
            Abstract{solver_name}::PatchSize,
            Abstract{solver_name}::HaloSize,
            Abstract{solver_name}::NumberOfUnknowns,
            Abstract{solver_name}::NumberOfAuxiliaryVariables
            >
        >,
        "host, AoS, serial, cpp",
        launchingThreads,
        ::tarch::accelerator::Device::HostDevice,
        localPatches
    );
    """
 
 
if args.gpu:
    function_calls += f"""
    assessKernel(wrapDeviceKernel<
        {solver_name}_timeStepWithRusanovStateless<
            Abstract{solver_name},
            Abstract{solver_name}::PatchSize,
            Abstract{solver_name}::HaloSize,
            Abstract{solver_name}::NumberOfUnknowns,
            Abstract{solver_name}::NumberOfAuxiliaryVariables
            >
        >,
        "device, AoS, parallel, omp",
        launchingThreads,
        0,
        localPatches
    );
    """
 
project.mainfile_generator = KernelMicrobenchmarkMain(project, function_calls, addition_includes)
 
unknowns = 4
auxiliary_variables = 0
dimensions = 2
verbosity = peano4.api.LoggerVerbosity.OFF
 
thesolver = exahype2.api.Solver(
    solver_name,
    unknowns,
    auxiliary_variables,
)
 
thesolver.add_user_solver_includes("#include <cmath>\n" "#include <algorithm>\n")
 
for stateless in [True, False]:
    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=stateless,
        )
    )
    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 = gamma * p * irho;
 
        switch (normal) {
            case 0:
            lambdaMax[0] = std::fmax(std::fabs(Q[1] * irho + c), std::fabs(Q[1] * irho - c));
            break;
            case 1:
            lambdaMax[0] = std::fmax(std::fabs(Q[2] * irho + c), std::fabs(Q[2] * irho - c));
            break;
        }
        """,
            stateless=stateless,
        )
    )
thesolver.add_pde_term(exahype2.api.InitPointWisely(implementation = """
    constexpr double gamma       = 1.4;
    constexpr double p           = 1.0;
    constexpr double u           = 1.0;
    constexpr double domainLeft  = -0.1;
    constexpr double domainWidth = 0.2;
    const double     rho         = 1.0 + 0.2 * std::sin(2.0 * M_PI * (x[0] - domainLeft) / domainWidth);
    Q[0] = rho;
    Q[1] = rho * u;
    Q[2] = 0.0;
    Q[3] = p / (gamma - 1) + 0.5 * rho * u * u;
    """,
        stateless=False,
    )
)
 
thesolver.numerical_scheme = exahype2.api.fv.RusanovWithExplicitEuler(thesolver, args.patch_size)
solver_variants = []
if args.cpu == True or args.gpu == False:
    solver_variants.append((exahype2.api.Target.Serial, True))
    solver_variants.append((exahype2.api.Target.Serial, False))
if args.gpu == True:
    solver_variants.append((exahype2.api.Target.Accelerator, True))
 
thesolver.numerical_scheme._generator = exahype2.api.fv.RusanovWithExplicitEulerGenerator(thesolver, solver_variants)
thesolver.grid_construction = exahype2.api.fv.SetupSimulationOverRegularGrid(
    thesolver, volume_size=0.001
)
 
project.add_solver(thesolver)
 
build_mode = peano4.api.output.CompileMode.RELEASE
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=[True, True, True],
)
 
 
project.set_Peano4_installation("../../../../../", build_mode)
peano4_project = project.generate_Peano4_project(verbosity=verbosity)
 
peano4_project.constants.export_constexpr_with_type("Accuracy", str(args.accuracy), "double")
peano4_project.constants.export_constexpr_with_type("NumberOfSamples", str(args.samples), "int")
 
formatted_num_patches = "{{{}}}".format(", ".join(str(val) for val in args.num_patches))
peano4_project.constants.export_const_with_type("NumberOfPatchesToStudy", str(formatted_num_patches), "tarch::la::Vector<%s, int>" % len(args.num_patches))
peano4_project.constants.export_constexpr_with_type("NumberOfLaunchingThreads", str(1), "int")
 
peano4_project.generate()
peano4_project.build(make_clean_first=True)