Composition verification for external sources

.github/workflows/ci.yml

@@ -84,7 +84,7 @@ jobs:
     runs-on: ${{matrix.os}}
     strategy:
       matrix:
-        os: [ubuntu-24.04, macos-latest]
+        os: [ubuntu-22.04, macos-latest]
         build_type: [Release, Debug]
         enable_eb: [EB-OFF, EB-ON]
         include:
@@ -93,7 +93,7 @@ jobs:
             comp: llvm
             procs: $(sysctl -n hw.ncpu)
             ccache_cache: /Users/runner/Library/Caches/ccache
-          - os: ubuntu-24.04
+          - os: ubuntu-22.04
             install_deps: sudo apt-get update && sudo apt-get install mpich libmpich-dev
             comp: gnu
             procs: $(nproc)

Exec/RegTests/EB_ODEQty/CMakeLists.txt

@@ -1,6 +1,6 @@
 set(PELE_PHYSICS_EOS_MODEL Fuego)
-set(PELE_PHYSICS_CHEMISTRY_MODEL air)
-set(PELE_PHYSICS_TRANSPORT_MODEL Constant)
+set(PELE_PHYSICS_CHEMISTRY_MODEL Davis)
+set(PELE_PHYSICS_TRANSPORT_MODEL Simple)
 set(PELE_PHYSICS_ENABLE_SPRAY OFF)
 set(PELE_PHYSICS_SPRAY_FUEL_NUM 0)
 set(PELE_PHYSICS_ENABLE_SOOT OFF)

Exec/RegTests/EB_ODEQty/GNUmakefile

@@ -30,9 +30,9 @@ CEXE_sources += PeleLMeX_ProblemSpecificFunctions.cpp
 PELELM_NUM_ODE = 3
 
 # PelePhysics
-Chemistry_Model = air
+Chemistry_Model = Davis
 Eos_Model = Fuego
-Transport_Model = Constant
+Transport_Model = Simple
 
 PELE_HOME ?= ../../..
 include $(PELE_HOME)/Exec/Make.PeleLMeX

Exec/RegTests/EB_ODEQty/PeleLMeX_ProblemSpecificFunctions.cpp

