sedov.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
 
 
 
 
parser = argparse.ArgumentParser(description="SPH test 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_directories",
    action="store_true",
    help="Remove temporary directories before script starts",
)
parser.add_argument(
    "-np",
    "--particle-number",
    dest="particle_number",
    required=True,
    type=int,
    help="Particle number (1D)",
)
parser.add_argument(
    "-et",
    "--end-time",
    dest="end_time",
    type=float,
    default=0.2,
    help="End of simulation",
)
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(
    "-csmax",
    "--cell-size-max",
    dest="cell_size_max",
    type=float,
    default=0.3,
    help="Cell size (default 0.3)",
)
parser.add_argument(
    "-csmin",
    "--cell-size-min",
    dest="cell_size_min",
    type=float,
    default=0.3,
    help="Cell size (default 0.3)",
)
parser.add_argument(
    "-pd",
    "--peano-dir",
    dest="peanodir",
    default="../../../..",
    help="Peano4 directory",
)
parser.add_argument(
    "-asserts",
    "--asserts",
    dest="asserts",
    action="store_true",
    default=False,
    help="Switch on assertions",
)
parser.add_argument(
    "-debug",
    "--debug",
    dest="debug",
    action="store_true",
    default=False,
    help="Switch on full debug checks",
)
parser.add_argument(
    "-k",
    "--kernel",
    dest="kernel",
    type=str,
    default="quartic_spline",
    choices=sphtools.sph_kernel_list,
    help="SPH kernel function to use.",
)
parser.add_argument(
    "-t",
    "--trees",
    dest="trees",
    default=0,
    type=int,
    help="Number of trees (subdomains) (0 uses all threads)",
)
 
args = parser.parse_args()
 
 
 
 
project_namespace = ["tests", "exaparty", "sedov1D"]
 
project = exaparty.api.Project(
    namespace=project_namespace, project_name="Sedov1D", executable="sedov1D"
)
 
# 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 = 1  # dimensions for the hydro
 
particle = exaparty.api.particle.SPHLeapfrogFixedSearchRadius(
    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,
    min_h=args.cell_size_min,
    max_h=args.cell_size_max,
)
 
 
 
 
offset = [0, 0, 0]
domain_size = [1, 1, 1]
periodic_boundary_conditions = [True, 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]
 
 
# Additional Particle Parameters -----------------------------------------------
depth = 0
cellsize = domain_size[0]
while cellsize > args.cell_size_min:
    cellsize /= 3.0
    depth += 1
    if depth == 100:
        raise ValueError("Suspicious depth")
 
dx = domain_size[0] / args.particle_number
hydro_H_max = 0.499 * cellsize
if hydro_H_max < 5.0 * dx:
    raise ValueError("Too small hydro_H_max", 5.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
 
# Optional: Reduce the mantissa size for particle attributes -------------------
# Native size is 53 (double precision)
# particle.mantissa_size = 23
 
 
# deploy particle properties
# This needs to be done if you're changing any of the default values.
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.debug:
    build_mode = peano4.api.output.CompileMode.DEBUG
elif args.asserts:
    build_mode = peano4.api.output.CompileMode.ASSERTS
else:
    build_mode = peano4.api.output.CompileMode.RELEASE
 
# 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,
)
 
# Choose realisation
parallelisation = "domain-decomposition"
#  parallelisation = "task-graph"
#  parallelisation = "multisweep-task-graph"
storage = "scattered"
#  storage = "continuous-per-vertex"
#  storage = "global-continuous"
#  sorting = "multiscale-sort"
sorting = "bucket-sort"
kernel_optimisation = "no-optimisation"
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})",
)
 
 
 
 
#  Sedov-Taylor Blast 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.SedovBlastIC()
initialisation_snippet = ICs.initialisation_snippet
 
coords = np.ones((args.particle_number, 2)) * 0.5
for i in range(args.particle_number):
    coords[i, 0] = (i + 0.5) / (args.particle_number + 1)
 
# Finally, pass this configuration to the initialisation script
# Note that Peano's internal cartesian grid particle generation can't do 1D grids at the moment.
project.algorithm_step_initial_conditions.add_action_set(
    exaparty.api.input.InsertParticlesByCoordinates(
        particle_set=particle_set,
        coordinates=coords,
        initialisation_call=initialisation_snippet,
        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),
)
# for some reason, using hydro_dimensions here writes wrong output.
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.asserts or args.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_directories:
    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)
 
# 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)