Working 0-1D Chemistry (among other things) (#653)

Co-authored-by: Dimitrios Adam <mitsos.gr@hotmail.gr>
Co-authored-by: Henry Le Berre <hberre3@gatech.edu>

.github/CODEOWNERS

@@ -13,3 +13,5 @@ CMakeLists.txt     @sbryngelson @henryleberre
 .vscode/           @sbryngelson @henryleberre
 .github/workflows/ @sbryngelson @henryleberre
 
+# Spencer H. Bryngelson & Ben Wilfong
+.typos.toml        @sbryngelson @wilfonba

.github/file-filter.yml

@@ -7,7 +7,7 @@ fortran_src: &fortran_src
 
 python_src: &python_src
   - '**/*.py'
-  - 'toolchain/requirements.txt'
+  - 'toolchain/pyproject.toml'
 
 cmakelist: &cmakelist
   - 'CMakeLists.txt'

CMakeLists.txt

@@ -133,6 +133,8 @@ if (CMAKE_Fortran_COMPILER_ID STREQUAL "GNU")
             -fimplicit-none
             #-ffpe-trap=invalid,zero,denormal,overflow
             -fsignaling-nans
+            -finit-real=snan
+            -finit-integer=-99999999
 	    )
     endif()
 

README.md

@@ -68,7 +68,7 @@ cd MFC
 ```
 and be sure MFC knows where to find Boost by appending to your dotfiles and sourcing them again
 ```shell
-echo -e 'export BOOST_INCLUDE=/opt/homebrew/' | tee -a ~/.bash_profile ~/.zshrc
+echo -e "export BOOST_INCLUDE='$(brew --prefix --installed boost)/include'" | tee -a ~/.bash_profile ~/.zshrc
 . ~/.bash_profile 2>/dev/null || . ~/.zshrc 2>/dev/null
 ! [ -z "${BOOST_INCLUDE+x}" ] && echo 'Environment is ready!' || echo 'Error: $BOOST_INCLUDE is unset. Please adjust the previous commands to fit with your environment.'
 ```

benchmarks/5eq_rk3_weno3_hllc/case.py

@@ -14,11 +14,13 @@
     description="This MFC case was created for the purposes of benchmarking MFC.",
     formatter_class=argparse.ArgumentDefaultsHelpFormatter)
 
-parser.add_argument("dict", type=str, metavar="DICT", help=argparse.SUPPRESS)
-parser.add_argument("gbpp", type=int, metavar="MEM", default=16, help="Adjusts the problem size per rank to fit into [MEM] GB of GPU memory per GPU.")
+parser.add_argument("--mfc", type=json.loads, default='{}', metavar="DICT",
+                    help="MFC's toolchain's internal state.")
+parser.add_argument("--gbpp", type=int, metavar="MEM", default=16,
+                    help="Adjusts the problem size per rank to fit into [MEM] GB of GPU memory per GPU.")
 
 ARGS = vars(parser.parse_args())
-DICT = json.loads(ARGS["dict"])
+DICT = ARGS["mfc"]
 
 size = 1 if DICT["gpu"] else 0
 

benchmarks/hypo_hll/case.py

@@ -12,11 +12,13 @@
     description="This MFC case was created for the purposes of benchmarking MFC.",
     formatter_class=argparse.ArgumentDefaultsHelpFormatter)
 
-parser.add_argument("dict", type=str, metavar="DICT", help=argparse.SUPPRESS)
-parser.add_argument("gbpp", type=int, metavar="MEM", default=16, help="Adjusts the problem size per rank to fit into [MEM] GB of GPU memory per GPU.")
+parser.add_argument("--mfc", type=json.loads, default='{}', metavar="DICT",
+                    help="MFC's toolchain's internal state.")
+parser.add_argument("--gbpp", type=int, metavar="MEM", default=16,
+                    help="Adjusts the problem size per rank to fit into [MEM] GB of GPU memory per GPU.")
 
 ARGS = vars(parser.parse_args())
-DICT = json.loads(ARGS["dict"])
+DICT = ARGS["mfc"]
 
 size = 1 if DICT["gpu"] else 0
 

benchmarks/ibm/case.py

@@ -11,11 +11,13 @@
     description="This MFC case was created for the purposes of benchmarking MFC.",
     formatter_class=argparse.ArgumentDefaultsHelpFormatter)
 
-parser.add_argument("dict", type=str, metavar="DICT", help=argparse.SUPPRESS)
-parser.add_argument("gbpp", type=int, metavar="MEM", default=16, help="Adjusts the problem size per rank to fit into [MEM] GB of GPU memory per GPU.")
+parser.add_argument("--mfc", type=json.loads, default='{}', metavar="DICT",
+                    help="MFC's toolchain's internal state.")
+parser.add_argument("--gbpp", type=int, metavar="MEM", default=16,
+                    help="Adjusts the problem size per rank to fit into [MEM] GB of GPU memory per GPU.")
 
 ARGS = vars(parser.parse_args())
-DICT = json.loads(ARGS["dict"])
+DICT = ARGS["mfc"]
 
 size = 1 if DICT["gpu"] else 0
 

benchmarks/viscous_weno5_sgb_acoustic/case.py

@@ -16,11 +16,13 @@
     description="This MFC case was created for the purposes of benchmarking MFC.",
     formatter_class=argparse.ArgumentDefaultsHelpFormatter)
 
-parser.add_argument("dict", type=str, metavar="DICT", help=argparse.SUPPRESS)
-parser.add_argument("gbpp", type=int, metavar="MEM", default=16, help="Adjusts the problem size per rank to fit into [MEM] GB of GPU memory per GPU.")
+parser.add_argument("--mfc", type=json.loads, default='{}', metavar="DICT",
+                    help="MFC's toolchain's internal state.")
+parser.add_argument("--gbpp", type=int, metavar="MEM", default=16,
+                    help="Adjusts the problem size per rank to fit into [MEM] GB of GPU memory per GPU.")
 
 ARGS = vars(parser.parse_args())
-DICT = json.loads(ARGS["dict"])
+DICT = ARGS["mfc"]
 
 size = 1 if DICT["gpu"] else 0
 

docs/documentation/case.md

@@ -34,24 +34,26 @@ Input files can accept command line arguments, forwarded by `mfc.sh run`.
 Consider this example from the `scaling` case:
 
 ```python