@@ -25,7 +25,7 @@ set_ode_names(Vector<std::string>& a_ode_names)
 
 void
 problem_modify_ext_sources(
-  Real /*time*/,
+  Real time,
   Real /*dt*/,
   const MultiFab& state_old,
   const MultiFab& /*state_new*/,
@@ -45,10 +45,27 @@ problem_modify_ext_sources(
 
   ParallelFor(
     *ext_src, [=] AMREX_GPU_DEVICE(int box_no, int i, int j, int k) noexcept {
-      for (int n = 0; n < NUM_ODE; n++) {
-        Real B_n = state_old_arr[box_no](i, j, k, FIRSTODE + n);
-        Real src = prob_parm.ode_srcstrength * pow(10.0, n + 1) * B_n;
-        ext_src_arr[box_no](i, j, k, FIRSTODE + n) += src;
+      if (prob_parm.ode_qty_test) {
+        for (int n = 0; n < NUM_ODE; n++) {
+          // Source terms for ODE qty test
+
+          Real B_n = state_old_arr[box_no](i, j, k, FIRSTODE + n);
+          Real src = prob_parm.ode_srcstrength * pow(10.0, n + 1) * B_n;
+          ext_src_arr[box_no](i, j, k, FIRSTODE + n) += src;
+        }
+      }
+
+      // Source terms for composition test
+      if (prob_parm.composition_test) {
+        Real src = 0.0;
+
+        if (time >= prob_parm.extRhoYCO2_ts) {
+          Real CO2 = state_old_arr[box_no](i, j, k, FIRSTSPEC + CO2_ID);
+          src = prob_parm.extRhoYCO2 * CO2;
+        }
+
+        ext_src_arr[box_no](i, j, k, FIRSTSPEC + CO2_ID) += src;
+        ext_src_arr[box_no](i, j, k, DENSITY) += src;
       }
     });
   Gpu::streamSynchronize();

Exec/RegTests/EB_ODEQty/README.md

@@ -1,2 +1,4 @@
 ## EB\_ODEQty
-A 2D inflow/outflow setup with an optional EB cylinder in the middle of the flow. Demonstrates how to use ProblemSpecificFunctions and the ODE quantities.  The ODE quantities experience simple exponential decay that gets stiffer for each quantity.  Specifically, $\frac{\partial B_k}{\partial t} = \gamma \cdot 10^{k+1} B_k$, for $k = 0, 1,\dots,$ NUM_ODE and $\gamma < 0$.
+This RegTest contains two test cases that demonstrate the use of user defined external source terms:
+1. `prob_parm.ode_qty_test = true`: A 2D inflow/outflow case with an optional EB cylinder in the middle of the flow. This demonstrates how to use `ProblemSpecificFunctions` and the ODE quantities.  The ODE quantities experience simple exponential decay that gets stiffer for each quantity.  Specifically, $\frac{\partial B_k}{\partial t} = \gamma \cdot 10^{k+1} B_k$, for $k = 0, 1,\dots,$ NUM_ODE and $\gamma < 0$.
+2. `prob_parm.composition_test = true`:  This case tests the implementation of the user defined external sources for density, $\rho$, and the species $Y_m$.  The flow is set to zero with initial mass fractions $Y_{CO_2} = 0.15$, $Y_{N_2} = 0.25$ and $Y_{AR} = 0.6$.  After $t_s$ seconds, CO2 is added via an external source at a rate of $S_{ext,\rho Y_{CO_2}}$ kg/m$^3$/s.  Since CO2 is being added to the domain, density will increase at a rate of $S_{ext,\rho} = S_{ext,\rho Y_{CO_2}}$.  Results are compared to analytical solutions in `plot-comps.py` and tested in the CI via `test.py`.

Exec/RegTests/EB_ODEQty/composition-test-2d.inp

@@ -0,0 +1,58 @@
+#----------------------DOMAIN DEFINITION-------------------------
+geometry.is_periodic =  1 1                # For each dir, 0: non-perio, 1: periodic
+geometry.coord_sys   =  0                  # 0 => cart, 1 => RZ
+geometry.prob_lo     =  0 0
+geometry.prob_hi     =  0.01 0.01
+
+#--------------------------BC FLAGS------------------------------
+# Interior, Inflow, Outflow, Symmetry,
+# SlipWallAdiab, NoSlipWallAdiab, SlipWallIsotherm, NoSlipWallIsotherm
+peleLM.lo_bc = Interior  Interior
+peleLM.hi_bc = Interior  Interior
+
+#-------------------------AMR CONTROL----------------------------
+amr.n_cell          = 1 1          # Level 0 number of cells in each direction
+amr.v               = 0                # AMR verbose
+amr.max_level       = 0                # maximum level number allowed
+amr.n_error_buf     = 2 2 2 2          # number of buffer cells in error est
+amr.grid_eff        = 0.7              # what constitutes an efficient grid
+amr.blocking_factor = 1               # block factor in grid generation (min box size)
+amr.max_grid_size   = 64               # max box size
+
+#--------------------------- Problem ----------------------------
+prob.T_mean = 300.0            
+prob.P_mean = 101325.0    
+prob.meanFlowMag = 0.0
+prob.meanFlowDir = 1       
+prob.composition_test = true
+prob.Y_CO2_0  = 0.15  # Initial CO2 mass frac contribution
+prob.Y_N2_0   = 0.25 # Initial N2 mass frac contribution
+prob.Y_AR_0   = 0.6 # Initial AR mass frac contribution
+prob.extRhoYCO2 = 145.0  # Source strength for rhoYCO2
+prob.extRhoYCO2_ts = 0.003 # time when source starts
+
+#---------------------- Temporal CONTROL -------------------------
+peleLM.do_temporals = 1   
+peleLM.temporal_int = 1                 
+peleLM.do_extremas = 1                      
+
+#---------------------------TIME STEPPING------------------------
+amr.stop_time = 0.01
+amr.fixed_dt = 0.001
+
+#-------------------------PELE CONTROLS----------------------
+peleLM.v = 1
+peleLM.deltaT_verbose = 1
+peleLM.print_chi_convergence = 1
+peleLM.sdc_iterMax = 2
+peleLM.user_defined_ext_sources = true
+
+#---------------------------IO CONTROL---------------------------
+amr.check_file = "chk"
+amr.check_per = -1
+amr.plot_file = "plt"
+amr.plot_per = 0.1
+amr.derive_plot_vars = mass_fractions
+
+#------------------------- EB SETUP -----------------------------
+eb2.geom_type = "all_regular"

Exec/RegTests/EB_ODEQty/eb-odeqty-2d.inp

@@ -26,19 +26,13 @@ prob.T_mean = 300.0
 prob.P_mean = 101325.0
 prob.meanFlowMag = 4.255
 prob.meanFlowDir = 1
+prob.ode_qty_test = true
 prob.ode_IC = 10.0
 prob.ode_xy_lo = -0.01
 prob.ode_length = 0.01
 prob.ode_height = 0.02
 prob.ode_srcstrength = -10.0
 
-#-----------------INPUTS TO CONSTANT TRANSPORT ------------------
-transport.units                  = MKS
-transport.const_viscosity        = 0.0001
-transport.const_bulk_viscosity   = 0.0
-transport.const_conductivity     = 0.0
-transport.const_diffusivity      = 0.0
-
 #---------------------------TIME STEPPING------------------------
 #amr.max_step = 1
 amr.stop_time = 0.02
@@ -51,7 +45,6 @@ peleLM.v = 1
 peleLM.user_defined_ext_sources = 1
 
 #---------------------------IO CONTROL---------------------------
-#amr.restart = chk_07136
 amr.check_file = "chk_"
 amr.check_per = -1
 amr.plot_file = "plt_"

Exec/RegTests/EB_ODEQty/eb-odeqty-3d.inp

@@ -26,19 +26,13 @@ prob.T_mean = 300.0
 prob.P_mean = 101325.0
 prob.meanFlowMag = 4.255
 prob.meanFlowDir = 1
+prob.ode_qty_test = true
 prob.ode_IC = 10.0
 prob.ode_xy_lo = -0.01
 prob.ode_length = 0.01
 prob.ode_height = 0.02
 prob.ode_srcstrength = -10.0
 
-#-----------------INPUTS TO CONSTANT TRANSPORT ------------------
-transport.units                  = MKS
-transport.const_viscosity        = 0.0001
-transport.const_bulk_viscosity   = 0.0
-transport.const_conductivity     = 0.0
-transport.const_diffusivity      = 0.0
-
 #---------------------------TIME STEPPING------------------------
 #amr.max_step = 1
 amr.stop_time = 0.02
@@ -51,7 +45,6 @@ peleLM.v = 1
 peleLM.user_defined_ext_sources = 1
 
 #---------------------------IO CONTROL---------------------------
-#amr.restart = chk_07136
 amr.check_file = "chk_"
 amr.check_per = -1
 amr.plot_file = "plt_"

Exec/RegTests/EB_ODEQty/pelelmex_prob.H

@@ -27,8 +27,9 @@ pelelmex_initdata(
   auto eos = pele::physics::PhysicsType::eos();
 
   amrex::Real massfrac[NUM_SPECIES] = {0.0};
-  massfrac[O2_ID] = 0.233;
-  massfrac[N2_ID] = 0.767;
+  massfrac[CO2_ID] = prob_parm.Y_CO2_0;
+  massfrac[N2_ID] = prob_parm.Y_N2_0;
+  massfrac[AR_ID] = prob_parm.Y_AR_0;
 
   switch (prob_parm.meanFlowDir) {
   case 1:
@@ -90,29 +91,31 @@ pelelmex_initdata(
 #if NUM_ODE > 0
   amrex::Real ode_qty[NUM_ODE] = {0.0};
 
-  // Initial state for ode quantities
-  amrex::Real ode_qty_0[NUM_ODE] = {0.0};
-  for (int n = 0; n < NUM_ODE; n++) {
-    ode_qty_0[n] = (n + 1) * prob_parm.ode_IC;
-  }
-
-  // Current x,y,z locations
-  const amrex::Real* prob_lo = geomdata.ProbLo();
-  const amrex::Real* dx = geomdata.CellSize();
-  const amrex::Real x[AMREX_SPACEDIM] = {AMREX_D_DECL(
-    prob_lo[0] + static_cast<amrex::Real>(i + 0.5) * dx[0],
-    prob_lo[1] + static_cast<amrex::Real>(j + 0.5) * dx[1],
-    prob_lo[2] + static_cast<amrex::Real>(k + 0.5) * dx[2])};
+  if (prob_parm.ode_qty_test) {
+    // Initial state for ode quantities
+    amrex::Real ode_qty_0[NUM_ODE] = {0.0};
+    for (int n = 0; n < NUM_ODE; n++) {
+      ode_qty_0[n] = (n + 1) * prob_parm.ode_IC;
+    }
 
-  // Create NUM_ODE boxes from (x_strt,y_strt) to (x_end,y_end)
-  for (int n = 0; n < NUM_ODE; n++) {
-    amrex::Real x_strt = prob_parm.ode_xy_lo + n * 4. * prob_parm.ode_length;
-    amrex::Real x_end = x_strt + prob_parm.ode_length;
-    amrex::Real y_strt = prob_parm.ode_xy_lo;
-    amrex::Real y_end = y_strt + prob_parm.ode_height;
-    if ((x[0] > x_strt) && (x[1] > y_strt)) {
-      if ((x[0] < x_end) && (x[1] < y_end)) {
-        ode_qty[n] = ode_qty_0[n];
+    // Current x,y,z locations
+    const amrex::Real* prob_lo = geomdata.ProbLo();
+    const amrex::Real* dx = geomdata.CellSize();
+    const amrex::Real x[AMREX_SPACEDIM] = {AMREX_D_DECL(
+      prob_lo[0] + static_cast<amrex::Real>(i + 0.5) * dx[0],
+      prob_lo[1] + static_cast<amrex::Real>(j + 0.5) * dx[1],
+      prob_lo[2] + static_cast<amrex::Real>(k + 0.5) * dx[2])};
+
+    // Create NUM_ODE boxes from (x_strt,y_strt) to (x_end,y_end)
+    for (int n = 0; n < NUM_ODE; n++) {
+      amrex::Real x_strt = prob_parm.ode_xy_lo + n * 4. * prob_parm.ode_length;
+      amrex::Real x_end = x_strt + prob_parm.ode_length;
+      amrex::Real y_strt = prob_parm.ode_xy_lo;
+      amrex::Real y_end = y_strt + prob_parm.ode_height;
+      if ((x[0] > x_strt) && (x[1] > y_strt)) {
+        if ((x[0] < x_end) && (x[1] < y_end)) {
+          ode_qty[n] = ode_qty_0[n];
+        }
       }
     }
   }
@@ -139,6 +142,9 @@ bcnormal(
   pele::physics::PMF::PmfData::DataContainer const* /*pmf_data*/)
 {
   amrex::Real massfrac[NUM_SPECIES] = {0.0};
+  massfrac[CO2_ID] = prob_parm.Y_CO2_0;
+  massfrac[N2_ID] = prob_parm.Y_N2_0;
+  massfrac[AR_ID] = prob_parm.Y_AR_0;
 
   auto eos = pele::physics::PhysicsType::eos();
 
@@ -159,9 +165,6 @@ bcnormal(
     break;
   }
 
-  massfrac[O2_ID] = 0.333;
-  massfrac[N2_ID] = 0.667;
-
   s_ext[TEMP] = prob_parm.T_mean;
 
   amrex::Real rho_cgs, P_cgs, RhoH_temp;

Exec/RegTests/EB_ODEQty/pelelmex_prob.cpp

@@ -10,9 +10,16 @@ PeleLM::readProbParm() // NOLINT(readability-make-member-function-const)
   pp.query("P_mean", prob_parm->P_mean);
   pp.query("meanFlowDir", prob_parm->meanFlowDir);
   pp.query("meanFlowMag", prob_parm->meanFlowMag);
+  pp.query("ode_qty_test", prob_parm->ode_qty_test);
   pp.query("ode_IC", prob_parm->ode_IC);
   pp.query("ode_xy_lo", prob_parm->ode_xy_lo);
   pp.query("ode_length", prob_parm->ode_length);
   pp.query("ode_height", prob_parm->ode_height);
   pp.query("ode_srcstrength", prob_parm->ode_srcstrength);
+  pp.query("composition_test", prob_parm->composition_test);
+  pp.query("Y_CO2_0", prob_parm->Y_CO2_0);
+  pp.query("Y_N2_0", prob_parm->Y_N2_0);
+  pp.query("Y_AR_0", prob_parm->Y_AR_0);
+  pp.query("extRhoYCO2", prob_parm->extRhoYCO2);
+  pp.query("extRhoYCO2_ts", prob_parm->extRhoYCO2_ts);
 }

Exec/RegTests/EB_ODEQty/pelelmex_prob_parm.H

@@ -12,11 +12,18 @@ struct ProbParm : amrex::Gpu::Managed
   amrex::Real T_mean = 298.0_rt;
   amrex::Real P_mean = 101325.0_rt;
   amrex::Real meanFlowMag = 0.0;
-  amrex::Real ode_IC = 10.0;
-  amrex::Real ode_xy_lo = -0.01;
-  amrex::Real ode_length = 0.01;
-  amrex::Real ode_height = 0.020;
-  amrex::Real ode_srcstrength = -10.0;
   int meanFlowDir = 1;
+  bool ode_qty_test = false;
+  amrex::Real ode_IC = 0.0;
+  amrex::Real ode_xy_lo = 0.0;
+  amrex::Real ode_length = 0.0;
+  amrex::Real ode_height = 0.0;
+  amrex::Real ode_srcstrength = 0.0;
+  bool composition_test = false;
+  amrex::Real Y_CO2_0 = 0.0;
+  amrex::Real Y_N2_0 = 0.25;
+  amrex::Real Y_AR_0 = 0.75;
+  amrex::Real extRhoYCO2 = 0.0;
+  amrex::Real extRhoYCO2_ts = 0.0;
 };
 #endif