Merge branch 'boltzmann-integration' into main

Merging PR #242

src/M2ulPhyS.cpp

@@ -405,7 +405,14 @@ void M2ulPhyS::initVariables() {
     serial_mesh->GetNodes(coordinates);
 
     // Evaluate the distance function
-    evaluateDistanceSerial(*serial_mesh, wall_patch_list, coordinates, *serial_distance);
+    if (!config.read_distance || !config.GetRestartCycle()) {
+      if (rank0_) grvy_printf(ginfo, "Computing distance function\n");
+      evaluateDistanceSerial(*serial_mesh, wall_patch_list, coordinates, *serial_distance);
+    } else {
+      // If distance function is read from restart, this will be overwritten later
+      if (rank0_) grvy_printf(ginfo, "Distance function to be read from restart\n");
+      *serial_distance = 0.0;
+    }
     delete tmp_dfes;
   }
 
@@ -617,6 +624,11 @@ void M2ulPhyS::initVariables() {
 #endif
   initSolutionAndVisualizationVectors();
 
+  if (distance_ != NULL) {
+    ioData.registerIOFamily("Distance function", "/distance", distance_, false, config.read_distance);
+    ioData.registerIOVar("/distance", "distance", 0, config.read_distance);
+  }
+
   average = new Averaging(config.avg_opts_, config.GetOutputName());
   average->registerField(std::string("primitive_state"), Up, true, 1, nvel);
   average->initializeVizForM2ulPhyS(fes, dfes, nvel);
@@ -1076,6 +1088,12 @@ void M2ulPhyS::initIndirectionArrays() {
   // See #199 for more info.
   const int NumBCelems = fes->GetNBE();
 
+  // NB: *Must* call this here, as otherwise some faces are
+  // erroneously included as boundary faces and asserts below may
+  // fail
+  mesh->ExchangeFaceNbrNodes();
+  mesh->ExchangeFaceNbrData();
+
   if (NumBCelems > 0) {
     bdry_face_data.shape.UseDevice(true);
     bdry_face_data.shape.SetSize(NumBCelems * maxIntPoints * maxDofs);
@@ -1119,12 +1137,6 @@ void M2ulPhyS::initIndirectionArrays() {
     FaceElementTransformations *tr;
     // Mesh *mesh = fes->GetMesh();
 
-    // NB: *Must* call this here, as otherwise some faces are
-    // erroneously included as boundary faces and asserts below may
-    // fail
-    mesh->ExchangeFaceNbrNodes();
-    mesh->ExchangeFaceNbrData();
-
     std::vector<int> uniqueElems;
     uniqueElems.clear();
 
@@ -1941,6 +1953,10 @@ void M2ulPhyS::projectInitialSolution() {
 
   initGradUp();
 
+  // Exchange before computing primitives
+  U->ParFESpace()->ExchangeFaceNbrData();
+  U->ExchangeFaceNbrData();
+
   updatePrimitives();
 
   // update pressure grid function
@@ -2662,6 +2678,7 @@ void M2ulPhyS::parseFlowOptions() {
   }
   tpsP->getInput("flow/refinement_levels", config.ref_levels, 0);
   tpsP->getInput("flow/computeDistance", config.compute_distance, false);
+  tpsP->getInput("flow/readDistance", config.read_distance, false);
 
   std::string type;
   tpsP->getInput("flow/sgsModel", type, std::string("none"));
@@ -3298,6 +3315,11 @@ void M2ulPhyS::parseReactionInputs() {
       config.reactionModels[r - 1] = TABULATED_RXN;
       std::string inputPath(basepath + "/tabulated");
       readTable(inputPath, config.chemistryInput.reactionInputs[r - 1].tableInput);
+    } else if (model == "bte") {
+      config.reactionModels[r - 1] = GRIDFUNCTION_RXN;
+      int index;
+      tpsP->getRequiredInput((basepath + "/bte/index").c_str(), index);
+      config.chemistryInput.reactionInputs[r - 1].indexInput = index;
     } else {
       grvy_printf(GRVY_ERROR, "\nUnknown reaction_model -> %s", model.c_str());
       exit(ERROR);
@@ -3416,7 +3438,7 @@ void M2ulPhyS::parseReactionInputs() {
             config.equilibriumConstantParams[p + r * gpudata::MAXCHEMPARAMS];
       }
 
-      if (config.reactionModels[r] != TABULATED_RXN) {
+      if (config.reactionModels[r] == ARRHENIUS || config.reactionModels[r] == HOFFERTLIEN) {
         assert(rxn_param_idx < config.rxnModelParamsHost.size());
         config.chemistryInput.reactionInputs[r].modelParams = config.rxnModelParamsHost[rxn_param_idx].Read();
         rxn_param_idx += 1;
@@ -3976,8 +3998,7 @@ void M2ulPhyS::checkSolverOptions() const {
 }
 
 void M2ulPhyS::updatePrimitives() {
-  // U.V.: should this be U->HostRead() instead? U->HostWrite() does not sync memory before returning the pointer.
-  double *data = U->HostWrite();
+  const double *data = U->HostRead();
   double *dataUp = Up->HostWrite();
   int dof = vfes->GetNDofs();
 
@@ -4175,7 +4196,7 @@ void M2ulPhyS::updateVisualizationVariables() {
       Th = prim[1 + _nvel];
       Te = (in_mix->IsTwoTemperature()) ? prim[_num_equation - 1] : Th;
       double kfwd[gpudata::MAXREACTIONS], kC[gpudata::MAXREACTIONS];
-      in_chem->computeForwardRateCoeffs(Th, Te, kfwd);
+      in_chem->computeForwardRateCoeffs(Th, Te, n, kfwd);
       in_chem->computeEquilibriumConstants(Th, Te, kC);
       // get reaction rates
       double progressRates[gpudata::MAXREACTIONS];

src/M2ulPhyS.hpp

@@ -250,6 +250,9 @@ class M2ulPhyS : public TPS::PlasmaSolver {
   // ParNonlinearForm *gradUp_A;
   GradNonLinearForm *gradUp_A;
 
+  // Auxiliary grid function to store external reaction rates
+  std::unique_ptr<ParGridFunction> externalReactionRates;
+
   // Average handler
   Averaging *average;
 

src/M2ulPhyS2Boltzmann.cpp

@@ -78,9 +78,24 @@ void M2ulPhyS::push(TPS::Tps2Boltzmann &interface) {
   interface.interpolateFromNativeFES(*heavyTemperature, TPS::Tps2Boltzmann::Index::HeavyTemperature);
   interface.interpolateFromNativeFES(*electronTemperature, TPS::Tps2Boltzmann::Index::ElectronTemperature);
 
+  interface.setTimeStep(this->dt);
+  interface.setCurrentTime(this->time);
+
   delete species;
   delete heavyTemperature;
   delete electronTemperature;
 }
 
-void M2ulPhyS::fetch(TPS::Tps2Boltzmann &interface) { return; }
+void M2ulPhyS::fetch(TPS::Tps2Boltzmann &interface) {
+  mfem::ParFiniteElementSpace *reaction_rates_fes(&(interface.NativeFes(TPS::Tps2Boltzmann::Index::ReactionRates)));
+  externalReactionRates.reset(new mfem::ParGridFunction(reaction_rates_fes));
+  interface.interpolateToNativeFES(*externalReactionRates, TPS::Tps2Boltzmann::Index::ReactionRates);
+  #if defined(_CUDA_) || defined(_HIP_)
+    const double * data(externalReactionRates->Read() );
+    int size(externalReactionRates->FESpace()->GetNDofs() );
+    assert(externalReactionRates->FESpace()->GetOrdering() == mfem::Ordering::byNODES);
+    gpu::deviceSetChemistryReactionData<<<1, 1>>>(data, size, chemistry_);
+  #else
+    chemistry_->setGridFunctionRates(*externalReactionRates);
+  #endif
+}

src/chemistry.cpp

@@ -86,6 +86,9 @@ MFEM_HOST_DEVICE Chemistry::Chemistry(GasMixture *mixture, const ChemistryInput
       case TABULATED_RXN: {
         reactions_[r] = new Tabulated(inputs.reactionInputs[r].tableInput);
       } break;
+      case GRIDFUNCTION_RXN: {
+        reactions_[r] = new GridFunctionReaction(inputs.reactionInputs[r].indexInput);
+      } break;
       default:
         printf("Unknown reactionModel.");
         assert(false);
@@ -106,6 +109,25 @@ MFEM_HOST_DEVICE Chemistry::~Chemistry() {
   }
 }
 
+MFEM_HOST_DEVICE void Chemistry::setRates(const double * data, int size) {
+  for (int r = 0; r < numReactions_; r++) {
+    if (reactions_[r]->reactionModel == GRIDFUNCTION_RXN) {
+      GridFunctionReaction *rx = (GridFunctionReaction *)(reactions_[r]);
+      rx->setData(data, size);
+    }
+  }
+}
+
+void Chemistry::setGridFunctionRates(mfem::GridFunction &f) {
+  for (int r = 0; r < numReactions_; r++) {
+    if (reactions_[r]->reactionModel == GRIDFUNCTION_RXN) {
+      GridFunctionReaction *rx = dynamic_cast<GridFunctionReaction *>(reactions_[r]);
+      rx->setGridFunction(f);
+    }
+  }
+}
+
+#if 0
 void Chemistry::computeForwardRateCoeffs(const double &T_h, const double &T_e, Vector &kfwd) {
   kfwd.SetSize(numReactions_);
 
@@ -122,8 +144,10 @@ void Chemistry::computeForwardRateCoeffs(const double &T_h, const double &T_e, V
 
   return;
 }
+#endif
 
-MFEM_HOST_DEVICE void Chemistry::computeForwardRateCoeffs(const double &T_h, const double &T_e, double *kfwd) {
+MFEM_HOST_DEVICE void Chemistry::computeForwardRateCoeffs(const double &T_h, const double &T_e, const int &dofindex,
+                                                          double *kfwd) {
   // kfwd.SetSize(numReactions_);
   for (int r = 0; r < numReactions_; r++) kfwd[r] = 0.0;
 
@@ -132,12 +156,13 @@ MFEM_HOST_DEVICE void Chemistry::computeForwardRateCoeffs(const double &T_h, con
 
   for (int r = 0; r < numReactions_; r++) {
     bool isElectronInvolved = isElectronInvolvedAt(r);
-    kfwd[r] = reactions_[r]->computeRateCoefficient(Thlim, Telim, isElectronInvolved);
+    kfwd[r] = reactions_[r]->computeRateCoefficient(Thlim, Telim, dofindex, isElectronInvolved);
   }
 
   return;
 }
 
+#if 0
 // NOTE: if not detailedBalance, equilibrium constant is returned as zero, though it cannot be used.
 void Chemistry::computeEquilibriumConstants(const double &T_h, const double &T_e, Vector &kC) {
   kC.SetSize(numReactions_);
@@ -159,6 +184,7 @@ void Chemistry::computeEquilibriumConstants(const double &T_h, const double &T_e
 
   return;
 }
+#endif
 
 MFEM_HOST_DEVICE void Chemistry::computeEquilibriumConstants(const double &T_h, const double &T_e, double *kC) {
   for (int r = 0; r < numReactions_; r++) kC[r] = 0.0;

src/chemistry.hpp

@@ -95,13 +95,18 @@ class Chemistry {
 
   MFEM_HOST_DEVICE ~Chemistry();
 
+  // Set the grid function rates for GRIDFUNCTION_RXN reaction types
+  void setGridFunctionRates(mfem::GridFunction &f);
+  MFEM_HOST_DEVICE void setRates(const double * data, int size);
+
   // return Vector of reaction rate coefficients, with the size of numReaction_.
   // WARNING(marc) I have removed "virtual" qualifier here assuming these functions will not
   // change for child classes. Correct if wrong
-  void computeForwardRateCoeffs(const double &T_h, const double &T_e, Vector &kfwd);
-  MFEM_HOST_DEVICE void computeForwardRateCoeffs(const double &T_h, const double &T_e, double *kfwd);
+  // void computeForwardRateCoeffs(const double &T_h, const double &T_e, Vector &kfwd);
+  MFEM_HOST_DEVICE void computeForwardRateCoeffs(const double &T_h, const double &T_e, const int &dofindex,
+                                                 double *kfwd);
 
-  void computeEquilibriumConstants(const double &T_h, const double &T_e, Vector &kC);
+  // void computeEquilibriumConstants(const double &T_h, const double &T_e, Vector &kC);
   MFEM_HOST_DEVICE void computeEquilibriumConstants(const double &T_h, const double &T_e, double *kC);
 
   // return rate coefficients of (reactionIndex)-th reaction. (start from 0)

src/dataStructures.hpp

@@ -74,7 +74,7 @@ enum TransportModel { ARGON_MINIMAL, ARGON_MIXTURE, CONSTANT, LTE_TRANSPORT, MIX
 
 enum ChemistryModel { /* CANTERA, */ NUM_CHEMISTRYMODEL };
 
-enum ReactionModel { ARRHENIUS, HOFFERTLIEN, TABULATED_RXN, NUM_REACTIONMODEL };
+enum ReactionModel { ARRHENIUS, HOFFERTLIEN, TABULATED_RXN, GRIDFUNCTION_RXN, NUM_REACTIONMODEL };
 
 enum RadiationModel { NONE_RAD, NET_EMISSION, NUM_RADIATIONMODEL };
 
@@ -632,6 +632,7 @@ struct ReactionInput {
   TableInput tableInput;
   // NOTE(kevin): with gpu, this pointer is only valid on the device.
   const double *modelParams;
+  int indexInput;
 };
 
 struct ChemistryInput {

src/forcing_terms.cpp

@@ -540,7 +540,6 @@ SpongeZone::SpongeZone(const int &_dim, const int &_num_equation, const int &_or
   ParFiniteElementSpace fes(mesh, fec);
 
   sigma = new ParGridFunction(&fes);
-  *sigma = 0.;
   double *hSigma = sigma->HostWrite();
 
   ParGridFunction coords(&dfes);
@@ -560,6 +559,7 @@ SpongeZone::SpongeZone(const int &_dim, const int &_num_equation, const int &_or
     Vector Xn(dim);
     for (int d = 0; d < dim; d++) Xn[d] = coords[n + d * ndofs];
 
+    hSigma[n] = 0.0;
     if (szData.szType == SpongeZoneType::PLANAR) {
       // distance to the mix-out plane
       double distInit = 0.;

src/gpu_constructor.cpp

@@ -121,5 +121,16 @@ __global__ void freeDeviceRadiation(Radiation *radiation) {
   if (radiation != NULL) radiation->~Radiation();
 }
 
+//---------------------------------------------------
+// And finally device setters
+//---------------------------------------------------
+__global__ void deviceSetGridFunctionReactionData(const double * data, int size, GridFunctionReaction * reaction) {
+  reaction->setData(data, size);
+}
+
+__global__ void deviceSetChemistryReactionData(const double * data, int size, Chemistry * chem) {
+  chem->setRates(data, size);
+}
+
 #endif  // cuda or hip
 }  // namespace gpu

src/gpu_constructor.hpp

@@ -160,6 +160,10 @@ __global__ void freeDeviceChemistry(Chemistry *chem);
 //! Explicit call to Radiation destructor on the device
 __global__ void freeDeviceRadiation(Radiation *radiation);
 
+//! Set the data to a GridFunctionReaction
+__global__ void deviceSetGridFunctionReactionData(const double * data, int size, GridFunctionReaction * reaction);
+__global__ void deviceSetChemistryReactionData(const double * data, int size, Chemistry * chem);
+
 #endif  // cuda or hip
 }  // namespace gpu
 

src/outletBC.cpp

@@ -74,6 +74,12 @@ OutletBC::OutletBC(MPI_Groups *_groupsMPI, Equations _eqSystem, RiemannSolver *_
   hmeanUp[1 + nvel_] = 300.0;  // 101300;
   if (eqSystem == NS_PASSIVE) hmeanUp[num_equation_ - 1] = 0.;
 
+  if (mixture->GetNumActiveSpecies() > 0) {
+    for (int sp = 0; sp < mixture->GetNumActiveSpecies(); sp++) {
+      hmeanUp[nvel_ + 2 + sp] = 0.0;
+    }
+  }
+
   area_ = 0.;
   parallelAreaComputed = false;
 

src/reaction.cpp

@@ -36,27 +36,29 @@ using namespace mfem;
 using namespace std;
 
 MFEM_HOST_DEVICE Arrhenius::Arrhenius(const double &A, const double &b, const double &E)
-    : Reaction(), A_(A), b_(b), E_(E) {}
+    : Reaction(ARRHENIUS), A_(A), b_(b), E_(E) {}
 
 MFEM_HOST_DEVICE double Arrhenius::computeRateCoefficient(const double &T_h, const double &T_e,
+                                                          [[maybe_unused]] const int &dofindex,
                                                           const bool isElectronInvolved) {
   double temp = (isElectronInvolved) ? T_e : T_h;
 
   return A_ * pow(temp, b_) * exp(-E_ / UNIVERSALGASCONSTANT / temp);
 }
 
 MFEM_HOST_DEVICE HoffertLien::HoffertLien(const double &A, const double &b, const double &E)
-    : Reaction(), A_(A), b_(b), E_(E) {}
+    : Reaction(HOFFERTLIEN), A_(A), b_(b), E_(E) {}
 
 MFEM_HOST_DEVICE double HoffertLien::computeRateCoefficient(const double &T_h, const double &T_e,
+                                                            [[maybe_unused]] const int &dofindex,
                                                             const bool isElectronInvolved) {
   double temp = (isElectronInvolved) ? T_e : T_h;
   double tempFactor = E_ / BOLTZMANNCONSTANT / temp;
 
   return A_ * pow(temp, b_) * (tempFactor + 2.0) * exp(-tempFactor);
 }
 
-MFEM_HOST_DEVICE Tabulated::Tabulated(const TableInput &input) : Reaction() {
+MFEM_HOST_DEVICE Tabulated::Tabulated(const TableInput &input) : Reaction(TABULATED_RXN) {
   switch (input.order) {
     case 1: {
       table_ = new LinearTable(input);
@@ -71,7 +73,41 @@ MFEM_HOST_DEVICE Tabulated::Tabulated(const TableInput &input) : Reaction() {
 MFEM_HOST_DEVICE Tabulated::~Tabulated() { delete table_; }
 
 MFEM_HOST_DEVICE double Tabulated::computeRateCoefficient(const double &T_h, const double &T_e,
+                                                          [[maybe_unused]] const int &dofindex,
                                                           const bool isElectronInvolved) {
   double temp = (isElectronInvolved) ? T_e : T_h;
   return table_->eval(temp);
 }
+
+MFEM_HOST_DEVICE GridFunctionReaction::GridFunctionReaction(int comp)
+    : Reaction(GRIDFUNCTION_RXN), data_(nullptr), comp_(comp), size_(0) {}
+
+MFEM_HOST_DEVICE GridFunctionReaction::~GridFunctionReaction() {}
+
+MFEM_HOST_DEVICE void GridFunctionReaction::setData(const double * data, int size) {
+  data_ = data + comp_ * size_;
+  size_ = size;
+}
+
+void GridFunctionReaction::setGridFunction(const mfem::GridFunction & f) {
+  size_ = f.FESpace()->GetNDofs();
+  assert(comp_ < f.FESpace()->GetVDim());
+  assert(f.FESpace()->GetOrdering() == mfem::Ordering::byNODES);
+#if defined(_CUDA_) || defined(_HIP_)
+  data_ = f.Read() + comp_ * size_;
+#else
+  data_ = f.HostRead() + comp_ * size_;
+#endif
+}
+
+MFEM_HOST_DEVICE double GridFunctionReaction::computeRateCoefficient([[maybe_unused]] const double &T_h,
+                                                                     [[maybe_unused]] const double &T_e,
+                                                                     const int &dofindex,
+                                                                     [[maybe_unused]] const bool isElectronInvolved) {
+  if (data_) {
+    assert(dofindex < size_);
+    return data_[dofindex];
+  } else {
+    return 0.;
+  }
+}