noh.py

Go to the documentation of this file.

import os
import argparse
import numpy as np
 
import peano4.api
import exaparty
import exaparty.api.sphtools as sphtools
import dastgen2
 
import shutil
 
 
 
 
 
def parse_args():
    parser = argparse.ArgumentParser(description="SPH benchmarking script")
    parser.add_argument(
        "-j",
        "--parallel-builds",
        dest="j",
        type=int,
        default=-1,
        help="Parallel builds",
    )
    parser.add_argument(
        "-v",
        "--verbose",
        dest="verbose",
        action="store_true",
        default=False,
        help="Verbose",
    )
    parser.add_argument(
        "-c",
        "--clean",
        dest="clean_temporary_files",
        action="store_true",
        help="Remove temporary directories and output files before script starts",
    )
    parser.add_argument(
        "-np",
        "--particle-number",
        dest="particle_number",
        required=True,
        type=int,
        help="Particle number (1D)",
    )
    parser.add_argument(
        "-t",
        "--trees",
        dest="trees",
        default=0,
        type=int,
        help="Number of trees (subdomains) rank rank (0 uses all threads)",
    )
    parser.add_argument(
        "-et",
        "--end-time",
        dest="end_time",
        type=float,
        default=0.2,
        help="End of simulation",
    )
    parser.add_argument(
        "-o",
        "--output",
        dest="output_file",
        default="noh2D",
        help="Executable file name",
    )
    parser.add_argument(
        "-plot",
        "--plot-delta",
        dest="plot_delta",
        type=float,
        default=0.01,
        help="Plot delta (0 switches off)",
    )
    parser.add_argument(
        "-dt",
        "--timestep-size",
        dest="timestep_size",
        type=float,
        default=1e-4,
        help="Global timestep size (initial value)",
    )
    parser.add_argument(
        "-cfl",
        "--cfl-factor",
        dest="cfl_factor",
        type=float,
        default=0.1,
        help="Value of the CFL constant",
    )
    parser.add_argument(
        "-ppc",
        "--particles-per-cell",
        dest="particles_per_cell",
        type=int,
        default=100,
        help="Particles per cell. If we have x particles, we will have 2^d*x particle comparisons on average per cell, as we compare to neighbouring half cells as well",
    )
    parser.add_argument(
        "-pd",
        "--peano-dir",
        dest="peanodir",
        default="../../../..",
        help="Peano4 directory",
    )
    parser.add_argument(
        "-m",
        "--mode",
        dest="mode",
        choices=["release", "stats", "trace", "asserts", "debug"],
        default="release",
        help="Pick build type",
    )
    parser.add_argument(
        "-rea",
        "--realisation",
        dest="realisation",
        default="domain-decomposition.scattered.multiscale-sorting.no-optimisation",
        help="Realisation to use. See Peano/src/exaparty/api/ProjectRealisation.py for options",
    )
    parser.add_argument(
        "-ms",
        "--mantissa-size",
        dest="mantissa_size",
        type=int,
        default=0,
        help="Mantissa size for clang truncation (if supported by compiler). Default of 0 disables any compression",
    )
    parser.add_argument(
        "-k",
        "--kernel",
        dest="kernel",
        type=str,
        default="quartic_spline",
        choices=sphtools.sph_kernel_list,
        help="SPH kernel function to use.",
    )
 
    return parser.parse_args()
 
 
if __name__ == "__main__":
    args = parse_args()
 
 
 
project_namespace = ["benchmarks", "exaparty", "noh2D"]
 
project = exaparty.api.Project(
    namespace=project_namespace, project_name="Noh2D", executable=args.output_file
)
 
# Initialise the SPH particle --------------------------------------------------
# If multiple particle species are present, we can specify the parameters
# independently for each case.
 
name = "HydroPart"
 
dimensions = 2
hydro_dimensions = 2  # dimensions for the hydro
 
