add the case to PeleLMeX

Exec/RegTests/SootRadTest/GNUmakefile

@@ -0,0 +1,33 @@
+# AMReX
+DIM             = 2
+DEBUG           = FALSE
+PRECISION       = DOUBLE
+VERBOSE         = FALSE
+TINY_PROFILE    = TRUE
+
+# Compilation
+COMP            = llvm
+USE_MPI         = TRUE
+USE_OMP         = FALSE
+USE_CUDA        = FALSE
+USE_HIP         = FALSE
+
+# PeleLMeX
+USE_EFIELD      = FALSE
+
+# PelePhysics
+Chemistry_Model = SootReaction
+Eos_Model	= Fuego
+Transport_Model = Simple
+
+USE_SOOT = TRUE
+NUM_SOOT_MOMENTS = 3
+
+USE_PELERAD = TRUE
+ifeq ($(USE_PELERAD), TRUE)
+ifeq ($(USE_HIP), TRUE)
+LIBRARIES += -lstdc++fs
+endif
+endif
+
+include $(PELELMEX_HOME)/Exec/Make.PeleLMeX

Exec/RegTests/SootRadTest/README.md

@@ -0,0 +1,4 @@
+## SootRadTest
+Testing of the coupling between PeleLMeX and PeleMP soot and radiation modules.
+Please specify the radiation database path in the input file before running.
+echo "pelerad.kppath = "$PELERAD_HOME/data/kpDB/"" >> first-input

Exec/RegTests/SootRadTest/first-input

