Add CheMFC test cases

MFlowCode · Nov 2, 2024 · 06a757a · 06a757a
1 parent d811942
commit 06a757a
Show file tree

Hide file tree

Showing 35 changed files with 918 additions and 118 deletions.
diff --git a/benchmarks/5eq_rk3_weno3_hllc/case.py b/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("--mfc", 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["mfc"])
+DICT = ARGS["mfc"]
 
 size = 1 if DICT["gpu"] else 0
 

diff --git a/benchmarks/hypo_hll/case.py b/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("--mfc", 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["mfc"])
+DICT = ARGS["mfc"]
 
 size = 1 if DICT["gpu"] else 0
 

diff --git a/benchmarks/ibm/case.py b/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("--mfc", 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["mfc"])
+DICT = ARGS["mfc"]
 
 size = 1 if DICT["gpu"] else 0
 

diff --git a/benchmarks/viscous_weno5_sgb_acoustic/case.py b/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("--mfc", 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["mfc"])
+DICT = ARGS["mfc"]
 
 size = 1 if DICT["gpu"] else 0
 

diff --git a/docs/documentation/case.md b/docs/documentation/case.md
@@ -34,16 +34,17 @@ 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("--mfc",           type=str, metavar="DICT", default='{}')
+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.")
+                    help="Whether weak- or strong-scaling is being exercised.")
 
 # Your parsed arguments are here
 args = parser.parse_args()

diff --git a/examples/1D_inert_shocktube/case.py b/examples/1D_inert_shocktube/case.py
@@ -4,8 +4,21 @@
 # + 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'
@@ -18,7 +31,6 @@
 dt   = 5e-9
 Tend = 40e-6
 
-chemistry  = True
 NT         = int(Tend / dt)
 SAVE_COUNT = 200
 NS         = NT // SAVE_COUNT
@@ -62,8 +74,7 @@
     # ==========================================================================
 
     # Chemistry ================================================================
-    'chemistry'                    : 'F' if not chemistry else 'T',
-    'chem_params%advection'        : 'T',
+    'chemistry'                    : 'F' if not args.chemistry else 'T',
     'chem_params%diffusion'        : 'F',
     'chem_params%reactions'        : 'T',
     # ==========================================================================
@@ -104,7 +115,7 @@
     # ==========================================================================
 }
 
-if chemistry:
+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]

diff --git a/examples/1D_reactive_shocktube/case.py b/examples/1D_reactive_shocktube/case.py
@@ -4,9 +4,21 @@
 # + https://doi.org/10.1016/j.ijhydene.2023.03.190:  Verification of numerical method
 # + https://doi.org/10.1016/j.compfluid.2013.10.014: 4.7. Multi-species reactive shock tube
 
-import json
+import json, argparse
 import cantera as ct
 
+parser = argparse.ArgumentParser(
+    prog="1D_reactive_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='chem', default=True, action="store_false",
+                                 help="Disable chemistry.")
+parser.add_argument("--scale",   type=float,  default=1,    help="Scale.")
+
+args = parser.parse_args()
+
 ctfile    = 'h2o2.yaml'
 sol_L     = ct.Solution(ctfile)
 sol_L.DPX =   0.072,  7173, 'H2:2,O2:1,AR:7'
@@ -18,7 +30,7 @@
 u_r = -487.34
 
 L  = 0.12
-Nx = 800
+Nx = 800 * args.scale
 dx = L/Nx
 dt = dx/abs(u_r)*0.01
 Tend=230e-6
@@ -27,8 +39,6 @@
 SAVE_COUNT=100
 NS=NT//SAVE_COUNT
 