particles_per_cell = args.particles_per_cell
total_number_of_particles = args.particle_number**hydro_dimensions
total_number_of_cells = total_number_of_particles / particles_per_cell
cells_per_dimension = total_number_of_cells ** (1.0 / hydro_dimensions)
max_cell_size = 2.0 / cells_per_dimension * 3.0
 
print(
    "{} particles in total fit into {} cells, i.e. {} cells per dimension. This gives h={}".format(
        total_number_of_particles,
        total_number_of_cells,
        cells_per_dimension,
        max_cell_size,
    )
)
 
particle = exaparty.api.particle.SPHLeapfrogFixedSearchRadius(
    name=name,
    dimensions_hydro=hydro_dimensions,
    cfl_factor=args.cfl_factor,
    initial_time_step_size=args.timestep_size,
    constant_time_step_size=True,
    exaparty_project_namespace=project_namespace,
    particles_per_cell=particles_per_cell,
    min_h=max_cell_size,
    max_h=max_cell_size,
)
 
 
 
offset = [0, 0, 0]
domain_size = [2, 2, 2]
periodic_boundary_conditions = [False, False, False]
 
# HYDRO_DIMENSION can be smaller than the Peano4 dimensions
HYDRO_DIMENSIONS = hydro_dimensions
 
# Particle number along 1 dimension
HYDRO_PART_NUMBER = args.particle_number
 
# Time integration parameters
GLOBAL_TIME_STEP_SIZE = args.timestep_size
CFL_CONSTANT = args.cfl_factor
 
# EoS parameter ---------------------------------------------------------------
gamma_hydro_list = [str(5 / 3), str(7 / 5), str(4 / 3), str(2 / 1)]
gamma_hydro_list_symbols = [
    "HYDRO_GAMMA_5_3",
    "HYDRO_GAMMA_7_5",
    "HYDRO_GAMMA_4_3",
    "HYDRO_GAMMA_2_1",
]
GAMMA = gamma_hydro_list[0]
GAMMA_HYDRO_SYMBOL = gamma_hydro_list_symbols[0]
 
# SPH kernel -------------------------------------------------------
 
kernel = args.kernel
KERNEL_SYMBOL = sphtools.sph_kernel_macro_name[kernel]
KERNEL_SUPPORT_RAD = sphtools.sph_kernel_H_over_h[HYDRO_DIMENSIONS - 1][kernel]
 
hydro_H_max = max_cell_size / 2
 
dx = domain_size[0] / args.particle_number
if hydro_H_max < 3.0 * dx:
    raise ValueError("Too small hydro_H_max", 3.0 * dx, hydro_H_max)
hydro_h_max = exaparty.api.sphtools.sph_kernel_get_smoothing_length_from_interaction_radius(
    hydro_H_max, kernel, hydro_dimensions
)
particle.h_hydro_max = hydro_h_max
 
# particle.h_hydro_min = 1e-6
# particle.h_max_iterations = 10
# particle.h_tolerance = 1e-4
# particle.eta_factor = 1.2348
# particle.beta_av = 3.0
# particle.alpha_av = 0.8
 