@@ -0,0 +1,112 @@
+#----------------------DOMAIN DEFINITION------------------------
+geometry.is_periodic = 0 1 1               # For each dir, 0: non-perio, 1: periodic
+geometry.coord_sys   = 0                  # 0 => cart, 1 => RZ
+geometry.prob_lo     = 0.0 0.0 0.0        # x_lo y_lo (z_lo)
+geometry.prob_hi     = 0.04 0.0025 0.0025        # x_hi y_hi (z_hi)
+
+# >>>>>>>>>>>>>  BC FLAGS <<<<<<<<<<<<<<<<
+# Interior, Inflow, Outflow, Symmetry,
+# SlipWallAdiab, NoSlipWallAdiab, SlipWallIsotherm, NoSlipWallIsotherm
+peleLM.lo_bc = Inflow Interior Interior
+peleLM.hi_bc = Outflow Interior Interior
+
+#-------------------------AMR CONTROL----------------------------
+amr.n_cell          = 128 8 8      # Level 0 number of cells in each direction   
+amr.v               = 1                # AMR verbose
+amr.max_level       = 2                # maximum level number allowed
+amr.regrid_int      = 4                # how often to regrid
+amr.n_error_buf     = 1 1 2 2          # number of buffer cells in error est
+amr.grid_eff        = 0.7              # what constitutes an efficient grid
+amr.ref_ratio = 2 2 2
+amr.blocking_factor = 8
+amr.max_grid_size = 128
+
+#--------------------------- Problem -------------------------------
+prob.P_mean = 98700.
+prob.standoff = 0.0
+pmf.datafile = "datafile_init/mueller_burner.dat"
+pmf.do_cellAverage = 0
+
+#--------------------SOOT MODELING------------------------
+peleLM.do_soot_solve = 1
+soot.incept_pah = A2 # Soot inception species
+soot.v = 0
+soot.temp_cutoff = 290.
+soot.conserve_mass = false
+soot.num_subcycles = 10
+soot.max_subcycles = 1000
+
+#-------------------------PeleLM CONTROL----------------------------
+peleLM.v = 1
+peleLM.incompressible = 0
+peleLM.use_wbar = 0
+peleLM.sdc_iterMax = 1
+peleLM.floor_species = 0
+peleLM.advection_scheme = Godunov_BDS
+
+peleLM.do_temporals = 0
+peleLM.temporal_int = 2
+peleLM.mass_balance = 1
+peleLM.num_init_iter = 1
+peleLM.plot_react = 0
+
+#amr.restart = chk00005
+#amr.check_int = 2000
+amr.plot_per = 1.E-3
+amr.dt_shrink = 0.1
+amr.max_step = 10000
+amr.stop_time = 0.022
+amr.cfl = 0.3
+amr.derive_plot_vars = rhoRT mass_fractions rhominsumrhoY
+#amr.fixed_dt = 0.008
+#amr.fixed_dt = 1.E-6
+
+# --------------- INPUTS TO CHEMISTRY REACTOR ---------------
+peleLM.chem_integrator = "ReactorCvode"
+peleLM.use_typ_vals_chem = 0          # Use species/temp typical values in CVODE
+# ode.rtol = 1.0e-6                     # Relative tolerance of the chemical solve
+# ode.atol = 1.0e-5                     # Absolute tolerance factor applied on typical values
+#cvode.solve_type = GMRES
+#cvode.solve_type = denseAJ_direct     # CVODE Linear solve type (for Newton direction) 
+cvode.solve_type = magma_direct     # CVODE Linear solve type (for Newton direction) 
+#cvode.max_order  = 4                  # CVODE max BDF order.
+#ode.atol = 1.E-12
+
+mac_proj.verbose = 0
+mac_proj.atol = 1.E-14
+mac_proj.rtol = 1.E-11
+nodal_proj.verbose = 0
+nodal_proj.atol        = 6.0e-14     # tolerence for projections
+nodal_proj.rtol        = 6.0e-11
+
+#--------------------REFINEMENT CONTROL------------------------
+amr.refinement_indicators = gradT
+amr.gradT.max_level     = 2
+amr.gradT.adjacent_difference_greater = 30.
+amr.gradT.field_name    = temp
+
+amrex.regtest_reduction = 1
+amrex.fpe_trap_invalid = 1
+amrex.fpe_trap_zero = 1
+amrex.fpe_trap_overflow = 1
+
+#--------------------RADIATION MODELING------------------------
+peleLM.do_rad_solve = 1
+pelerad.composite_solve = 1
+pelerad.use_hypre = 0
+pelerad.verbose = 0
+pelerad.max_iter = 200
+pelerad.max_coarsening_level = 10
+pelerad.reltol = 1.0e-3
+pelerad.abstol = 1.0e-3
+pelerad.bottom_reltol = 1.0e-6
+pelerad.bottom_abstol = 1.0e-6
+pelerad.agglomeration = 1
+pelerad.consolidation = 0
+pelerad.bottom_verbose = 0
+pelerad.maxorder = 2
+pelerad.linop_maxorder = 2
+pelerad.max_fmg_iter = 0
+pelerad.lo_bc = Robin Periodic Periodic
+pelerad.hi_bc = Robin Periodic Periodic
+#please set the pelerad.kppath

Exec/RegTests/SootRadTest/pelelmex_prob.H