-import argparse
+import json, argparse
 
 parser = argparse.ArgumentParser(
     prog="scaling",
     description="Weak- and strong-scaling benchmark case.",
     formatter_class=argparse.ArgumentDefaultsHelpFormatter)
 
-parser.add_argument("dict", type=str, metavar="DICT")
-parser.add_argument("-s", "--scaling",  type=str, metavar="SCALING",  choices=["weak", "strong"], help="Whether weak- or strong-scaling is being exercised.")
+parser.add_argument("--mfc", type=json.loads, default='{}', metavar="DICT",
+                    help="MFC's toolchain's internal state.")
+parser.add_argument("-s", "--scaling", type=str, metavar="SCALING", choices=["weak", "strong"],
+                    help="Whether weak- or strong-scaling is being exercised.")
 
 # Your parsed arguments are here
 args = parser.parse_args()
 ```
 
-The first argument is always a JSON string representing `mfc.sh run`'s internal
-state.
+The `--mfc` argument is a JSON string representing `mfc.sh run`'s internal
+state, passed in when MFC runs your input file.
 It contains all the runtime information you might want from the build/run system.
-You can add as many additional arguments as you may need.
+You can add as many additional arguments and options as you may need.
 
 To run such a case, use the following format:
 

docs/documentation/getting-started.md

@@ -106,7 +106,7 @@ before configuring your environment:
 
 ```shell
 brew install coreutils python cmake fftw hdf5 gcc boost open-mpi
