Preliminary implementation of bte coupling for low-mach

docker/test/Dockerfile

@@ -29,7 +29,6 @@ RUN yum -y install valgrind-ohpc
 RUN yum -y install boost-gnu9-mpich-ohpc
 
 RUN yum -y install python38-pip
-RUN pip3 install matplotlib
 
 RUN yum -y install procps-ng
 RUN yum -y install diffutils
@@ -143,14 +142,18 @@ RUN . /etc/profile.d/lmod.sh \
 RUN rm /tmp/doxygen-1.9.5.src.tar.gz
 
 RUN yum -y install graphviz
-RUN pip3 install gcovr
+RUN . /etc/profile.d/lmod.sh \ && pip3 install gcovr
 
 RUN yum -y install gsl-gnu9-ohpc
 
-RUN pip3 install "pybind11[global]"
-
 RUN . /etc/profile.d/lmod.sh \
-    && pip3 install mpi4py
+    && pip3 install "numpy<2" \
+    && pip3 install scipy \
+    && pip3 install matplotlib \
+    && pip3 install h5py \
+    && pip3 install "mpi4py<4"
+
+RUN pip3 install "pybind11[global]"
 
 # Register new libs installed into /usr/local/lib with linker
 RUN echo "/usr/local/lib" > /etc/ld.so.conf.d/class.conf

src/M2ulPhyS2Boltzmann.cpp

@@ -36,7 +36,6 @@
 #include "M2ulPhyS.hpp"
 #include "tps2Boltzmann.hpp"
 
