Regression tests for solver examples

examples/H3LAPD/.gitignore

@@ -1,4 +1,3 @@
 # Ignore xml files generated from .geos
-cuboid*.xml
 low_res.xml
 full_res.xml

examples/H3LAPD/README.md

@@ -45,22 +45,22 @@ Then run the example with
 
 This script expects to find mpirun on the path and executes with four MPI ranks by default. It looks for a solver executable in the most recently modified spack-build* directory, but this can be overridden using the '-b' option.
 
-### 2Din3D-hw
+### 2Din3DHW
 
 Solves equations (1) and (2), as in the previous example, but also enables a system of neutral particles that are coupled to the fluid solver. Particles deposit density into the (plasma) fluid via ionization.
 
 Generate the mesh with
 
-    ./scripts/geo_to_xml.sh examples/H3LAPD/2Din3D-hw/cuboid_periodic_8x8x16.geo -x 1,2 -y 3,4 -z 5,6 -o cuboid.xml
+    ./scripts/geo_to_xml.sh examples/H3LAPD/2Din3DHW/cuboid_periodic_8x8x16.geo -x 1,2 -y 3,4 -z 5,6 -o cuboid.xml
 
 Then run the example with
 
-    ./scripts/run_eg.sh H3LAPD 2Din3D-hw
+    ./scripts/run_eg.sh H3LAPD 2Din3DHW
 
 This script expects to find mpirun on the path and executes with four MPI ranks by default. It looks for a solver executable in the most recently modified spack-build* directory, but this can be overridden using the '-b' option.
 
 ## Diagnostics
-For the '2Din3DHW' equation system (used in the `2Din3D-hw` and `2Din3D-hw_fluid-only` examples), the solver can be made to output the total fluid energy ($E$) and enstrophy ($W$), which are defined as:  
+For the '2Din3DHW' equation system (used in the `2Din3DHW` and `2Din3D-hw_fluid-only` examples), the solver can be made to output the total fluid energy ($E$) and enstrophy ($W$), which are defined as:  
 
 $$
 \begin{align}

examples/SimpleSOL/2DWithParticles/2DWithParticles_config.xml

@@ -15,7 +15,23 @@
             <P> pInf           = 1.0         </P>
             <P> rhoInf         = 1.0         </P>
             <P> uInf           = 1.0         </P>
-
+            <P> num_particle_steps_per_fluid_step = 1 </P>
+             <!--
+                Note that, for the SOLWithParticles equation system, num_particles_total sets the
+                number of particles added to the simulation per timestep per MPI rank, rather than
+                the total number of particles that will be added.
+            -->
+            <P> num_particles_total = 100 </P>
+            <P> num_particles_per_cell = -1 </P>
+            <P> particle_output_freq = 4 </P>
+            <P> particle_thermal_velocity = 1.0 </P>
+            <P> particle_number_density = 3e18 </P>
+            <P> particle_source_region_count = 2 </P>
+            <P> particle_source_region_offset = 0.01 </P>
+            <P> particle_source_region_gaussian_width = 0.00001 </P>
+            <P> particle_source_lines_per_gaussian = 7 </P>
+            <P> unrotated_x_max = 110.0 </P>
+            <P> unrotated_y_max = 1.0 </P>
         </PARAMETERS>
         <SOLVERINFO>
             <I PROPERTY="EQTYPE"                VALUE="SOLWithParticles"     />
@@ -99,7 +115,7 @@
 
     <PARTICLES>
         <INFO>
-            <I PROPERTY="PARTTYPE" VALUE="NeutralParticleSystem"/>
+            <I PROPERTY="PARTTYPE" VALUE="SOLParticleSystem"/>
         </INFO>
 
         <PARAMETERS>
@@ -121,29 +137,5 @@
             <P> unrotated_x_max = 110.0 </P>
             <P> unrotated_y_max = 1.0 </P>
         </PARAMETERS>
-
-        <SPECIES>
-            <S ID="0" NAME="Neutral">
-                <P PROPERTY="Mass" VALUE="1.0"/>
-                <P PROPERTY="Charge" VALUE="0"/>
-                <P PROPERTY="Number" VALUE="100"/>
-            </S>
-        </SPECIES>
-
-        <BOUNDARYINTERACTION>
-            <REGION REF="0">
-                <C SPECIES="Neutral"  VALUE="0.0" />
-            </REGION>
-            <REGION REF="1">
-                <C SPECIES="Neutral"  VALUE="0.0" />
-            </REGION>
-	    <REGION REF="2">
-                <P SPECIES="Neutral"  VALUE="[3]" />
-            </REGION>
-            <REGION REF="3">
-                <P SPECIES="Neutral"  VALUE="[2]" />
-            </REGION>
-        </BOUNDARYINTERACTION>
-
     </PARTICLES>
 </NEKTAR>

python/regression_data_generator/gen_eg_regression_data.py

