Adding the radiation module

Exec/Make.PeleLMeX

@@ -54,6 +54,15 @@ ifeq ($(USE_SOOT), TRUE)
    DEFINES += -DNUM_SOOT_MOMENTS=$(NUM_SOOT_MOMENTS)
 endif
 
+ifeq ($(USE_PELERAD), TRUE)
+   $(info ******** Radiation model is activated ********)
+   DEFINES += -DPELELM_USE_RAD
+   ifeq ($(USE_HIP), TRUE)
+       $(info ******** Radiation module is using hip ********)
+       DEFINES += -DPELERAD_USE_HIP
+   endif
+endif
+
 Bpack += $(foreach dir, $(LMdirs), $(PELELMEX_HOME)/$(dir)/Make.package)
 Blocs += $(foreach dir, $(LMdirs), $(PELELMEX_HOME)/$(dir))
 
@@ -142,6 +151,9 @@ ifeq ($(USE_SOOT), TRUE)
    Bpack += $(PELEMP_HOME)/Source/Soot_Models/Make.package
    Blocs += $(PELEMP_HOME)/Source/Soot_Models
 endif
+ifeq ($(USE_PELERAD), TRUE)
+   Blocs += $(PELERAD_HOME)/src
+endif
 
 #---------------
 # Includes

Source/Make.package

@@ -46,6 +46,12 @@ ifeq ($(USE_SOOT), TRUE)
   CEXE_sources += PeleLMeX_Soot.cpp
 endif
 
+ifeq ($(USE_PELERAD), TRUE)
+  CEXE_sources += PeleLMeX_Rad.cpp
+endif
+
 ifeq ($(USE_PARTICLES), TRUE)
   CEXE_sources += PeleLMeX_SprayParticles.cpp
 endif
+
+

Source/PeleLMeX.H

@@ -32,6 +32,10 @@
 #include <AMReX_ErrorList.H>
 #include <AMReX_VisMF.H>
 
+#ifdef PELELM_USE_RAD
+#include <PeleLMRad.hpp>
+#endif
+
 // Forward declarations
 #ifdef PELELM_USE_SPRAY
 class SprayParticleContainer;
@@ -857,6 +861,15 @@ public:
   //-----------------------------------------------------------------------------
 #endif
 
+#ifdef PELELM_USE_RAD
+  // Initialize the radiation module
+  void RadInit();
+
+  // Compute the radiation source term
+  void computeRadSource(
+    const PeleLM::TimeStamp& a_timestamp, const amrex::Real a_dt);
+#endif
+
   //-----------------------------------------------------------------------------
   // EOS
 
@@ -1679,6 +1692,12 @@ public:
   bool do_soot_solve;
 #endif
 
+#ifdef PELELM_USE_RAD
+  // Radiation parameters
+  std::unique_ptr<PeleRad::Radiation> rad_model;
+  bool do_rad_solve;
+#endif
+
   // Active controller
   int m_ctrl_active{0};
   int m_ctrl_useTemp{0};

Source/PeleLMeX_Advance.cpp

@@ -120,6 +120,13 @@ PeleLM::Advance(int is_initIter)
     computeSootSource(AmrOldTime, m_dt);
   }
 #endif