-chemistry = True
-
 case = {
     # Logistics ================================================================
     'run_time_info'                : 'T',
@@ -45,7 +55,7 @@
     't_step_stop'                  : NT,
     't_step_save'                  : NS,
     't_step_print'                 : NS,
-    'parallel_io'                  : 'F',
+    'parallel_io'                  : 'T' if args.mfc.get("mpi", True) else 'F',
 
     # Simulation Algorithm Parameters ==========================================
     'model_eqns'                   : 2,
@@ -66,8 +76,7 @@
     'bc_x%end'                     :-3,
 
     # Chemistry ================================================================
-    'chemistry'                    : 'F' if not chemistry else 'T',
-    'chem_params%advection'        : 'T',
+    'chemistry'                    : 'F' if not args.chemistry else 'T',
     'chem_params%diffusion'        : 'F',
     'chem_params%reactions'        : 'T',
     'cantera_file'                 : ctfile,
@@ -77,6 +86,7 @@
     'format'                       : 1,
     'precision'                    : 2,
     'prim_vars_wrt'                : 'T',
+    'chem_wrt_T'                   : 'T',
     # ==========================================================================
 
     # ==========================================================================
@@ -105,8 +115,9 @@
     # ==========================================================================
 }
 
-if chemistry:
+if args.chemistry:
     for i in range(len(sol_L.Y)):
+        case[f'chem_wrt_Y({i + 1})']    = 'T'
         case[f'patch_icpp(1)%Y({i+1})'] = sol_L.Y[i]
         case[f'patch_icpp(2)%Y({i+1})'] = sol_R.Y[i]
 

diff --git a/examples/nD_perfect_reactor/analyze.py b/examples/nD_perfect_reactor/analyze.py
@@ -4,7 +4,7 @@
 import matplotlib.pyplot as plt
 
 import mfc.viz
-from case import dt, NS, Tend, SAVE_COUNT, sol
+from case import dt, Tend, SAVE_COUNT, sol
 
 
 case = mfc.viz.Case(".", dt)

diff --git a/examples/nD_perfect_reactor/case.py b/examples/nD_perfect_reactor/case.py
@@ -9,14 +9,15 @@
 import cantera as ct
 
 parser = argparse.ArgumentParser(
-    prog="nD_Reactor",
+    prog="nD_perfect_reactor",
     formatter_class=argparse.ArgumentDefaultsHelpFormatter)
 
-parser.add_argument("--mfc",     type=str,    default='{}', metavar="DICT")
-parser.add_argument("--no-chem", dest='chem', default=True, action="store_false",
-                                 help="Disable chemistry.")
-parser.add_argument("--scale",   type=float,  default=1,    help="Scale.")
-parser.add_argument("--ndim",    type=int,    default=1,    help="Number of dimensions.")
+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.")
+parser.add_argument("--scale",   type=float,       default=1,    help="Scale.")
+parser.add_argument("--ndim",    type=int,         default=1,    help="Number of dimensions.")
 
 args = parser.parse_args()
 
@@ -28,10 +29,10 @@
 Nx   = 25 * args.scale
 Tend = 1e-4
 s    = 1e-2
-dt   = 1e-9
+dt   = 1e-7
 
 NT         = int(Tend / dt)
-SAVE_COUNT = 2000
+SAVE_COUNT = 20
 NS         = NT // SAVE_COUNT
 
 case = {
@@ -54,7 +55,7 @@
     't_step_stop'                  : NT,
     't_step_save'                  : NS,
     't_step_print'                 : NS,
-    'parallel_io'                  : 'T' if args.ndim > 1 else 'F',
+    'parallel_io'                  : 'T' if args.ndim > 1 and args.mfc.get("mpi", True) else 'F',
 
     # Simulation Algorithm Parameters ==========================================
     'model_eqns'                   : 2,
@@ -82,6 +83,7 @@
     'format'                       : 1,
     'precision'                    : 2,
     'prim_vars_wrt'                : 'T',
+    'chem_wrt_T'                   : 'T',
     # ==========================================================================
 
     # Patch 1 ==================================================================
@@ -106,18 +108,18 @@
     # ==========================================================================
 }
 