@@ -0,0 +1,140 @@
+#ifndef PELELM_PROB_H
+#define PELELM_PROB_H
+
+#include <AMReX_Geometry.H>
+#include <AMReX_FArrayBox.H>
+#include <AMReX_ParmParse.H>
+
+#include <pelelmex_prob_parm.H>
+#include <PMF.H>
+#include <PMFData.H>
+
+#include <PeleLMeX_Index.H>
+#include <PelePhysics.H>
+#include "SootModel.H"
+
+AMREX_GPU_DEVICE
+AMREX_FORCE_INLINE
+void
+pelelmex_initdata(
+  int i,
+  int j,
+  int k,
+  int /*is_incompressible*/,
+  amrex::Array4<amrex::Real> const& state,
+  amrex::Array4<amrex::Real> const& aux,
+  amrex::GeometryData const& geomdata,
+  ProbParm const& prob_parm,
+  pele::physics::PMF::PmfData::DataContainer const* pmf_data)
+{
+
+  const amrex::Real* prob_lo = geomdata.ProbLo();
+  const amrex::Real* prob_hi = geomdata.ProbHi();
+  const amrex::Real* dx = geomdata.CellSize();
+
+  AMREX_D_TERM(const amrex::Real x = prob_lo[0] + (i + 0.5) * dx[0];
+               , const amrex::Real y = prob_lo[1] + (j + 0.5) * dx[1];
+               , const amrex::Real z = prob_lo[2] + (k + 0.5) * dx[2];);
+  amrex::GpuArray<amrex::Real, NUM_SPECIES> massfrac = {{0.0}};
+  amrex::GpuArray<amrex::Real, NUM_SPECIES + 4> pmf_vals = {{0.0}};
+  amrex::Real x1 = (x - prob_parm.standoff - 0.5 * dx[0]) * 100.;
+  amrex::Real x2 = (x - prob_parm.standoff + 0.5 * dx[0]) * 100.;
+  pele::physics::PMF::pmf(pmf_data, x1, x2, pmf_vals);
+  state(i, j, k, TEMP) = pmf_vals[1];
+  amrex::Real norm = 0.;
+  for (int n = 0; n < NUM_SPECIES; n++) {
+    massfrac[n] = amrex::max(0., amrex::min(1., pmf_vals[3 + n]));
+    norm += massfrac[n];
+  }
+  for (int n = 0; n < NUM_SPECIES; ++n) {
+    massfrac[n] = massfrac[n] / norm;
+  }
+  AMREX_D_TERM(state(i, j, k, VELX) = pmf_vals[0] * 1.E-2;
+               , state(i, j, k, VELY) = 0.;, state(i, j, k, VELZ) = 0.;);
+  amrex::Real rho_cgs, P_cgs;
+  P_cgs = prob_parm.P_mean * 10.;
+
+  auto eos = pele::physics::PhysicsType::eos();
+  eos.PYT2R(P_cgs, massfrac.data(), state(i, j, k, TEMP), rho_cgs);
+  state(i, j, k, DENSITY) = rho_cgs * 1.0e3; // CGS -> MKS conversion
+
+  eos.TY2H(state(i, j, k, TEMP), massfrac.data(), state(i, j, k, RHOH));
+  state(i, j, k, RHOH) = state(i, j, k, RHOH) * 1.0e-4 *
+                         state(i, j, k, DENSITY); // CGS -> MKS conversion
+
+  for (int n = 0; n < NUM_SPECIES; n++) {
+    state(i, j, k, FIRSTSPEC + n) = massfrac[n] * state(i, j, k, DENSITY);
+  }
+  for (int is = 0; is < NUM_SOOT_MOMENTS + 1; ++is) {
+    state(i, j, k, FIRSTSOOT + is) = prob_parm.soot_vals[is];
+  }
+}
+
+AMREX_GPU_DEVICE
+AMREX_FORCE_INLINE
+void
+bcnormal(
+  const amrex::Real x[AMREX_SPACEDIM],
+  const int /*m_nAux*/,
+  amrex::Real s_ext[NVAR],
+  const int idir,
+  const int sgn,
+  const amrex::Real time,
+  amrex::GeometryData const& geomdata,
+  ProbParm const& prob_parm,
+  pele::physics::PMF::PmfData::DataContainer const* pmf_data)
+{
+  const amrex::Real* prob_lo = geomdata.ProbLo();
+  amrex::GpuArray<amrex::Real, NUM_SPECIES + 4> pmf_vals = {{0.0}};
+  amrex::Real massfrac[NUM_SPECIES] = {0.0};
+  if (sgn == 1) {
+    pele::physics::PMF::pmf(pmf_data, prob_lo[idir], prob_lo[idir], pmf_vals);
+    AMREX_D_TERM(s_ext[VELX] = pmf_vals[0] * 1.E-2;, s_ext[VELY] = 0.0;
+                 , s_ext[VELZ] = 0.0;);
+
+    s_ext[TEMP] = pmf_vals[1];
+
+    for (int n = 0; n < NUM_SPECIES; n++) {
+      massfrac[n] = pmf_vals[3 + n];
+    }
+
+    amrex::Real rho_cgs, P_cgs, RhoH_temp;
+    P_cgs = prob_parm.P_mean * 10.0;
+
+    auto eos = pele::physics::PhysicsType::eos();
+    eos.PYT2R(P_cgs, massfrac, s_ext[TEMP], rho_cgs);
+    s_ext[DENSITY] = rho_cgs * 1.0e3;
+
+    eos.TY2H(s_ext[TEMP], massfrac, RhoH_temp);
+    s_ext[RHOH] = RhoH_temp * 1.0e-4 * s_ext[DENSITY]; // CGS -> MKS conversion
+
+    for (int n = 0; n < NUM_SPECIES; n++) {
+      s_ext[FIRSTSPEC + n] = massfrac[n] * s_ext[DENSITY];
+    }
+    for (int is = 0; is < NUM_SOOT_MOMENTS + 1; ++is) {
+      s_ext[FIRSTSOOT + is] = prob_parm.soot_vals[is];
+    }
+  }
+}
+
+AMREX_GPU_DEVICE
+AMREX_FORCE_INLINE
+void
+zero_visc(
+  int i,
+  int j,
+  int k,
+  amrex::Array4<amrex::Real> const& beta,
+  amrex::GeometryData const& geomdata,
+  amrex::Box const& domainBox,
+  const int dir,
+  const int beta_comp,
+  const int nComp)
+{
+  amrex::ignore_unused(
+    i, j, k, beta, geomdata, domainBox, dir, beta_comp, nComp);
+  // We treat species when beta_comp == 0 and nComp == NUM_SPECIES
+  // otherwise this routine could be called for other face diffusivity (Temp,
+  // velocity, ...)
+}
+#endif

