test_periodic_BC-1D.py

Go to the documentation of this file.

import os
import argparse
import numpy as np
 
import peano4.api
import exaparty
import dastgen2
 
import shutil
 
 
 
parser = argparse.ArgumentParser(description="Test periodic boundary conditions")
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(
    "-et",
    "--end-time",
    dest="end_time",
    type=float,
    default=1.0,
    help="End of simulation",
)
parser.add_argument(
    "-dt",
    "--timestep-size",
    dest="timestep_size",
    type=float,
    default=1e-3,
    help="Global timestep size",
)
parser.add_argument(
    "-plot",
    "--plot-delta",
    dest="plot_delta",
    type=float,
    default=0.01,
    help="Plot delta (0 switches off)",
)
parser.add_argument(
    "-cs",
    "--cell-size",
    dest="cell_size",
    type=float,
    default=0.1,
    help="Cell size (default 0.1)",
)
parser.add_argument(
    "-pd",
    "--peano-dir",
    dest="peanodir",
    default="../../../../",
    help="Peano4 directory",
)
parser.add_argument(
    "-a",
    "--asserts",
    dest="asserts",
    action="store_true",
    default=False,
    help="Switch on assertions",
)
parser.add_argument(
    "-d",
    "--debug",
    dest="debug",
    action="store_true",
    default=False,
    help="Switch on debug assertions",
)
args = parser.parse_args()
 
 
 
 
# peano doesn't like 1D grids
dimensions = 2
hydro_dimensions = 1
 
 
project = exaparty.api.Project(
    namespace=["tests", "exaparty", "periodicBoundaryConditions1D"],
    project_name="periodicBoundaryConditions1D",
    executable="periodicBCTest1D",
)
 
 
particle = exaparty.api.particle.LeapfrogFixedSearchRadius(
    name="GenericParticle",
    cfl_factor=0.999,
    initial_time_step_size=args.timestep_size,
    particles_per_cell=0,
    min_h=args.cell_size,
    max_h=0.3,
    add_partid=(not args.asserts),
)
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
 
offset = [0, 0, 0]
domain_size = [1, 1, 1]
periodic_boundary_conditions = [True, False, False]  # Periodic BC
 
 
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,  # plotter precision
)
 
project.set_load_balancing(
    "toolbox::loadbalancing::strategies::SpreadOutHierarchically",
    "new toolbox::loadbalancing::DefaultConfiguration()",
)
 
project.set_Peano4_installation(args.peanodir, build_mode)
 
 
 
 
particle_force = ""
 
#  particle.cell_kernel = particle_force
 
 
 
search_radius = args.cell_size / 100.0
initialisation_snippet = (
    """
  /* Internal search radius used by Peano.
   * Not really used here. */
  particle->setSearchRadius("""
    + str(search_radius)
    + """);
  particle->setA( {0.} );
 
  if( particle->getX(0) > 0.5 ){
    particle->setV( {-1.} );
  } else {
    particle->setV( {1.} );
  }
 
  // hack to get a somewhat flexible, but always unique
  // and reproducible particle ID: Assume we have a 1600^3
  // uniform grid particle distribution. (1600^3 is close)
  // to the numeric limits of integers.
 
  double boxsize_x = """
    + str(domain_size[0])
    + """;
  double offset_x = """
    + str(offset[0])
    + """;
 
  double dx = boxsize_x / 1600.;
 
  int id = (int) ((particle->getX(0) - offset_x) / dx + 0.5);
  particle->setPartid(id);
 
"""
)
 
coords = np.array([[0.21, 0.5, 0.0], [0.81, 0.5, 0.0]])
 
project.algorithm_step_initial_conditions.add_action_set(
    exaparty.api.input.InsertParticlesByCoordinates(
        particle_set=particle_set,
        initialisation_call=initialisation_snippet,
        coordinates=coords,
    )
)
 
 
 
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.partid, 1)
project.algorithm_step_plot.add_action_set(particle_plotter)
 
 
 
peano4_project = project.generate_Peano4_project(verbose=args.verbose)
peano4_project.output.makefile.add_CXX_flag("-D" + "HYDRO_DIMENSION=1")
 
peano4_project.constants.export_const_with_type(
    "GLOBAL_TIME_STEP_SIZE", str(args.timestep_size), "double"
)
 
# 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)
 
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)