-echo -e 'export BOOST_INCLUDE=/opt/homebrew/' | tee -a ~/.bash_profile ~/.zshrc
+echo -e "export BOOST_INCLUDE='$(brew --prefix --installed boost)/include'" | tee -a ~/.bash_profile ~/.zshrc
 . ~/.bash_profile 2>/dev/null || . ~/.zshrc 2>/dev/null
 ! [ -z "${BOOST_INCLUDE+x}" ] && echo 'Environment is ready!' || echo 'Error: $BOOST_INCLUDE is unset. Please adjust the previous commands to fit with your environment.'
 ```

examples/1D_inert_shocktube/README.md

@@ -0,0 +1,12 @@
+# 1D Multi-Component Inert Shock Tube
+
+References:
+> P. J. Martínez Ferrer, R. Buttay, G. Lehnasch, and A. Mura, “A detailed verification procedure for compressible reactive multicomponent Navier–Stokes solvers”, Comput. & Fluids, vol. 89, pp. 88–110, Jan. 2014. Accessed: Oct. 13, 2024. [Online]. Available: https://doi.org/10.1016/j.compfluid.2013.10.014
+
+## Initial Condition
+
+![Initial Condition](initial.png)
+
+## Results
+
+![Results](result.png)

examples/1D_inert_shocktube/case.py

@@ -0,0 +1,124 @@
+#!/usr/bin/env python3
+
+# References:
+# + https://doi.org/10.1016/j.compfluid.2013.10.014: 4.3. Multi-component inert shock tube
+
+import json
+import argparse
+
+import cantera as ct
+
+parser = argparse.ArgumentParser(
+    prog="nD_inert_shocktube",
+    formatter_class=argparse.ArgumentDefaultsHelpFormatter)
+
+parser.add_argument("--mfc", type=json.loads, default='{}', metavar="DICT",
+                    help="MFC's toolchain's internal state.")
+parser.add_argument("--no-chem", dest='chemistry', default=True, action="store_false",
+                    help="Disable chemistry.")
+
+args = parser.parse_args()
+
+ctfile    = 'h2o2.yaml'
+sol_L     = ct.Solution(ctfile)
+sol_L.TPX =  400,  8000, 'H2:2,O2:1,AR:7'
+sol_R     = ct.Solution(ctfile)
+sol_R.TPX = 1200, 80000, 'H2:2,O2:1,AR:7'
+
+L    = 0.10
+Nx   = 400
+dx   = L / Nx
+dt   = 20e-8
+Tend = 40e-6
+
+NT         = int(Tend / dt)
+SAVE_COUNT = 200
+NS         = NT // SAVE_COUNT
+
+case = {
+    # Logistics ================================================================
+    'run_time_info'                : 'T',
+    # ==========================================================================
+
+    # Computational Domain Parameters ==========================================
+    'x_domain%beg'                 : -L/2,
+    'x_domain%end'                 : +L/2,
+    'm'                            : Nx,
+    'n'                            : 0,
+    'p'                            : 0,
+    'dt'                           : float(dt),
+    't_step_start'                 : 0,
+    't_step_stop'                  : NT,
+    't_step_save'                  : NS,
+    't_step_print'                 : NS,
+    'parallel_io'                  : 'F',
+    # ==========================================================================
+
+    # Simulation Algorithm Parameters ==========================================
+    'model_eqns'                   : 2,
+    'num_fluids'                   : 1,
+    'num_patches'                  : 2,
+    'mpp_lim'                      : 'F',
+    'mixture_err'                  : 'F',
+    'time_stepper'                 : 3,
+    'weno_order'                   : 5,
+    'weno_eps'                     : 1E-16,
+    'weno_avg'                     : 'F',
+    'mapped_weno'                  : 'T',
+    'mp_weno'                      : 'T',
+    'riemann_solver'               : 2,
+    'wave_speeds'                  : 2,
+    'avg_state'                    : 1,
+    'bc_x%beg'                     :-2,
+    'bc_x%end'                     :-3,
+    # ==========================================================================
+
+    # Chemistry ================================================================
+    'chemistry'                    : 'F' if not args.chemistry else 'T',
+    'chem_params%diffusion'        : 'F',
+    'chem_params%reactions'        : 'T',
+    # ==========================================================================
+
+    # Formatted Database Files Structure Parameters ============================
+    'format'                       : 1,
+    'precision'                    : 2,
+    'prim_vars_wrt'                : 'T',
+    # ==========================================================================
+
+    # ==========================================================================
+    'patch_icpp(1)%geometry'       : 1,
+    'patch_icpp(1)%x_centroid'     : -L/4,
+    'patch_icpp(1)%length_x'       : L/2,
+    'patch_icpp(1)%vel(1)'         : 0,
+    'patch_icpp(1)%pres'           : sol_L.P,
+    'patch_icpp(1)%alpha(1)'       : 1,
+    'patch_icpp(1)%alpha_rho(1)'   : sol_L.density,
+    # ==========================================================================
+
+    # ==========================================================================
+    'patch_icpp(2)%geometry'       : 1,
+    'patch_icpp(2)%x_centroid'     : L/4,
+    'patch_icpp(2)%length_x'       : L/2,
+    'patch_icpp(2)%vel(1)'         : 0,
+    'patch_icpp(2)%pres'           : sol_R.P,
+    'patch_icpp(2)%alpha(1)'       : 1,
+    'patch_icpp(2)%alpha_rho(1)'   : sol_R.density,
+    # ==========================================================================
+
+    # Fluids Physical Parameters ===============================================
+    'fluid_pp(1)%gamma'            : 1.0E+00/(1.55E+00-1.0E+00),
+    'fluid_pp(1)%pi_inf'           : 0,
+    # ==========================================================================
+
+    # Chemistry ================================================================
+    'cantera_file'                 : ctfile,
+    # ==========================================================================
+}
+
+if args.chemistry:
+    for i in range(len(sol_L.Y)):
+        case[f'patch_icpp(1)%Y({i+1})'] = sol_L.Y[i]
+        case[f'patch_icpp(2)%Y({i+1})'] = sol_R.Y[i]
+
+if __name__ == '__main__':
+    print(json.dumps(case))

examples/1D_inert_shocktube/viz.py

@@ -0,0 +1,65 @@
+import mfc.viz
+
+import os
+
+import subprocess
+import seaborn as sns
+import matplotlib.pyplot as plt
+from tqdm import tqdm
+
+from case import sol_L as sol
+
+case = mfc.viz.Case(".")
+
+os.makedirs("viz", exist_ok=True)
+
+#sns.set_theme(style=mfc.viz.generate_cpg_style())
+
+Y_VARS = ["H2", "O2", "H2O", "N2"]
+
+variables = [
+    ("rho", "prim.1"),
+    ("u_x", "prim.2"),
+    ("p",   "prim.3"),
+    ("E",   "cons.3"),
+    *[(f"Y_{name}", f"prim.{5 + sol.species_index(name)}") for name in Y_VARS],
+    ("T",   "prim.15"),
+]
+
+for variable in tqdm(variables, desc="Loading Variables"):
+    case.load_variable(*variable)
+
+for step in tqdm(case.get_timesteps(), desc="Rendering Frames"):
+    fig, axes = plt.subplots(2, 3, figsize=(16, 9))
+
+    def pad_ylim(ylim, pad=0.1):
+        return (ylim[0] - pad*(ylim[1] - ylim[0]), ylim[1] + pad*(ylim[1] - ylim[0]))
+
+    case.plot_step(step, "rho", ax=axes[0, 0])
+    axes[0, 0].set_ylim(*pad_ylim(case.get_minmax_time("rho")))
+    axes[0, 0].set_ylabel("$\\rho$")
+    case.plot_step(step, "u_x", ax=axes[0, 1])
+    axes[0, 1].set_ylim(*pad_ylim(case.get_minmax_time("u_x")))
+    axes[0, 1].set_ylabel("$u_x$")
+    case.plot_step(step, "p",   ax=axes[1, 0])
+    axes[1, 0].set_ylim(*pad_ylim(case.get_minmax_time("p")))
+    axes[1, 0].set_ylabel("$p$")
+    for y in Y_VARS:
+        case.plot_step(step, f"Y_{y}", ax=axes[1, 1], label=y)
+    axes[1, 1].set_ylim(0, 1.1*max(case.get_minmax_time(f"Y_{y}")[1] for y in Y_VARS))
+    axes[1, 1].set_ylabel("$Y_k$")
+    case.plot_step(step, "T",   ax=axes[1, 2])
+    axes[1, 2].set_ylim(*pad_ylim(case.get_minmax_time("T")))
+    axes[1, 2].set_ylabel("$T$")
+    case.plot_step(step, "E",   ax=axes[0, 2])
+    axes[0, 2].set_ylim(*pad_ylim(case.get_minmax_time("E")))
+    axes[0, 2].set_ylabel("$E$")
+
+    plt.tight_layout()
+    plt.savefig(f"viz/{step:06d}.png")
+    plt.close()
+
+subprocess.run([
+    "ffmpeg", "-y", "-framerate", "60", "-pattern_type", "glob", "-i",
+    "viz/*.png", "-c:v", "libx264", "-pix_fmt", "yuv420p", "viz.mp4"
+])

examples/1D_reactive_shocktube/README.md

@@ -0,0 +1,14 @@
+# 1D Multi-Component Reactive Shock Tube
+
+References:
+> P. J. Martínez Ferrer, R. Buttay, G. Lehnasch, and A. Mura, “A detailed verification procedure for compressible reactive multicomponent Navier–Stokes solvers”, Comput. & Fluids, vol. 89, pp. 88–110, Jan. 2014. Accessed: Oct. 13, 2024. [Online]. Available: https://doi.org/10.1016/j.compfluid.2013.10.014
+
+> H. Chen, C. Si, Y. Wu, H. Hu, and Y. Zhu, “Numerical investigation of the effect of equivalence ratio on the propagation characteristics and performance of rotating detonation engine”, Int. J. Hydrogen Energy, Mar. 2023. Accessed: Oct. 13, 2024. [Online]. Available: https://doi.org/10.1016/j.ijhydene.2023.03.190
+
+## Initial Condition
+
+![Initial Condition](initial.png)
+
+## Results
+
+![Results](result.png)