# IMPORTANT: Reduce the mantissa size for particle attributes -------------------
# Native size is 53 (double precision)
if args.mantissa_size > 0:
    particle.mass.valid_mantissa_bits = args.mantissa_size
    particle.velocity.valid_mantissa_bits = args.mantissa_size
    particle.acceleration.valid_mantissa_bits = args.mantissa_size
    particle.density.valid_mantissa_bits = args.mantissa_size
    particle.pressure.valid_mantissa_bits = args.mantissa_size
    particle.smoothL.valid_mantissa_bits = args.mantissa_size
    particle.u.valid_mantissa_bits = args.mantissa_size
    particle.uDot.valid_mantissa_bits = args.mantissa_size
    particle.f.valid_mantissa_bits = args.mantissa_size
    particle.wcount_dh.valid_mantissa_bits = args.mantissa_size
    particle.rho_dh.valid_mantissa_bits = args.mantissa_size
    particle.hleft.valid_mantissa_bits = args.mantissa_size
    particle.hright.valid_mantissa_bits = args.mantissa_size
    particle.smoothL_prev.valid_mantissa_bits = args.mantissa_size
    particle.wcount.valid_mantissa_bits = args.mantissa_size
    particle.hDot.valid_mantissa_bits = args.mantissa_size
    particle.residual.valid_mantissa_bits = args.mantissa_size
    particle.balsara.valid_mantissa_bits = args.mantissa_size
    particle.rot_v.valid_mantissa_bits = args.mantissa_size
    particle.div_v.valid_mantissa_bits = args.mantissa_size
    particle.v_sig_AV.valid_mantissa_bits = args.mantissa_size
    particle.soundSpeed.valid_mantissa_bits = args.mantissa_size
    particle.v_full.valid_mantissa_bits = args.mantissa_size
    particle.u_full.valid_mantissa_bits = args.mantissa_size
 
particle.data.expose_all_attributes_in_header_file()
 
# deploy particle properties
# This needs to be done if you're changing any of the default values.
# Otherwise, the changes won't propagate into dastgen.
particle.set_parameters()
 
# Add particle to the ExaParty project --------------------------------------------
# Notice that you any change to the particle will not take effect after adding
# the particle_set to the project
particle_set = project.add_particle_species(particle)
 
if args.mode == "release":
    build_mode = peano4.api.output.CompileMode.RELEASE
if args.mode == "stats":
    build_mode = peano4.api.output.CompileMode.STATS
if args.mode == "asserts":
    build_mode = peano4.api.output.CompileMode.ASSERTS
if args.mode == "trace":
    build_mode = peano4.api.output.CompileMode.TRACE
if args.mode == "debug":
    build_mode = peano4.api.output.CompileMode.DEBUG
 
# Set Peano4 simulation parameters ---------------------------------------------
project.set_global_simulation_parameters(
    dimensions=dimensions,
    offset=offset,
    domain_size=domain_size,
    min_end_time=args.end_time,
    max_end_time=0.0,
    first_plot_time_stamp=0.0,
    time_in_between_plots=args.plot_delta,
    periodic_BC=periodic_boundary_conditions,
    plotter_precision=8,
)
 
# Set realisation
realisation = args.realisation.split(".")
parallelisation = realisation[0]
storage = realisation[1]
sorting = realisation[2]
kernel_optimisation = realisation[3]
 
project.set_realisation(
    parallelisation=parallelisation,
    storage=storage,
    sorting=sorting,
    kernel_optimisation=kernel_optimisation,
)
 
project.set_Peano4_installation(args.peanodir, build_mode)
 
project.set_load_balancing(
    "toolbox::loadbalancing::strategies::SpreadOutOnceGridStagnates",
    f"new toolbox::loadbalancing::DefaultConfiguration({args.trees})",
)
 
 
 
 
#  Noh problem (2D) snippet  --------------------------------------------------
#  The particle attributes are set via an initalisation snippet that sets
#  the initial conditions specific for this problem, i.e. the velocity field
#  that converges towards the center of the domain and that generates the
#  collapse.
 
ICs = exaparty.api.scenario.NohProblemIC(dimensions=HYDRO_DIMENSIONS)
initialisation_snippet = ICs.initialisation_snippet
 