Exec/RegTests/SootRadTest/pelelmex_prob.cpp

@@ -0,0 +1,20 @@
+#include <PeleLMeX.H>
+#include <pelelmex_prob.H>
+
+void
+PeleLM::readProbParm()
+{
+  // Parse params
+  amrex::ParmParse pp("prob");
+  pp.query("P_mean", PeleLM::prob_parm->P_mean);
+  pp.query("standoff", PeleLM::prob_parm->standoff);
+  PeleLM::pmf_data.initialize();
+  amrex::Real moments[NUM_SOOT_MOMENTS + 1] = {0.0};
+  if (PeleLM::do_soot_solve) {
+    SootData* const sd = PeleLM::soot_model->getSootData();
+    sd->initialSmallMomVals(moments);
+  }
+  for (int n = 0; n < NUM_SOOT_MOMENTS + 1; ++n) {
+    PeleLM::prob_parm->soot_vals[n] = moments[n];
+  }
+}

Exec/RegTests/SootRadTest/pelelmex_prob_parm.H

@@ -0,0 +1,16 @@
+#ifndef PELELM_PROB_PARM_H
+#define PELELM_PROB_PARM_H
+
+#include <AMReX_REAL.H>
+#include <AMReX_GpuMemory.H>
+
+using namespace amrex::literals;
+
+struct ProbParm
+{
+  amrex::GpuArray<amrex::Real, NUM_SOOT_MOMENTS + 1> soot_vals;
+  amrex::Real P_mean = 101325.0;
+  amrex::Real standoff = 0.;
+};
+
+#endif

Source/PeleLMeX_Rad.cpp

@@ -73,7 +73,6 @@ PeleLM::computeRadSource(
 
   BL_PROFILE_VAR("PeleLM::advance::rad::solve", PLM_RAD_SOLV);
   rad_model->evaluateRad();
-  BL_PROFILE_VAR_STOP(PLM_RAD_SOLV);
 
   for (int lev = 0; lev <= finest_level; lev++) {
 #ifdef AMREX_USE_OMP
@@ -85,4 +84,5 @@ PeleLM::computeRadSource(
       rad_model->calcRadSource(mfi, source_arr, lev);
     }
   }
+  BL_PROFILE_VAR_STOP(PLM_RAD_SOLV);
 }

Source/PeleLMeX_Setup.cpp

@@ -15,10 +15,6 @@
 #include "SootModel.H"
 #endif
 
-#ifdef PELELM_USE_RAD
-#include <PeleLMRad.hpp>
-#endif
-
 using namespace amrex;
 
 static Box