Latest updates (#1)

* DivU constraint with arbitrary sources (AMReX-Combustion#428)

* Included additional terms in divu constraint to account for external forces on density.

* Updated Docs to include external source terms in model equations and divergence constraint.

* Formatting

* fix a bug in divu calc

* Update deriviation of divu in Docs

* Restored code to development

* Reverted text in line 168

* Minor fix in Docs

* Comment in divu about extRho

---------

Co-authored-by: Bruce Perry <[email protected]>

* "ODE" quantities and user defined external sources (AMReX-Combustion#427)

* Initial commit.  Needs debugging for source terms.

* Added ode_var to typical values. Still debugging initialization and source.

* Minor edit in set_ode_names function.

* Fixed misalignment in plotting. IC's appear to be correct now for multiple ode vars.  Still working on source terms.

* Fix indexing.

* Updated problem_modify_ext_sources.  Works with local case definition.

* Comment arguments in problem_modify_ext_sources in Source code.

* Added compilation variable USE_MODIFIED_SOURCES to avoid unnecessary loops through levels.

* Tested plotting for Efield and Spray.  Tested all args of ProblemSpecificFunctions.

* Added ode_srcstrength to test case, fixed minor bugs.

* Changed flag for user defined external sources.  Added future flag for including external sources in SDC loop.

* Set peleLM.user_defined_ext_sources to true in input file.

* Added Euler step for advancing the ODE quantities in time.

* Changed name of advanceODEQty() to predictODEQty() for consistency with SDC iterations.

* Updated input.2d for validation

* Added .cpp file for ProblemSpecificFunctions to avoid need for static functions.

* Added convergence test for ODE Qty.

* Updated the test case to include 3 odes with increasing stiffness.

* Renamed case

* Combined all external forces in single function call to clean up PeleLMeX_Advance.cpp

* Update README for case.

* Update README for case.

* Update README.md

* Removed NUM_AUX for now.

* Cleaning up

* Documentation for ode quantities.

* Clang-format

* Remove validation tests from case directory.

* Minor fix for soot and spray

* Spelling and unused variable handling.

* Unused variable handling

* Formatting

* Removed unwanted print statements, renamed some files/vars for consistency.

* Added case to CI

* Add USE_EB to EB_ODEQty/CmakeLists.txt

* Add EB_ODEQty to RegTests/CMakeLists.txt

* CI fix

* More CI

* Add ODEQty.cpp to cmake

* More CI...

* Only define predictODEQty if NUM_ODE > 0

* Updated CMake to handle user defined EB's and ProblemSpecificFunctions

* Now passes in MultiFabs for state_old and state_new. This provides access to state_old.Factory() etc.

* Update PeleLMeX_ProblemSpecificFunctions.cpp

* Removed lev argument in problem_modify_ext_sources

* Update comment

---------

Co-authored-by: Bruce Perry <[email protected]>

* pin intel-oneapi version in CI (AMReX-Combustion#434)

* Bump Submodules/AMReX-Hydro from `3ab9864` to `3788361` (AMReX-Combustion#433)

Bumps [Submodules/AMReX-Hydro](https://github.com/AMReX-Fluids/AMReX-Hydro) from `3ab9864` to `3788361`.
- [Commits](AMReX-Fluids/AMReX-Hydro@3ab9864...3788361)

---
updated-dependencies:
- dependency-name: Submodules/AMReX-Hydro
  dependency-type: direct:production
...

Signed-off-by: dependabot[bot] <[email protected]>
Co-authored-by: dependabot[bot] <49699333+dependabot[bot]@users.noreply.github.com>
Co-authored-by: Bruce Perry <[email protected]>

* Add full manifold model capability (AMReX-Combustion#415)

* basics of having eos_parm in PeleLMeX

* add eosparm to PeleLMeX_K functions

* more eosparm: now everywhere except Efield & BCs

* eosparm in BC stuff

* re-order functions for template

* remove manifold-specfic changes, will go in later PR

* gnu make changes for manifold

* create manifold version of adjust fluxes function

* Add abort for Closed Chamber + Manifold

* Better setup chackes for Manifold + wbar, closed chamber, mixfrac/prog

* advection of rho corrections for manifold models

* manifold-specific transport, derives, divu

* remove divide by zero

* go to templated EOS functions to use cellData

* fix oops in pelephysics

* fix compile errors for Manifold derives

* fix compile condition for transparm initialization

* fix missing RY2RRinvY conversion

* use PelePhysics that free manfunc shared pointer

* clang-formatting

* fux derives for manifold

* add flamesheet test case for manifold

* remove extraneous changes in file

* Update Source/PeleLMeX_Setup.cpp

Co-authored-by: Marc T. Henry de Frahan <[email protected]>

* address Marc HdF's comments with minor formatting

* remove copied line

* update PP for merged RY2R stuff

---------

Co-authored-by: Marc T. Henry de Frahan <[email protected]>

* Update PelePhysics and add path to new utility  (AMReX-Combustion#435)

* Update PelePhysics and add new utility to Make.PeleLMeX

* Add PelePhysics utilities to CMake

* Test PelePhysics:tidy against clang-tidy.

* Update PelePhysics submodule

* Bump Submodules/AMReX-Hydro from `3ab9864` to `3788361` (AMReX-Combustion#438)

Bumps [Submodules/AMReX-Hydro](https://github.com/AMReX-Fluids/AMReX-Hydro) from `3ab9864` to `3788361`.
- [Commits](AMReX-Fluids/AMReX-Hydro@3ab9864...3788361)

---
updated-dependencies:
- dependency-name: Submodules/AMReX-Hydro
  dependency-type: direct:production
...

Signed-off-by: dependabot[bot] <[email protected]>
Co-authored-by: dependabot[bot] <49699333+dependabot[bot]@users.noreply.github.com>

* Bump Submodules/PelePhysics from `327aadb` to `f54b084` (AMReX-Combustion#439)

Bumps [Submodules/PelePhysics](https://github.com/AMReX-Combustion/PelePhysics) from `327aadb` to `f54b084`.
- [Release notes](https://github.com/AMReX-Combustion/PelePhysics/releases)
- [Commits](AMReX-Combustion/PelePhysics@327aadb...f54b084)

---
updated-dependencies:
- dependency-name: Submodules/PelePhysics
  dependency-type: direct:production
...

Signed-off-by: dependabot[bot] <[email protected]>
Co-authored-by: dependabot[bot] <49699333+dependabot[bot]@users.noreply.github.com>

---------

Signed-off-by: dependabot[bot] <[email protected]>
Co-authored-by: Dave Montgomery <[email protected]>
Co-authored-by: Bruce Perry <[email protected]>
Co-authored-by: dependabot[bot] <49699333+dependabot[bot]@users.noreply.github.com>
Co-authored-by: Marc T. Henry de Frahan <[email protected]>

.github/workflows/ci.yml

@@ -363,7 +363,7 @@ jobs:
           echo "deb https://apt.repos.intel.com/oneapi all main" \
               | sudo tee /etc/apt/sources.list.d/oneAPI.list
           sudo apt-get update
-          sudo apt-get install -y intel-oneapi-compiler-dpcpp-cpp intel-oneapi-mkl-devel
+          sudo apt-get install -y intel-oneapi-compiler-dpcpp-cpp=2024.2.1-1079 intel-oneapi-mkl-devel=2024.2.2-15
       - name: Install Ccache
         run: |
           wget https://github.com/ccache/ccache/releases/download/v4.8/ccache-4.8-linux-x86_64.tar.xz

CMake/BuildPeleExe.cmake

@@ -25,7 +25,6 @@ function(build_pele_exe pele_exe_name pele_physics_lib_name)
        ${SRC_DIR}/PeleLMeX_DeriveUserDefined.cpp
        ${SRC_DIR}/PeleLMeX_DiffusionOp.H
        ${SRC_DIR}/PeleLMeX_DiffusionOp.cpp
-       ${SRC_DIR}/PeleLMeX_EBUserDefined.H
        ${SRC_DIR}/PeleLMeX_FlowControllerData.H
        ${SRC_DIR}/PeleLMeX.H
        ${SRC_DIR}/PeleLMeX.cpp
@@ -46,10 +45,12 @@ function(build_pele_exe pele_exe_name pele_physics_lib_name)
        ${SRC_DIR}/PeleLMeX_Evolve.cpp
        ${SRC_DIR}/PeleLMeX_FlowController.cpp
        ${SRC_DIR}/PeleLMeX_Forces.cpp
+       ${SRC_DIR}/PeleLMeX_Init.cpp
+       ${SRC_DIR}/PeleLMeX_ODEQty.cpp
        ${SRC_DIR}/PeleLMeX_PatchFlowVariables.H
        ${SRC_DIR}/PeleLMeX_PatchFlowVariables.cpp
-       ${SRC_DIR}/PeleLMeX_Init.cpp
        ${SRC_DIR}/PeleLMeX_Plot.cpp
+       ${SRC_DIR}/PeleLMeX_ProblemSpecificFunctions.H
        ${SRC_DIR}/PeleLMeX_Projection.cpp
        ${SRC_DIR}/PeleLMeX_Reactions.cpp
        ${SRC_DIR}/PeleLMeX_Regrid.cpp
@@ -69,6 +70,26 @@ function(build_pele_exe pele_exe_name pele_physics_lib_name)
        ${SRC_DIR}/main.cpp
   )
 
+  if(PELE_EB_USER_DEFINED)
+    target_sources(${pele_exe_name}
+      PRIVATE
+      ${CMAKE_CURRENT_SOURCE_DIR}/PeleLMeX_EBUserDefined.H)
+  else()
+    target_sources(${pele_exe_name}
+      PRIVATE
+        ${SRC_DIR}/PeleLMeX_EBUserDefined.H)
+  endif()
+
+  if(PELELM_USER_DEFINED_EXT_SRC)
+    target_sources(${pele_exe_name}
+      PRIVATE
+      ${CMAKE_CURRENT_SOURCE_DIR}/PeleLMeX_ProblemSpecificFunctions.cpp)
+  else()
+    target_sources(${pele_exe_name}
+      PRIVATE
+        ${SRC_DIR}/PeleLMeX_ProblemSpecificFunctions.cpp)
+  endif()
+
   if(PELE_PHYSICS_ENABLE_SOOT)
     target_sources(${pele_exe_name}
       PRIVATE
@@ -88,6 +109,10 @@ function(build_pele_exe pele_exe_name pele_physics_lib_name)
         ${SRC_DIR}/PeleLMeX_Radiation.cpp)
   endif()
 
+  if (PELELM_NUM_ODE GREATER 0)
+    target_compile_definitions(${pele_exe_name} PRIVATE NUM_ODE=${PELELM_NUM_ODE})
+  endif()
+
   if(NOT "${pele_exe_name}" STREQUAL "${PROJECT_NAME}-UnitTests")
     target_sources(${pele_exe_name}
        PRIVATE

CMake/BuildPelePhysicsLib.cmake

@@ -69,6 +69,13 @@ function(build_pele_physics_lib pele_physics_lib_name)
     )
     target_include_directories(${pele_physics_lib_name} PUBLIC ${PELE_PHYSICS_UTILITY_DIR}/PMF)
 
+    target_sources(${pele_physics_lib_name}
+      PRIVATE
+        ${PELE_PHYSICS_UTILITY_DIR}/Utilities/Utilities.H
+        ${PELE_PHYSICS_UTILITY_DIR}/Utilities/UnitConversions.H
+    )
+    target_include_directories(${pele_physics_lib_name} PUBLIC ${PELE_PHYSICS_UTILITY_DIR}/Utilities)
+
     target_sources(${pele_physics_lib_name}
       PRIVATE
         ${AMREX_SUNDIALS_DIR}/AMReX_Sundials.H

Docs/sphinx/manual/LMeXControls.rst

@@ -269,6 +269,8 @@ PeleLMeX algorithm
     peleLM.spark1.radius = 1e-3            # [OPT] Radius of the spark [m]
     peleLM.spark1.duration = 1e-3          # [OPT] Duration of the spark [s]
     peleLM.spark1.time = 1e-2              # [OPT] Time when spark starts [s]
+
+    peleLM.user_defined_ext_sources = 0    # [OPT, DEF=0] Enable user defined source terms. Requires local ProblemSpecificFunctions.cpp.
     
     
 Transport coefficients and LES

Docs/sphinx/manual/Model.rst

@@ -63,16 +63,16 @@ The set of conservation equations specialized to the low Mach number regime is a
 
     &\frac{\partial (\rho \boldsymbol{u})}{\partial t} +
     \nabla \cdot \left(\rho  \boldsymbol{u} \boldsymbol{u} + \tau \right)
-    = -\nabla \pi + \rho \boldsymbol{F},\\
+    = -\nabla \pi + \boldsymbol{S}_{\text{ext},\rho\boldsymbol{u}},\\
     &\frac{\partial (\rho Y_m)}{\partial t} +
     \nabla \cdot \left( \rho Y_m \boldsymbol{u}
     + \boldsymbol{\mathcal{F}}_{m} \right)
-    = \rho \dot{\omega}_m,\\
+    = \rho \dot{\omega}_m + S_{\text{ext},\rho\boldsymbol{Y_m}},\\
     &\frac{ \partial (\rho h)}{ \partial t} +
     \nabla \cdot \left( \rho h \boldsymbol{u}
-    + \boldsymbol{\mathcal{Q}} \right) = 0 ,
+    + \boldsymbol{\mathcal{Q}} \right) = S_{\text{ext},\rho\boldsymbol{h}} ,
 
-where :math:`\rho` is the density, :math:`\boldsymbol{u}` is the velocity, :math:`h` is the mass-weighted enthalpy, :math:`T` is temperature and :math:`Y_m` is the mass fraction of species :math:`m`. :math:`\dot{\omega}_m` is the molar production rate for species :math:`m`, the modeling of which will be described later in this section. :math:`\tau` is the stress tensor, :math:`\boldsymbol{\mathcal{Q}}` is the heat flux and :math:`\boldsymbol{\mathcal{F}}_m` are the species diffusion fluxes. These transport fluxes require the evaluation of transport coefficients (e.g., the viscosity :math:`\mu`, the conductivity :math:`\lambda` and the diffusivity matrix :math:`D`) which are computed using the library EGLIB, as will be described in more depth in the diffusion section. The momentum source, :math:`\boldsymbol{F}`, is an external forcing term.  For example, we have used :math:`\boldsymbol{F}` to implement a long-wavelength time-dependent force to establish and maintain quasi-stationary turbulence.
+where :math:`\rho` is the density, :math:`\boldsymbol{u}` is the velocity, :math:`h` is the mass-weighted enthalpy, :math:`T` is temperature and :math:`Y_m` is the mass fraction of species :math:`m`. :math:`\dot{\omega}_m` is the molar production rate for species :math:`m`, the modeling of which will be described later in this section. :math:`\tau` is the stress tensor, :math:`\boldsymbol{\mathcal{Q}}` is the heat flux and :math:`\boldsymbol{\mathcal{F}}_m` are the species diffusion fluxes. These transport fluxes require the evaluation of transport coefficients (e.g., the viscosity :math:`\mu`, the conductivity :math:`\lambda` and the diffusivity matrix :math:`D`) which are computed using the library EGLIB, as will be described in more depth in the diffusion section. The source terms, :math:`S_{\text{ext}}`, are user defined external forcing terms.  For example, we have used :math:`\boldsymbol{S}_{\text{ext},\rho\boldsymbol{u}} = \rho\boldsymbol{F}` to implement a long-wavelength time-dependent force to establish and maintain quasi-stationary turbulence.
 
 These evolution equations are supplemented by an equation of state for the thermodynamic pressure.  For example, the ideal gas law,
 
@@ -94,7 +94,7 @@ Neither species diffusion nor reactions redistribute the total mass, hence we ha
 
 .. math::
 
-    \frac{\partial \rho}{\partial t} + \nabla \cdot \rho \boldsymbol{u} = 0
+    \frac{\partial \rho}{\partial t} + \nabla \cdot \rho \boldsymbol{u} = S_{\text{ext},\rho}
 
 This, together with the conservation equations form a differential-algebraic equation (DAE) system that describes an evolution subject to a constraint.  A standard approach to attacking such a system computationally is to differentiate the constraint until it can be recast as an initial value problem.  Following this procedure, we set the thermodynamic pressure constant in the frame of the fluid,
 
@@ -107,16 +107,47 @@ and observe that if the initial conditions satisfy the constraint, an evolution
 .. math::
 
     \nabla \cdot \boldsymbol{u} = \frac{1}{T}\frac{DT}{Dt}
-    + W \sum_m \frac{1}{W_m} \frac{DY_m}{Dt} = S
+    + W \sum_m \frac{1}{W_m} \frac{DY_m}{Dt} + \frac{1}{\rho}S_{\text{ext},\rho} = S.
 
-The constraint here take the form of a condition on the divergence of the flow.  Note that the actual expressions to use here will depend upon the chosen models for evaluating the transport fluxes.
+The constraint here takes the form of a condition on the divergence of the flow.  Note that the actual expressions to use here will depend upon the chosen models for evaluating the transport fluxes.
 
-For the standard ideal gas EOS, the divergence constraint on velocity becomes:
+For the standard ideal gas EOS, 
 
 .. math::
 
-    \nabla \cdot \boldsymbol{u} &= \frac{1}{\rho c_p T} \left(\nabla \cdot \lambda\nabla T - \sum_m \boldsymbol{\Gamma_m} \cdot \nabla h_m \right) \\
-    &- \frac{1}{\rho} \sum_m \frac{W}{W_m} \nabla \cdot \boldsymbol{\Gamma_m} + \frac{1}{\rho}\sum_m \left(\frac{W}{W_m} - \frac{h_m}{c_p T} \right) \dot \omega \equiv S .
+    \frac{DT}{Dt} &= \frac{1}{\rho}\Big[ - \nabla \cdot \boldsymbol{Q} + S_{\text{ext},\rho h} - h S_{\text{ext},\rho}\Big] - \sum_m \frac{h_m}{ c_p} \frac{DY_m}{Dt}, \\
+    \frac{DY_m}{Dt} &= \frac{1}{\rho}\Big[ - \nabla \cdot \boldsymbol{\mathcal{F_m}} + \rho \dot \omega_m + S_{\text{ext},\rho Y_m} - Y_m S_{\text{ext},\rho}\Big].
+
+Therefore, the divergence constraint on velocity becomes:
+
+.. math::
+
+    \nabla \cdot \boldsymbol{u} &= \frac{1}{\rho c_p T} \Big(-\nabla \cdot \boldsymbol{Q} + S_{\text{ext},\rho h} - h S_{\text{ext},\rho}\Big) \\
+    &\;\;\;\; +  \sum_m \bigg( \frac{W}{\rho W_m} -  \frac{h_m}{\rho c_p T}\bigg)\bigg( - \nabla \cdot \boldsymbol{\mathcal{F}}_m + \rho \dot \omega_m + S_{\text{ext},\rho Y_m} - Y_m S_{\text{ext},\rho}\bigg) + \frac{1}{\rho} S_{\text{ext},\rho}\equiv S .
+
+However, it can be shown that 
+
+.. math::
+
+    \sum_m \frac{W}{\rho W_m} Y_m S_{\text{ext},\rho} = \frac{1}{\rho}S_{\text{ext},\rho}
+
+and 
+
+.. math::
+    \sum_m h_m Y_m S_{\text{ext},\rho} = h S_{\text{ext},\rho}.
+
+Thus, the terms containing :math:`S_{\text{ext},\rho}` cancel and the divergence constraint for the standard ideal gas EOS simplifies to: 
+
+.. math::
+
+    \nabla \cdot \boldsymbol{u} = \frac{1}{\rho c_p T} \Big(-\nabla \cdot \boldsymbol{Q} + S_{\text{ext},\rho h} \Big) 
+    +  \sum_m \bigg( \frac{W}{\rho W_m} -  \frac{h_m}{\rho c_p T}\bigg)\bigg( - \nabla \cdot \boldsymbol{\mathcal{F}}_m + \rho \dot \omega_m + S_{\text{ext},\rho Y_m} \bigg) \equiv S .
+
+In addition to the flow equations, `PeleLMeX` can also solve for a set of quantities that are neither advected nor diffused, satisfying:
+
+.. math::
+
+    \frac{\partial B_k}{\partial t} = S_{\text{ext},B_k}.
 
 Confined domain ambient pressure
 ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^

Exec/Make.PeleLMeX

@@ -57,6 +57,10 @@ ifeq ($(USE_RADIATION), TRUE)
    DEFINES += -DPELE_USE_RADIATION
 endif
 
+ifeq ($(shell test 0$(PELELM_NUM_ODE) -gt 0; echo $$?), 0)
+   DEFINES += -DNUM_ODE=$(PELELM_NUM_ODE)
+endif
+
 Bpack += $(foreach dir, $(LMdirs), $(PELE_HOME)/$(dir)/Make.package)
 Blocs += $(foreach dir, $(LMdirs), $(PELE_HOME)/$(dir))
 
@@ -73,6 +77,19 @@ endif
 ifeq ($(Eos_Model),$(filter $(Eos_Model),Soave-Redlich-Kwong))
    DEFINES += -DUSE_SRK_EOS
 endif
+ifeq ($(Eos_Model),$(filter $(Eos_Model),Manifold))
+   DEFINES += -DUSE_MANIFOLD_EOS
+   DEFINES += -DMANIFOLD_DIM=$(Manifold_Dim)
+   ifeq ($(Manifold_Type),Table)
+      DEFINES += -DMANIFOLD_EOS_TYPE=1
+   else
+      ifeq ($(Manifold_Type),Network)
+         DEFINES += -DMANIFOLD_EOS_TYPE=2
+      else
+         DEFINES += -DMANIFOLD_EOS_TYPE=0
+      endif
+   endif
+endif
 
 # Transport model switches
 ifeq ($(Transport_Model), Simple)
@@ -88,6 +105,9 @@ endif
 ifeq ($(Transport_Model), Sutherland)
   DEFINES += -DUSE_SUTHERLAND_TRANSPORT
 endif
+ifeq ($(Transport_Model), Manifold)
+   DEFINES += -DUSE_MANIFOLD_TRANSPORT
+endif
 
 ifeq ($(PELE_USE_KLU), TRUE)
   DEFINES += -DPELE_USE_KLU
@@ -107,6 +127,7 @@ ChemDir  = Mechanisms/$(Chemistry_Model)
 
 PPdirs := Source/Utility/PMF Source/Utility/TurbInflow Source/Utility/PltFileManager
 PPdirs += Source/Utility/Diagnostics Source/Utility/BlackBoxFunction
+PPdirs += Source/Utility/Utilities
 PPdirs += Source $(ChemDir) Source/Reactions Source/Eos Source/Transport
 Bpack += $(foreach dir, $(PPdirs), $(PELE_PHYSICS_HOME)/$(dir)/Make.package)
 Blocs += $(foreach dir, $(PPdirs), $(PELE_PHYSICS_HOME)/$(dir))

Exec/RegTests/CMakeLists.txt

@@ -3,6 +3,7 @@ if(PELE_ENABLE_EB)
   add_subdirectory(EB_EnclosedFlame)
   add_subdirectory(EB_EnclosedVortex)
   add_subdirectory(EB_FlowPastCylinder)
+  add_subdirectory(EB_ODEQty)
   if(PELE_DIM EQUAL 3)
     add_subdirectory(EB_PipeFlow)
   endif()

Exec/RegTests/EB_BackwardStepFlame/CMakeLists.txt

@@ -5,4 +5,5 @@ set(PELE_PHYSICS_ENABLE_SPRAY OFF)
 set(PELE_PHYSICS_SPRAY_FUEL_NUM 0)
 set(PELE_PHYSICS_ENABLE_SOOT OFF)
 set(PELE_PHYSICS_ENABLE_RADIATION OFF)
+set(PELE_EB_USER_DEFINED ON)
 include(BuildExeAndLib)

Exec/RegTests/EB_EnclosedFlame/CMakeLists.txt

@@ -5,4 +5,5 @@ set(PELE_PHYSICS_ENABLE_SPRAY OFF)
 set(PELE_PHYSICS_SPRAY_FUEL_NUM 0)
 set(PELE_PHYSICS_ENABLE_SOOT OFF)
 set(PELE_PHYSICS_ENABLE_RADIATION OFF)
+set(PELE_EB_USER_DEFINED OFF)
 include(BuildExeAndLib)

Exec/RegTests/EB_EnclosedVortex/CMakeLists.txt

@@ -5,4 +5,5 @@ set(PELE_PHYSICS_ENABLE_SPRAY OFF)
 set(PELE_PHYSICS_SPRAY_FUEL_NUM 0)
 set(PELE_PHYSICS_ENABLE_SOOT OFF)
 set(PELE_PHYSICS_ENABLE_RADIATION OFF)
+set(PELE_EB_USER_DEFINED OFF)
 include(BuildExeAndLib)

Exec/RegTests/EB_FlowPastCylinder/CMakeLists.txt

@@ -5,4 +5,5 @@ set(PELE_PHYSICS_ENABLE_SPRAY OFF)
 set(PELE_PHYSICS_SPRAY_FUEL_NUM 0)
 set(PELE_PHYSICS_ENABLE_SOOT OFF)
 set(PELE_PHYSICS_ENABLE_RADIATION OFF)
+set(PELE_EB_USER_DEFINED ON)
 include(BuildExeAndLib)

Exec/RegTests/EB_ODEQty/CMakeLists.txt

@@ -0,0 +1,11 @@
+set(PELE_PHYSICS_EOS_MODEL Fuego)
+set(PELE_PHYSICS_CHEMISTRY_MODEL air)
+set(PELE_PHYSICS_TRANSPORT_MODEL Constant)
+set(PELE_PHYSICS_ENABLE_SPRAY OFF)
+set(PELE_PHYSICS_SPRAY_FUEL_NUM 0)
+set(PELE_PHYSICS_ENABLE_SOOT OFF)
+set(PELE_PHYSICS_ENABLE_RADIATION OFF)
+set(PELE_EB_USER_DEFINED OFF)
+set(PELELM_USER_DEFINED_EXT_SRC ON)
+set(PELELM_NUM_ODE 3)
+include(BuildExeAndLib)

Exec/RegTests/EB_ODEQty/GNUmakefile

@@ -0,0 +1,38 @@
+# AMReX
+DIM = 2
+COMP = gnu
+PRECISION = DOUBLE
+USE_EB = TRUE
+USE_HYPRE = FALSE
+
+# Profiling
+PROFILE = FALSE
+TINY_PROFILE = FALSE
+COMM_PROFILE = FALSE
+TRACE_PROFILE = FALSE
+MEM_PROFILE = FALSE
+USE_GPROF = FALSE
+
+# Performance
+USE_MPI = TRUE
+USE_OMP = FALSE
+USE_CUDA = FALSE
+USE_HIP = FALSE
+USE_SYCL = FALSE
+
+# Debugging
+DEBUG = FALSE
+FSANITIZER = FALSE
+THREAD_SANITIZER = FALSE
+
+# PeleLMeX
+CEXE_sources += PeleLMeX_ProblemSpecificFunctions.cpp
+PELELM_NUM_ODE = 3
+
+# PelePhysics
+Chemistry_Model = air
+Eos_Model = Fuego
+Transport_Model = Constant
+
+PELE_HOME ?= ../../..
+include $(PELE_HOME)/Exec/Make.PeleLMeX

Exec/RegTests/EB_ODEQty/PeleLMeX_ProblemSpecificFunctions.cpp

@@ -0,0 +1,55 @@
+#include <PeleLMeX.H>
+#include <PeleLMeX_K.H>
+#include <PeleLMeX_ProblemSpecificFunctions.H>
+
+using namespace amrex;
+
+/*
+Problem specific functions:
+- This file must be copied locally to the case directory
+- Add the following to GNUmakefile:
+          CEXE_sources += PeleLMeX_ProblemSpecificFunctions.cpp
+- Modify as needed
+*/
+
+#if NUM_ODE > 0
+void
+set_ode_names(Vector<std::string>& a_ode_names)
+{
+  a_ode_names.resize(NUM_ODE);
+  for (int n = 0; n < NUM_ODE; n++) {
+    a_ode_names[n] = "MY_ODE_" + std::to_string(n);
+  }
+}
+#endif
+
+void
+problem_modify_ext_sources(
+  Real /*time*/,
+  Real /*dt*/,
+  const MultiFab& state_old,
+  const MultiFab& /*state_new*/,
+  std::unique_ptr<MultiFab>& ext_src,
+  const GeometryData& /*geomdata*/,
+  const ProbParm& prob_parm)
+{
+  /*
+  Notes:
+    1) ext_src contains sources from velocity forcing coming in.
+       This function should add to rather than overwrite ext_src.
+    2) Requires "peleLM.user_defined_ext_sources = true" in input file
+  */
+
+  auto ext_src_arr = ext_src->arrays();
+  auto const& state_old_arr = state_old.const_arrays();
+
+  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;
+      }
+    });
+  Gpu::streamSynchronize();
+}

Exec/RegTests/EB_ODEQty/README.md

@@ -0,0 +1,2 @@
+## 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$.