+#ifdef PELELM_USE_RAD
+  if (do_rad_solve) {
+    BL_PROFILE_VAR("PeleLM::advance::rad", PLM_RAD);
+    computeRadSource(AmrOldTime, m_dt);
+    BL_PROFILE_VAR_STOP(PLM_RAD);
+  }
+#endif
 
   if (m_incompressible == 0) {
     floorSpecies(AmrOldTime);

Source/PeleLMeX_Init.cpp

@@ -187,6 +187,11 @@ PeleLM::initData()
 #ifdef PELELM_USE_SPRAY
     SprayInit();
 #endif
+#ifdef PELELM_USE_RAD
+    if (do_rad_solve) {
+      RadInit();
+    }
+#endif
 
     //----------------------------------------------------------------
     // Set typical values
@@ -249,6 +254,11 @@ PeleLM::initData()
 #ifdef PELELM_USE_SPRAY
     SprayInit();
 #endif
+#ifdef PELELM_USE_RAD
+    if (do_rad_solve) {
+      RadInit();
+    }
+#endif
 #ifdef PELE_USE_EFIELD
     // If restarting from a non efield simulation
     if (m_restart_nonEF) {

Source/PeleLMeX_Plot.cpp

@@ -17,6 +17,10 @@
 #ifdef PELELM_USE_SOOT
 #include "SootModel.H"
 #endif
+#ifdef PELELM_USE_RAD
+#include "PeleLMRad.hpp"
+#endif
+
 using namespace amrex;
 
 namespace {
@@ -116,6 +120,10 @@ PeleLM::WritePlotFile()
   ncomp += 1;
 #endif
 
+#ifdef PELELM_USE_RAD
+  ncomp += 3;
+#endif
+
   // Derive
   int deriveEntryCount = 0;
   for (int ivar = 0; ivar < m_derivePlotVarCount; ivar++) {
@@ -177,6 +185,13 @@ PeleLM::WritePlotFile()
       std::string sootname = soot_model->sootVariableName(mom);
       plt_VarsName.push_back(sootname);
     }
+#endif
+#ifdef PELELM_USE_RAD
+    if (do_rad_solve) {
+      plt_VarsName.push_back("rad.G");
+      plt_VarsName.push_back("rad.kappa");
+      plt_VarsName.push_back("rad.emis");
+    }
 #endif
     if (m_has_divu != 0) {
       plt_VarsName.push_back("divu");
@@ -283,6 +298,16 @@ PeleLM::WritePlotFile()
         mf_plt[lev], m_leveldata_new[lev]->state, FIRSTSOOT, cnt, NUMSOOTVAR,
         0);
       cnt += NUMSOOTVAR;
+#endif
+#ifdef PELELM_USE_RAD
+      if (do_rad_solve) {
+        MultiFab::Copy(mf_plt[lev], rad_model->G()[lev], 0, cnt, 1, 0);
+        cnt += 1;
+        MultiFab::Copy(mf_plt[lev], rad_model->kappa()[lev], 0, cnt, 1, 0);
+        cnt += 1;
+        MultiFab::Copy(mf_plt[lev], rad_model->emis()[lev], 0, cnt, 1, 0);
+        cnt += 1;
+      }
 #endif
       if (m_has_divu != 0) {
         MultiFab::Copy(mf_plt[lev], m_leveldata_new[lev]->divu, 0, cnt, 1, 0);

Source/PeleLMeX_Rad.cpp

@@ -0,0 +1,88 @@
+#include <PeleLMeX.H>
+#include <PeleLMeX_Derive.H>
+
+#include <PeleLMRad.hpp>
+
+void
+PeleLM::RadInit()
+{
+  amrex::Print() << "Radiation model is activated "
+                 << "\n";
+  PeleRad::RadComps rc;
+  amrex::Vector<std::string> names;
+  pele::physics::eos::speciesNames<pele::physics::PhysicsType::eos_type>(names);
+  for (int i = 0; i < NUM_SPECIES; ++i) {
+    if (names[i] == "CO2") {
+      rc.co2Indx = i;
+      continue;
+    }
+    if (names[i] == "H2O") {
+      rc.h2oIndx = i;
+      continue;
+    }
+    if (names[i] == "CO")
+      rc.coIndx = i;
+  }
+  amrex::ParmParse mlmgpp("pelerad");
+
+#ifdef AMREX_USE_EB
+  amrex::Print() << "The EB support is activated in the radiation module \n";
+  rad_model = std::make_unique<PeleRad::Radiation>(
+    geom, grids, dmap, rc, mlmgpp, m_factory);
+#else
+  rad_model =
+    std::make_unique<PeleRad::Radiation>(geom, grids, dmap, rc, mlmgpp);
+#endif
+}
+
+void
+PeleLM::computeRadSource(
+  const PeleLM::TimeStamp& a_timestamp, const amrex::Real a_dt)
+{
+  int const co2Indx = rad_model->readRadIndices().co2Indx;
+  int const h2oIndx = rad_model->readRadIndices().h2oIndx;
+  int const coIndx = rad_model->readRadIndices().coIndx;
+
+  BL_PROFILE_VAR("PeleLM::advance::rad::spec", PLM_RAD_SPEC);
+  for (int lev = 0; lev <= finest_level; lev++) {
+    auto ldata_p = PeleLM::getLevelDataPtr(lev, a_timestamp);
+#ifdef AMREX_USE_OMP
+#pragma omp parallel if (Gpu::notInLaunchRegion())
+#endif
+    for (amrex::MFIter mfi(*(m_extSource[lev]), amrex::TilingIfNotGPU());
+         mfi.isValid(); ++mfi) {
+      auto const& q_yin_co2 =
+        ldata_p->state.const_array(mfi, FIRSTSPEC + co2Indx);
+      auto const& q_yin_h2o =
+        ldata_p->state.const_array(mfi, FIRSTSPEC + h2oIndx);
+      auto const& q_yin_co =
+        ldata_p->state.const_array(mfi, FIRSTSPEC + coIndx);
+      auto const& q_Tin = ldata_p->state.const_array(mfi, TEMP);
+      auto const& q_Pin = ldata_p->state.const_array(mfi, RHORT);
+#ifdef PELELM_USE_SOOT
+      auto const& q_fvin = ldata_p->state.const_array(mfi, FIRSTSOOT + 1);
+      rad_model->updateSpecProp(
+        mfi, q_yin_co2, q_yin_h2o, q_yin_co, q_Tin, q_Pin, q_fvin, lev);
+#else
+      rad_model->updateSpecProp(
+        mfi, q_yin_co2, q_yin_h2o, q_yin_co, q_Tin, q_Pin, lev);
+#endif
+    }
+  }
+  BL_PROFILE_VAR_STOP(PLM_RAD_SPEC);
+
+  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
+#pragma omp parallel if (Gpu::notInLaunchRegion())
+#endif
+    for (amrex::MFIter mfi(*(m_extSource[lev]), amrex::TilingIfNotGPU());
+         mfi.isValid(); ++mfi) {
+      auto const& source_arr = m_extSource[lev]->array(mfi, RHOH);
+      rad_model->calcRadSource(mfi, source_arr, lev);
+    }
+  }
+}

Source/PeleLMeX_Setup.cpp

@@ -14,6 +14,11 @@
 #ifdef PELELM_USE_SOOT
 #include "SootModel.H"
 #endif
+
+#ifdef PELELM_USE_RAD
+#include <PeleLMRad.hpp>
+#endif
+
 using namespace amrex;
 
 static Box
@@ -588,6 +593,13 @@ PeleLM::readParameters()
   }
   soot_model->readSootParams();
 #endif
+#ifdef PELELM_USE_RAD
+  do_rad_solve = false;
+  pp.query("do_rad_solve", do_rad_solve);
+  if (m_verbose && do_rad_solve) {
+    Print() << "Simulation performed with radiation modeling \n";
+  }
+#endif
 }
 
 void