@@ -0,0 +1,159 @@
+import argparse
+import glob
+import h5py
+import math
+import os.path
+import sys
+
+from NekPy.FieldUtils import Field, InputModule, ProcessModule
+
+examples_root = os.path.normpath(sys.path[0] + "/../../examples")
+eg_reg_tests_dir = os.path.normpath(sys.path[0] + "/../../test/regression/examples")
+
+
+def determine_config_paths(solver_name, eg_name):
+    cmd_template_path = os.path.join(
+        examples_root, solver_name, eg_name, "run_cmd_template.txt"
+    )
+    config_fpaths = []
+    if os.path.exists(cmd_template_path):
+        with open(cmd_template_path, "r") as cmd_template:
+            line = cmd_template.readline()
+            config_fpaths = [
+                os.path.join(examples_root, solver_name, eg_name, elmt)
+                for elmt in line.split()
+                if elmt.endswith(".xml")
+            ]
+
+    assert (
+        len(config_fpaths) > 0
+    ), f"Unable to determine xml file paths from command template at {cmd_template}"
+    return config_fpaths
+
+
+def get_nsteps(session, tol=1e-12):
+    """
+    Extract number of steps from the nektar session.
+    """
+    nsteps_str = "NumSteps"
+    if session.DefinesParameter(nsteps_str):
+        return int(session.GetParameter(nsteps_str))
+    else:
+        exit(f"Expected {nsteps_str} to be defined in the example's session file")
+
+
+def read_nektar_fld_and_gen_pts(example_runs_root, solver_name, eg_name):
+    """
+    1. Look for xml files and a .fld file in run_dir
+    2. Use nektar's FieldConvert modules to generate equispaced points from the .fld data
+    3. Write the points
+    """
+    config_fpaths = determine_config_paths(solver_name, eg_name)
+    run_dir = os.path.join(example_runs_root, solver_name, eg_name)
+    fld_fpath = glob.glob(os.path.join(run_dir, "*.fld"))
+    assert len(fld_fpath) > 0, f"No .fld file in {run_dir}"
+    assert len(fld_fpath) == 1, f"Found multiple .fld files in {run_dir}"
+
+    # Init Field object
+    field = Field(config_fpaths)
+
+    # Read config
+    InputModule.Create("xml", field, *config_fpaths).Run()
+
+    # Read fld file
+    InputModule.Create("fld", field, fld_fpath[0]).Run()
+
+    # Compute equi-spaced points
+    ProcessModule.Create("equispacedoutput", field).Run()
+
+    # Return points in a dict to simplify output later
+    ndims = field.graph.GetSpaceDimension()
+    nsteps = get_nsteps(field.session)
+    return nsteps, {
+        fld_name: field.GetPts(ndims + fld_idx)
+        for fld_idx, fld_name in enumerate(field.session.GetVariables())
+    }
+
+
+def user_confirms(msg: str) -> bool:
+    opts = ("y", "n")
+    opt_str = " or ".join(opts)
+    while (
+        answer := input(f"{msg} (Choose {opt_str})").strip().lower()[:1]
+    ) not in opts:
+        print(f"Invalid input; choose {opt_str}")
+
+    return answer == opts[0]
+
+
+def data_dir(solver_name: str):
+    return os.path.join(eg_reg_tests_dir, solver_name)
+
+
+def reg_data_fname(eg_name: str):
+    return f"{eg_name}.regression_data.h5"
+
+
+def gen_eg_regression_data(
+    solver_name: str, eg_name: str, attrs={}, example_runs_root=None
+):
+    eg_reg_data_dir = data_dir(solver_name)
+    os.makedirs(eg_reg_data_dir, exist_ok=True)
+
+    if example_runs_root is None:
+        example_runs_root = os.path.normpath(sys.path[0] + "/../../runs")
+    nsteps, fld_data = read_nektar_fld_and_gen_pts(
+        example_runs_root, solver_name, eg_name
+    )
+    attrs["nsteps"] = nsteps
+
+    pth = os.path.join(eg_reg_data_dir, reg_data_fname(eg_name))
+    if os.path.exists(pth):
+        if not user_confirms(f"Overwrite file at {pth}?"):
+            print("Aborted.")
+            return
+    with h5py.File(pth, "w") as fh:
+        for fld_name, fld_vals in fld_data.items():
+            fh.create_dataset(fld_name, data=fld_vals)
+        fh.attrs.update(attrs)
+    print(f"Wrote regression data to {pth}")
+
+
+def read_regression_data(solver_name: str, eg_name: str, nvals: int):
+    pth = os.path.join(data_dir(solver_name), reg_data_fname(eg_name))
+    if not os.path.exists(pth):
+        print(f"No regression data found at {pth}")
+        return
+
+    print(f"Reading regression data at {pth}")
+    with h5py.File(pth, "r") as data:
+        print(f" nsteps = {data.attrs['nsteps']}")
+        print(f" Field data:")
+        for fld_name, fld_data in data.items():
+            print(f"  {fld_name}: {fld_data}")
+            print(f"    {','.join([str(v) for v in [fld_data[:nvals]]])}")
+
+
+if __name__ == "__main__":
+    parser = argparse.ArgumentParser(
+        description=f"Extract data from a Nektar++ .fld file (generated by running one of the solver examples) and write it to an HDF5 file in {eg_reg_tests_dir}/[solver_name]/[example_name]",
+    )
+    parser.add_argument("solver_name")
+    parser.add_argument("example_name")
+    parser.add_argument(
+        "-r",
+        "--read",
+        action="store_true",
+        help="read regression data and output some of its properties",
+    )
+    parser.add_argument(
+        "-n",
+        "--nvals",
+        type=int,
+        help="number of values to print from stored data; ignored unless -r/--read was supplied",
+    )
+    args = parser.parse_args(args=None if sys.argv[1:] else ["--help"])
+    if args.read:
+        read_regression_data(args.solver_name, args.example_name, args.nvals)
+    else:
+        gen_eg_regression_data(args.solver_name, args.example_name)