project.algorithm_step_initial_conditions.add_action_set(
    exaparty.api.input.InsertParticlesAlongCartesianGrid(
        particle_set=particle_set,
        initialisation_call=initialisation_snippet,
        distance_between_particles_along_one_axis=(
            domain_size[0] / (args.particle_number)
        ),
        additional_includes="""
#include "exaparty/kernels/legacy/equation_of_state.h"
#include "exaparty/kernels/legacy/kernel_hydro.h"
""",
    )
)
 
 
particle_plotter = peano4.api.toolbox.particles.PlotParticlesInVTKFormat(
    "particles",
    particle_set,
    "globaldata::{}::getSpecies().getMinTimeStamp()".format(particle.name),
)
particle_plotter.add_attribute_to_plot(particle.velocity, dimensions)
particle_plotter.add_attribute_to_plot(particle.acceleration, dimensions)
particle_plotter.add_attribute_to_plot(particle.mass, 1)
particle_plotter.add_attribute_to_plot(particle.density, 1)
particle_plotter.add_attribute_to_plot(particle.pressure, 1)
particle_plotter.add_attribute_to_plot(particle.u, 1)
particle_plotter.add_attribute_to_plot(particle.smoothL, 1)
if args.mode == "asserts" or args.mode == "debug":
    particle_plotter.add_attribute_to_plot(particle.partid, 1)
    particle_plotter.add_attribute_to_plot(particle.density_neighbourcount, 1)
    particle_plotter.add_attribute_to_plot(particle.force_neighbourcount, 1)
project.algorithm_step_plot.add_action_set(particle_plotter)
 
 
 
 
peano4_project = project.generate_Peano4_project(verbose=args.verbose)
 
# Some sanity checks ----------------------------------------------------------
if GAMMA not in gamma_hydro_list:
    print(f"Please check the value of GAMMA. You have chosen: {GAMMA}")
 
print("-----------------------------------------------------------------------")
print("SPH parameters summary:")
print(f"Effective dimensions of the Simulation: {HYDRO_DIMENSIONS}")
print(f"Global time-step size: {GLOBAL_TIME_STEP_SIZE}")
print(f"Number of particles (1D): {HYDRO_PART_NUMBER}")
print(f"Adiabatic index: {GAMMA}")
print(f"Smoothing length free factor: {particle.eta_factor}")
print("-----------------------------------------------------------------------")
 
 
# Export Constants.h file ------------------------------------------------------
peano4_project.constants.export_const_with_type(
    "HYDRO_DIMENSIONS", str(HYDRO_DIMENSIONS), "double"
)
peano4_project.constants.export_const_with_type(
    "HYDRO_PART_NUMBER", str(HYDRO_PART_NUMBER), "double"
)
peano4_project.constants.export_const_with_type(
    "GLOBAL_TIME_STEP_SIZE", str(GLOBAL_TIME_STEP_SIZE), "double"
)
peano4_project.constants.export_const_with_type("GAMMA", GAMMA, "double")
 
# Export symbols for compilation (e.g. for SPH kernel) -----------------------
peano4_project.output.makefile.add_CXX_flag(
    "-D" + "HYDRO_DIMENSION=" + str(HYDRO_DIMENSIONS)
)
peano4_project.output.makefile.add_CXX_flag("-D" + GAMMA_HYDRO_SYMBOL)
peano4_project.output.makefile.add_CXX_flag("-D" + KERNEL_SYMBOL)
 
 
 
# If we want to start a project from scratch, we may want to remove some of the generated code
if args.clean_temporary_files:
    dirs = ["vertexdata/", "repositories/", "globaldata/", "observers/"]
    for dir in dirs:
        shutil.rmtree(dir)
 
# It is safe (and recommended) to overwrite any existing main file, since this
# is generated automatically by the ExaParty graph compiler.
peano4_project.generate(
    overwrite=peano4.api.output.Overwrite.Always, throw_away_data_after_generation=True
)
peano4_project.build(make_clean_first=True, number_of_parallel_builds=args.j)
 
if args.clean_temporary_files:
    # remove output files from previous runs
    output_files = [
        f
        for f in os.listdir(".")
        if f.endswith(".peano-patch-file") or f.endswith(".vtu") or f.endswith(".pvd")
    ]
    for f in output_files:
        os.remove(f)