-// CPU version (just for starting up)
 void M2ulPhyS::push(TPS::Tps2Boltzmann &interface) {
   assert(interface.IsInitialized());
 

src/cycle_avg_joule_coupling.cpp

@@ -348,7 +348,23 @@ void CycleAvgJouleCoupling::interpElectricFieldFromEMToFlow() {
   if (flow_fespace->IsDGSpace()) {
     efieldR_->SetFromTrueDofs(interp_vals);
   } else {
-    assert(false);
+    Array<int> vdofs;
+    Vector elem_dof_vals;
+    int n0 = 0;
+    const int NE = flow_solver_->getMesh()->GetNE();
+    for (int i = 0; i < NE; i++) {
+      flow_fespace->GetElementDofs(i, vdofs);
+      const int nsp = flow_fespace->GetFE(i)->GetNodes().GetNPoints();
+      assert(nsp == vdofs.Size());
+      elem_dof_vals.SetSize(nsp);
+      for (int j = 0; j < nsp; j++) {
+        elem_dof_vals(j) = interp_vals(n0 + j);
+      }
+      efieldR_->SetSubVector(vdofs, elem_dof_vals);
+      n0 += nsp;
+    }
+    efieldR_->SetTrueVector();
+    efieldR_->SetFromTrueVector();
   }
   efieldR_->HostRead();
 
@@ -357,7 +373,23 @@ void CycleAvgJouleCoupling::interpElectricFieldFromEMToFlow() {
   if (flow_fespace->IsDGSpace()) {
     efieldI_->SetFromTrueDofs(interp_vals);
   } else {
-    assert(false);
+    Array<int> vdofs;
+    Vector elem_dof_vals;
+    int n0 = 0;
+    const int NE = flow_solver_->getMesh()->GetNE();
+    for (int i = 0; i < NE; i++) {
+      flow_fespace->GetElementDofs(i, vdofs);
+      const int nsp = flow_fespace->GetFE(i)->GetNodes().GetNPoints();
+      assert(nsp == vdofs.Size());
+      elem_dof_vals.SetSize(nsp);
+      for (int j = 0; j < nsp; j++) {
+        elem_dof_vals(j) = interp_vals(n0 + j);
+      }
+      efieldI_->SetSubVector(vdofs, elem_dof_vals);
+      n0 += nsp;
+    }
+    efieldI_->SetTrueVector();
+    efieldI_->SetFromTrueVector();
   }
   efieldI_->HostRead();
 #else

src/loMach.hpp

@@ -57,6 +57,7 @@
 // forward-declaration for Tps support class
 namespace TPS {
 class Tps;
+class Tps2Boltzmann;
 }
 
 #include <tps_config.h>
@@ -258,6 +259,9 @@ class LoMachSolver : public TPS::PlasmaSolver {
   void evaluatePlasmaConductivityGF() override { thermo_->evaluatePlasmaConductivityGF(); }
 
   mfem::ParGridFunction *getJouleHeatingGF() override { return thermo_->getJouleHeatingGF(); }
+
+  void push(TPS::Tps2Boltzmann &interface) override { thermo_->push(interface); }
+  void fetch(TPS::Tps2Boltzmann &interface) override { thermo_->fetch(interface); }
 };
 
 #endif  // LOMACH_HPP_

src/pytps/__init__.py

@@ -0,0 +1,5 @@
+from .tabulated_reaction_rates import TabulatedSolver
+from .utils import master_print, resolve_runFile, libtps
+
+
+

src/pytps/tabulated_reaction_rates.py

@@ -0,0 +1,64 @@
+import h5py
+import scipy
+import scipy.interpolate
+import numpy as np
+from .utils import master_print, libtps
+
+
+class TabulatedSolver:
+    def __init__(self, comm, config):
+        self.comm = comm
+        self.config = config
+        self.species_densities = None
+        self.efield = None
+        self.heavy_temperature = None
+
+        self.tables = self._read_tables()
+        self.rates = []
+
+    def _findPythonReactions(self):
+        filenames = []
+        nreactions = self.config.getint("reactions","number_of_reactions",fallback=0)
+        for ir in range(nreactions):
+            sublist = self.config["reactions/reaction{0:d}".format(ir+1)]
+            rtype = sublist["model"]
+            if rtype == "bte":
+                filenames.append(sublist["tabulated/filename"].strip("'"))
+
+        return filenames
+
+    def _read_tables(self):
+        filenames = self._findPythonReactions()
+
+        tables = []
+        for filename in filenames:
+            with h5py.File(filename, 'r') as fid:
+                Tcoeff = fid['table'][:]
+
+            tables.append(scipy.interpolate.interp1d(Tcoeff[:,0], Tcoeff[:,1], kind='linear',
+                                                     bounds_error=False, fill_value='extrapolate'))
+
+        return tables
+
+    def fetch(self, interface):
+        n_reactions =interface.nComponents(libtps.t2bIndex.ReactionRates)
+        for r in range(n_reactions):
+            master_print(self.comm,"Reaction ", r+1, ": ", interface.getReactionEquation(r))
+        self.species_densities = np.array(interface.HostRead(libtps.t2bIndex.SpeciesDensities), copy=False)
+        self.efield = np.array(interface.HostRead(libtps.t2bIndex.ElectricField), copy=False)
+        self.heavy_temperature = np.array(interface.HostRead(libtps.t2bIndex.HeavyTemperature), copy=False)
+
+        efieldAngularFreq = interface.EfieldAngularFreq()
+        master_print(self.comm,"Electric field angular frequency: ", efieldAngularFreq)
+
+    def solve(self):
+        self.rates = []
+        for table in self.tables:
+            self.rates.append(table(self.heavy_temperature))
+
+    def push(self, interface):
+        rates =  np.array(interface.HostWrite(libtps.t2bIndex.ReactionRates), copy=False)
+        offset = 0
+        for rate in self.rates:
+            rates[offset:offset+rate.shape[0]] = rate
+            offset = offset+rate.shape[0]

src/pytps/utils.py

@@ -0,0 +1,24 @@
+from mpi4py import MPI
+import os
+import sys
+
+def master_print(comm: MPI.Comm, *args, **kwargs) -> None:
+    if comm.rank == 0:
+        print(*args, **kwargs)
+
+def resolve_runFile(cl_args):
+    ini_name = ''
+    if '-run' in cl_args:
+        ini_name = cl_args[cl_args.index('-run') + 1 ]
+    elif '--runFile' in sys.argv:
+        ini_name = cl_args[cl_args.index('--runFile') + 1 ]
+    else:
+        print("Could not parse command line in python. GOOD BYE!")
+        exit(-1)
+    return ini_name
+
+
+base = os.path.abspath(os.path.dirname( sys.argv[0] ) )
+print(os.path.join(base, ".libs"))
+sys.path.append(os.path.join(base, ".libs"))
+import libtps

src/reactingFlow.cpp

@@ -41,6 +41,7 @@
 #include "loMach.hpp"
 #include "loMach_options.hpp"
 #include "radiation.hpp"
+#include "tps2Boltzmann.hpp"
 
 using namespace mfem;
 using namespace mfem::common;
@@ -363,6 +364,11 @@ ReactingFlow::ReactingFlow(mfem::ParMesh *pmesh, LoMachOptions *loMach_opts, tem
       tpsP_->getRequiredInput((basepath + "/radius").c_str(), R);
       rxnModelParamsHost.push_back(Vector({R}));
 
+    } else if (model == "bte") {
+      reactionModels[r - 1] = GRIDFUNCTION_RXN;
+      int index;
+      tpsP_->getRequiredInput((basepath + "/bte/index").c_str(), index);
+      chemistryInput_.reactionInputs[r - 1].indexInput = index;
     } else {
       grvy_printf(GRVY_ERROR, "\nUnknown reaction_model -> %s", model.c_str());
       exit(ERROR);
@@ -421,7 +427,7 @@ ReactingFlow::ReactingFlow(mfem::ParMesh *pmesh, LoMachOptions *loMach_opts, tem
           equilibriumConstantParams[p + r * gpudata::MAXCHEMPARAMS];
     }
 
-    if (reactionModels[r] != TABULATED_RXN) {
+    if (reactionModels[r] != TABULATED_RXN && reactionModels[r] != GRIDFUNCTION_RXN) {
       assert(rxn_param_idx < rxnModelParamsHost.size());
       chemistryInput_.reactionInputs[r].modelParams = rxnModelParamsHost[rxn_param_idx].Read();
       rxn_param_idx += 1;
@@ -2835,6 +2841,71 @@ double species_uniform(const Vector &coords, double t) {
   return yn;
 }
 
+void ReactingFlow::push(TPS::Tps2Boltzmann &interface) {
+  assert(interface.IsInitialized());
+
+  // const int nscalardofs(vfes_->GetNDofs());
+
+  mfem::ParGridFunction *species =
+      new mfem::ParGridFunction(&interface.NativeFes(TPS::Tps2Boltzmann::Index::SpeciesDensities));
+
+  mfem::Vector speciesInt(sDofInt_*interface.Nspecies());
+  double *species_data = speciesInt.HostWrite();
+
+  const double *dataRho = rn_.HostRead();
+  const double *dataY = Yn_.HostRead();
+
+  double state_local[gpudata::MAXEQUATIONS];
+  double species_local[gpudata::MAXSPECIES];
+
+  for (int i = 0; i < gpudata::MAXEQUATIONS; ++i)
+    state_local[i] = 0.;
+
+  for (int i = 0; i < gpudata::MAXSPECIES; ++i)
+    species_local[i] = 0.;
+
+  for (int i = 0; i < sDofInt_; i++) {
+    state_local[0] = dataRho[i];
+    for (int asp = 0; asp < nActiveSpecies_; asp++)
+      state_local[dim_ + 2 + asp] = dataRho[i]*dataY[i+asp*sDofInt_];
+    mixture_->computeNumberDensities(state_local, species_local);
+
+    for (int sp = 0; sp < interface.Nspecies(); sp++)
+      species_data[i + sp * sDofInt_] = AVOGADRONUMBER * species_local[sp];
+  }
+
+  species->SetFromTrueDofs(speciesInt);
+  Tn_gf_.SetFromTrueDofs(Tn_);
+  interface.interpolateFromNativeFES(*species, TPS::Tps2Boltzmann::Index::SpeciesDensities);
+  interface.interpolateFromNativeFES(Tn_gf_, TPS::Tps2Boltzmann::Index::HeavyTemperature);
+  interface.interpolateFromNativeFES(Tn_gf_, TPS::Tps2Boltzmann::Index::ElectronTemperature);
+
+  interface.setTimeStep(this->dt_);
+  interface.setCurrentTime(this->time_);
+
+  delete species;
+}
+
+void ReactingFlow::fetch(TPS::Tps2Boltzmann &interface) {
+  // NOTE: Chemistry is only address by truedofs indexes.
+  mfem::ParFiniteElementSpace *reaction_rates_fes(&(interface.NativeFes(TPS::Tps2Boltzmann::Index::ReactionRates)));
+  externalReactionRates_gf_.reset(new mfem::ParGridFunction(reaction_rates_fes));
+  interface.interpolateToNativeFES(*externalReactionRates_gf_, TPS::Tps2Boltzmann::Index::ReactionRates);
+  externalReactionRates_gf_->GetTrueDofs(externalReactionRates_);
+  int size = sDofInt_;
+#if defined(_CUDA_) || defined(_HIP_)
+  const double *data(externalReactionRates_.Read());
+  // int size(externalReactionRates_->FESpace()->GetNDofs());
+  assert(externalReactionRates_gf_->FESpace()->GetOrdering() == mfem::Ordering::byNODES);
+  gpu::deviceSetChemistryReactionData<<<1, 1>>>(data, size, chemistry_);
+#else
+  // chemistry_->setGridFunctionRates(*externalReactionRates_gf_);
+  const double *data(externalReactionRates_.HostRead());
+  chemistry_->setRates(data, size);
+#endif
+}
+
+
 #if 0
 void ReactingFlow::uniformInlet() {
   int myRank;

src/reactingFlow.hpp

@@ -36,6 +36,7 @@
 // forward-declaration for Tps support class
 namespace TPS {
 class Tps;
+class Tps2Boltzmann;
 }
 
 #include <hdf5.h>
@@ -102,6 +103,9 @@ class ReactingFlow : public ThermoChemModelBase {
   PerfectMixture *mixture_ = NULL;
   ArgonMixtureTransport *transport_ = NULL;
   Chemistry *chemistry_ = NULL;
+  // External reaction rates when chemistry is implemented using the BTE option
+  std::unique_ptr<ParGridFunction> externalReactionRates_gf_;  // Has repeated interface dofs.
+  Vector externalReactionRates_;  // Only true data
 
   std::vector<std::string> speciesNames_;
   std::map<std::string, int> atomMap_;
@@ -415,5 +419,8 @@ class ReactingFlow : public ThermoChemModelBase {
 
   void evalSubstepNumber();
   void readTableWrapper(std::string inputPath, TableInput &result);
+
+  void push(TPS::Tps2Boltzmann &interface) override;
+  void fetch(TPS::Tps2Boltzmann &interface) override;
 };
 #endif  // REACTINGFLOW_HPP_

src/reaction.cpp

@@ -88,8 +88,8 @@ MFEM_HOST_DEVICE GridFunctionReaction::GridFunctionReaction(int comp)
 MFEM_HOST_DEVICE GridFunctionReaction::~GridFunctionReaction() {}
 
 MFEM_HOST_DEVICE void GridFunctionReaction::setData(const double *data, int size) {
-  data_ = data + comp_ * size_;
   size_ = size;
+  data_ = data + comp_ * size_;
 }
 
 void GridFunctionReaction::setGridFunction(const mfem::GridFunction &f) {
@@ -103,6 +103,7 @@ void GridFunctionReaction::setGridFunction(const mfem::GridFunction &f) {
 #endif
 }
 
+
 MFEM_HOST_DEVICE double GridFunctionReaction::computeRateCoefficient([[maybe_unused]] const double &T_h,
                                                                      [[maybe_unused]] const double &T_e,
                                                                      const int &dofindex,

src/thermo_chem_base.hpp

@@ -39,6 +39,7 @@
 
 namespace TPS {
 class Tps;
+class Tps2Boltzmann;
 }
 
 class IODataOrganizer;
@@ -197,6 +198,16 @@ class ThermoChemModelBase {
     std::cout << "ERROR: " << __func__ << " remains unimplemented" << std::endl;
     exit(1);
   }
+
+  virtual void push(TPS::Tps2Boltzmann &interface) {
+    std::cout << "ERROR: " << __func__ << " remains unimplemented" << std::endl;
+    exit(1);
+  }
+
+  virtual void fetch(TPS::Tps2Boltzmann &interface) {
+    std::cout << "ERROR: " << __func__ << " remains unimplemented" << std::endl;
+    exit(1);
+  }
 };
 
 /**