-if args.chem:
+if args.chemistry:
     case.update({
         # Chemistry ============================================================
         'chemistry'             : 'T',
-        'chem_params%advection' : 'F',
         'chem_params%diffusion' : 'F',
         'chem_params%reactions' : 'T',
         'cantera_file'          : ctfile,
         # ======================================================================
     })
 
     for i in range(len(sol.Y)):
+        case[f'chem_wrt_Y({i + 1})']    = 'T'
         case[f'patch_icpp(1)%Y({i+1})'] = sol.Y[i]
 
 case = remove_higher_dimensional_keys(case, args.ndim)

diff --git a/examples/nD_perfect_reactor/export.py b/examples/nD_perfect_reactor/export.py
@@ -0,0 +1,76 @@
+import csv
+import cantera as ct
+from tqdm import tqdm
+
+import mfc.viz
+from case import dt, NS, Tend, SAVE_COUNT, sol
+
+
+case = mfc.viz.Case(".", dt)
+
+for name in tqdm(sol.species_names, desc="Loading Variables"):
+    case.load_variable(f"Y_{name}", f"prim.{5 + sol.species_index(name)}")
+case.load_variable("rho", "prim.1")
+
+time_save = Tend/SAVE_COUNT
+
+oh_idx = sol.species_index("OH")
+def generate_ct_saves() -> tuple:
+    reactor         = ct.IdealGasReactor(sol)
+    reactor_network = ct.ReactorNet([reactor])
+
+    ct_time = 0.0
+    ct_ts, ct_Ys, ct_rhos = [0.0], [reactor.thermo.Y], [reactor.thermo.density]
+
+    while ct_time < Tend:
+        reactor_network.advance(ct_time + time_save)
+        ct_time += time_save
+        ct_ts.append(ct_time)
+        ct_Ys.append(reactor.thermo.Y)
+        ct_rhos.append(reactor.thermo.density)
+
+    return ct_ts, ct_Ys, ct_rhos
+
+ct_ts, ct_Ys, ct_rhos = generate_ct_saves()
+
+with open("mfc.csv", "w") as f:
+    writer = csv.writer(f)
+    keys = ["t"] + list(set(case.get_data()[0].keys()) - set(["x"]))
+    writer.writerow(keys)
+    for i, t_step in enumerate(sorted(case.get_timesteps())):
+        t = t_step * dt
+        row = [t] + [case.get_data()[t_step][key][0] for key in keys[1:]]
+        writer.writerow(row)
+
+with open("cantera.csv", "w") as f:
+    writer = csv.writer(f)
+    keys = ["t"] + [f"Y_{_}" for _ in list(sol.species_names)] + ["rho"]
+    writer.writerow(keys)
+    for step in range(len(ct_ts)):
+        row = [ct_ts[step]] + [ ct_Ys[step][i] for i in range(len(sol.species_names)) ] + [ct_rhos[step]]
+        print([ct_ts[step]], row)
+        writer.writerow(row)
+
+def find_induction_time(ts: list, Ys: list, rhos: list) -> float:
+    for t, y, rho in zip(ts, Ys, rhos):
+        if (y * rho / sol.molecular_weights[oh_idx]) >= 1e-6:
+            return t
+
+    return None
+
+skinner_induction_time = 0.052e-3
+ct_induction_time      = find_induction_time(
+    ct_ts,
+    [y[oh_idx] for y   in ct_Ys],
+    [rho       for rho in ct_rhos]
+)
+mfc_induction_time     = find_induction_time(
+    sorted(case.get_timestamps()),
+    [ case.get_data()[step]["Y_OH"][0] for step in sorted(case.get_timesteps()) ],
+    [ case.get_data()[step]["rho"][0]  for step in sorted(case.get_timesteps()) ]
+)
+
+print("Induction Times ([OH] >= 1e-6):")
+print(f" + Skinner et al.: {skinner_induction_time} s")
+print(f" + Cantera:        {ct_induction_time} s")
+print(f" + (Che)MFC:       {mfc_induction_time} s")