Initial low Mach reacting flow capabilities

src/Makefile.am

@@ -40,7 +40,8 @@ cxx_sources          = averaging.cpp faceGradientIntegration.cpp M2ulPhyS.cpp rh
                        pybindings.cpp gslib_interpolator.cpp loMach.cpp loMachIO.cpp loMach_options.cpp calorically_perfect.cpp \
 	               thermo_chem_base.cpp split_flow_base.cpp tomboulides.cpp dirichlet_bc_helper.cpp \
                        sponge_base.cpp geometricSponge.cpp gaussianInterpExtData.cpp turb_model_base.cpp \
-                       algebraicSubgridModels.cpp algebraic_rans.cpp mesh_base.cpp lte_thermo_chem.cpp cases.cpp
+                       algebraicSubgridModels.cpp algebraic_rans.cpp mesh_base.cpp lte_thermo_chem.cpp cases.cpp \
+                       reactingFlow.cpp
 
 mfem_extra_sources   = ../utils/mfem_extras/pfem_extras.cpp
 
@@ -57,7 +58,7 @@ headers              = averaging.hpp equation_of_state.hpp transport_properties.
 		       thermo_chem_base.hpp split_flow_base.hpp tomboulides.hpp dirichlet_bc_helper.hpp \
                        sponge_base.hpp geometricSponge.hpp externalData_base.hpp gaussianInterpExtData.hpp \
                        turb_model_base.hpp algebraicSubgridModels.hpp algebraic_rans.hpp mesh_base.hpp \
-                       lte_thermo_chem.hpp cases.hpp
+                       lte_thermo_chem.hpp cases.hpp reactingFlow.hpp
 
 mfem_extra_headers   = ../utils/mfem_extras/pfem_extras.hpp
 

src/argon_transport.cpp

@@ -334,10 +334,12 @@ MFEM_HOST_DEVICE double ArgonMinimalTransport::computeThirdOrderElectronThermalC
          (L11 - L12 * L12 / L22);
 }
 
-void ArgonMinimalTransport::computeMixtureAverageDiffusivity(const Vector &state, Vector &diffusivity) {
+void ArgonMinimalTransport::computeMixtureAverageDiffusivity(const Vector &state, const Vector &Efield,
+                                                             Vector &diffusivity, bool unused) {
   diffusivity.SetSize(3);
   diffusivity = 0.0;
-  computeMixtureAverageDiffusivity(&state[0], &diffusivity[0]);
+  computeMixtureAverageDiffusivity(&state[0], &Efield[0], &diffusivity[0], unused);
+
   // Vector primitiveState(num_equation);
   // mixture->GetPrimitivesFromConservatives(state, primitiveState);
   //
@@ -375,8 +377,8 @@ void ArgonMinimalTransport::computeMixtureAverageDiffusivity(const Vector &state
   // CurtissHirschfelder(X_sp, Y_sp, binaryDiff, diffusivity);
 }
 
-MFEM_HOST_DEVICE void ArgonMinimalTransport::computeMixtureAverageDiffusivity(const double *state,
-                                                                              double *diffusivity) {
+MFEM_HOST_DEVICE void ArgonMinimalTransport::computeMixtureAverageDiffusivity(const double *state, const double *Efield,
+                                                                              double *diffusivity, bool unused) {
   double primitiveState[gpudata::MAXEQUATIONS];
   mixture->GetPrimitivesFromConservatives(state, primitiveState);
 
@@ -550,6 +552,46 @@ MFEM_HOST_DEVICE void ArgonMinimalTransport::GetViscosities(const double *conser
   }
 }
 
+/*
+MFEM_HOST_DEVICE void ArgonMinimalTransport::GetThermalConductivities(const double *conserved, const double *primitive,
+double *kappa) { double n_sp[gpudata::MAXSPECIES], X_sp[gpudata::MAXSPECIES], Y_sp[gpudata::MAXSPECIES];
+  mixture->computeSpeciesPrimitives(conserved, X_sp, Y_sp, n_sp);
+  // double nTotal = 0.0;
+  // for (int sp = 0; sp < numSpecies; sp++) nTotal += n_sp[sp];
+  collisionInputs collInputs = computeCollisionInputs(primitive, n_sp);
+
+  double speciesViscosity[gpudata::MAXSPECIES], speciesHvyThrmCnd[gpudata::MAXSPECIES];
+  for (int sp = 0; sp < numSpecies; sp++) {
+    if (sp == electronIndex_) {
+      speciesViscosity[sp] = 0.0;
+      speciesHvyThrmCnd[sp] = 0.0;
+      continue;
+    }
+    speciesViscosity[sp] =
+        viscosityFactor_ * sqrt(mw_[sp] * collInputs.Th) / collisionIntegral(sp, sp, 2, 2, collInputs);
+    speciesHvyThrmCnd[sp] = speciesViscosity[sp] * kOverEtaFactor_ / mw_[sp];
+  }
+  // transportBuffer[FluxTrns::VISCOSITY] = linearAverage(X_sp, speciesViscosity);
+  // transportBuffer[FluxTrns::HEAVY_THERMAL_CONDUCTIVITY] = linearAverage(X_sp, speciesHvyThrmCnd);
+  //transportBuffer[FluxTrns::BULK_VISCOSITY] = 0.0;
+  kappa[0]  = linearAverage(X_sp, speciesHvyThrmCnd);
+
+  if (thirdOrderkElectron_) {
+    // transportBuffer[FluxTrns::ELECTRON_THERMAL_CONDUCTIVITY] = computeThirdOrderElectronThermalConductivity(X_sp,
+collInputs); kappa[1] = computeThirdOrderElectronThermalConductivity(X_sp, collInputs);
+
+  } else {
+    // transportBuffer[FluxTrns::ELECTRON_THERMAL_CONDUCTIVITY] =
+    //    viscosityFactor_ * kOverEtaFactor_ * sqrt(collInputs.Te / mw_[electronIndex_]) * X_sp[electronIndex_] /
+    //    collisionIntegral(electronIndex_, electronIndex_, 2, 2, collInputs);
+   kappa[1] = viscosityFactor_ * kOverEtaFactor_ * sqrt(collInputs.Te / mw_[electronIndex_]) * X_sp[electronIndex_] /
+              collisionIntegral(electronIndex_, electronIndex_, 2, 2, collInputs);
+
+  }
+
+}
+*/
+
 //////////////////////////////////////////////////////
 //////// Argon Mixture Transport
 //////////////////////////////////////////////////////
@@ -1016,3 +1058,155 @@ MFEM_HOST_DEVICE void ArgonMixtureTransport::GetViscosities(const double *conser
 
   return;
 }
+
+MFEM_HOST_DEVICE void ArgonMixtureTransport::GetThermalConductivities(const double *conserved, const double *primitive,
+                                                                      double *kappa) {
+  double n_sp[gpudata::MAXSPECIES], X_sp[gpudata::MAXSPECIES], Y_sp[gpudata::MAXSPECIES];
+  mixture->computeSpeciesPrimitives(conserved, X_sp, Y_sp, n_sp);
+  // double nTotal = 0.0;
+  // for (int sp = 0; sp < numSpecies; sp++) nTotal += n_sp[sp];
+  collisionInputs collInputs = computeCollisionInputs(primitive, n_sp);
+
+  double speciesViscosity[gpudata::MAXSPECIES], speciesHvyThrmCnd[gpudata::MAXSPECIES];
+  for (int sp = 0; sp < numSpecies; sp++) {
+    if (sp == electronIndex_) {
+      speciesViscosity[sp] = 0.0;
+      speciesHvyThrmCnd[sp] = 0.0;
+      continue;
+    }
+    speciesViscosity[sp] =
+        viscosityFactor_ * sqrt(mw_[sp] * collInputs.Th) / collisionIntegral(sp, sp, 2, 2, collInputs);
+    speciesHvyThrmCnd[sp] = speciesViscosity[sp] * kOverEtaFactor_ / mw_[sp];
+  }
+  // transportBuffer[FluxTrns::VISCOSITY] = linearAverage(X_sp, speciesViscosity);
+  // transportBuffer[FluxTrns::HEAVY_THERMAL_CONDUCTIVITY] = linearAverage(X_sp, speciesHvyThrmCnd);
+  // transportBuffer[FluxTrns::BULK_VISCOSITY] = 0.0;
+  kappa[0] = linearAverage(X_sp, speciesHvyThrmCnd);
+
+  if (thirdOrderkElectron_) {
+    // transportBuffer[FluxTrns::ELECTRON_THERMAL_CONDUCTIVITY] =
+    //   computeThirdOrderElectronThermalConductivity(X_sp, collInputs);
+    kappa[1] = computeThirdOrderElectronThermalConductivity(X_sp, collInputs);
+
+  } else {
+    // transportBuffer[FluxTrns::ELECTRON_THERMAL_CONDUCTIVITY] =
+    //    viscosityFactor_ * kOverEtaFactor_ * sqrt(collInputs.Te / mw_[electronIndex_]) * X_sp[electronIndex_] /
+    //    collisionIntegral(electronIndex_, electronIndex_, 2, 2, collInputs);
+    kappa[1] = viscosityFactor_ * kOverEtaFactor_ * sqrt(collInputs.Te / mw_[electronIndex_]) * X_sp[electronIndex_] /
+               collisionIntegral(electronIndex_, electronIndex_, 2, 2, collInputs);
+  }
+}
+
+void ArgonMixtureTransport::computeMixtureAverageDiffusivity(const Vector &state, const Vector &Efield,
+                                                             Vector &diffusivity, bool unused) {
+  diffusivity.SetSize(3);
+  diffusivity = 0.0;
+  computeMixtureAverageDiffusivity(&state[0], &Efield[0], &diffusivity[0], unused);
+}
+
+/*
+MFEM_HOST_DEVICE void ArgonMixtureTransport::computeMixtureAverageDiffusivity(const double *state, double *diffusivity)
+{ double primitiveState[gpudata::MAXEQUATIONS]; mixture->GetPrimitivesFromConservatives(state, primitiveState);
+
+  double n_sp[3], X_sp[3], Y_sp[3];
+  mixture->computeSpeciesPrimitives(state, X_sp, Y_sp, n_sp);
+
+  double nTotal = 0.0;
+  for (int sp = 0; sp < numSpecies; sp++) nTotal += n_sp[sp];
+
+  double Te = (twoTemperature_) ? primitiveState[num_equation - 1] : primitiveState[nvel_ + 1];
+  double Th = primitiveState[nvel_ + 1];
+
+  // Add Xeps to avoid zero number density case.
+  double nOverT = (n_sp[electronIndex_] + Xeps_) / Te + (n_sp[ionIndex_] + Xeps_) / Th;
+  // double debyeLength = sqrt(debyeFactor_ / AVOGADRONUMBER / nOverT);
+  double debyeLength;
+  debyeLength = debyeFactor_;
+  debyeLength /= AVOGADRONUMBER;
+  debyeLength /= nOverT;
+  debyeLength = std::sqrt(debyeLength);
+  double debyeCircle = PI_ * debyeLength * debyeLength;
+
+  double nondimTe = debyeLength * 4.0 * PI_ * debyeFactor_ * Te;
+  // double nondimTh = debyeLength * 4.0 * PI_ * debyeFactor_ * Th;
+
+  double binaryDiff[3 *3];
+  for (int i = 0; i < 3 * 3; i++) binaryDiff[i] = 0.0;
+  binaryDiff[electronIndex_ + neutralIndex_ * numSpecies] =
+    diffusivityFactor_ * std::sqrt(Te / getMuw(electronIndex_, neutralIndex_)) / nTotal / collision::argon::eAr11(Te);
+  binaryDiff[neutralIndex_ + electronIndex_ * numSpecies] = binaryDiff[electronIndex_ + neutralIndex_ * numSpecies];
+  binaryDiff[neutralIndex_ + ionIndex_ * numSpecies] =
+    diffusivityFactor_ * std::sqrt(Th / getMuw(neutralIndex_, ionIndex_)) / nTotal / collision::argon::ArAr1P11(Th);
+  binaryDiff[ionIndex_ + neutralIndex_ * numSpecies] = binaryDiff[neutralIndex_ + ionIndex_ * numSpecies];
+  binaryDiff[electronIndex_ + ionIndex_ * numSpecies] = diffusivityFactor_ *
+    std::sqrt(Te / getMuw(ionIndex_, electronIndex_)) / nTotal /
+                                                        (collision::charged::att11(nondimTe) * debyeCircle);
+  binaryDiff[ionIndex_ + electronIndex_ * numSpecies] = binaryDiff[electronIndex_ + ionIndex_ * numSpecies];
+
+  // diffusivity.SetSize(3);
+  for (int sp = 0; sp < 3; sp++) diffusivity[sp] = 0.0;
+  CurtissHirschfelder(X_sp, Y_sp, binaryDiff, diffusivity);
+  for (int sp = 0; sp < 3; sp++) {
+  }
+}
+*/
+
+MFEM_HOST_DEVICE void ArgonMixtureTransport::computeMixtureAverageDiffusivity(const double *state, const double *Efield,
+                                                                              double *diffusivity, bool unused) {
+  // double transportBuffer[FluxTrns::NUM_FLUX_TRANS];
+  // for (int p = 0; p < FluxTrns::NUM_FLUX_TRANS; p++) transportBuffer[p] = 0.0;
+
+  double primitiveState[gpudata::MAXEQUATIONS];
+  mixture->GetPrimitivesFromConservatives(state, primitiveState);
+
+  double n_sp[gpudata::MAXSPECIES], X_sp[gpudata::MAXSPECIES], Y_sp[gpudata::MAXSPECIES];
+  mixture->computeSpeciesPrimitives(state, X_sp, Y_sp, n_sp);
+  double nTotal = 0.0;
+  for (int sp = 0; sp < numSpecies; sp++) nTotal += n_sp[sp];
+
+  collisionInputs collInputs = computeCollisionInputs(primitiveState, n_sp);
+
+  /*
+  double speciesViscosity[gpudata::MAXSPECIES];
+  double speciesHvyThrmCnd[gpudata::MAXSPECIES];
+  for (int sp = 0; sp < numSpecies; sp++) {
+    if (sp == electronIndex_) {
+      speciesViscosity[sp] = 0.0;
+      speciesHvyThrmCnd[sp] = 0.0;
+      continue;
+    }
+    speciesViscosity[sp] =
+        viscosityFactor_ * sqrt(mw_[sp] * collInputs.Th) / collisionIntegral(sp, sp, 2, 2, collInputs);
+    speciesHvyThrmCnd[sp] = speciesViscosity[sp] * kOverEtaFactor_ / mw_[sp];
+  }
+  */
+
+  // transportBuffer[FluxTrns::VISCOSITY] = linearAverage(X_sp, speciesViscosity);
+  // transportBuffer[FluxTrns::HEAVY_THERMAL_CONDUCTIVITY] = linearAverage(X_sp, speciesHvyThrmCnd);
+  // transportBuffer[FluxTrns::BULK_VISCOSITY] = 0.0;
+
+  /*
+  if (thirdOrderkElectron_) {
+    transportBuffer[FluxTrns::ELECTRON_THERMAL_CONDUCTIVITY] =
+        computeThirdOrderElectronThermalConductivity(X_sp, collInputs);
+  } else {
+    transportBuffer[FluxTrns::ELECTRON_THERMAL_CONDUCTIVITY] =
+        viscosityFactor_ * kOverEtaFactor_ * sqrt(collInputs.Te / mw_[electronIndex_]) * X_sp[electronIndex_] /
+        collisionIntegral(electronIndex_, electronIndex_, 2, 2, collInputs);
+  }
+  */
+
+  double binaryDiff[gpudata::MAXSPECIES * gpudata::MAXSPECIES];
+  // binaryDiff = 0.0;
+  for (int spI = 0; spI < numSpecies - 1; spI++) {
+    for (int spJ = spI + 1; spJ < numSpecies; spJ++) {
+      double temp = ((spI == electronIndex_) || (spJ == electronIndex_)) ? collInputs.Te : collInputs.Th;
+      binaryDiff[spI + spJ * numSpecies] =
+          diffusivityFactor_ * sqrt(temp / getMuw(spI, spJ)) / nTotal / collisionIntegral(spI, spJ, 1, 1, collInputs);
+      binaryDiff[spJ + spI * numSpecies] = binaryDiff[spI + spJ * numSpecies];
+    }
+  }
+
+  for (int sp = 0; sp < 3; sp++) diffusivity[sp] = 0.0;
+  CurtissHirschfelder(X_sp, Y_sp, binaryDiff, diffusivity);
+}

src/argon_transport.hpp

@@ -68,7 +68,8 @@ class ArgonMinimalTransport : public MolecularTransport {
   const double eps0_ = VACUUMPERMITTIVITY;
   const double qe_ = ELECTRONCHARGE;
   const double R_ = UNIVERSALGASCONSTANT;
-  const double debyeFactor_ = kB_ * eps0_ / qe_ / qe_;
+  // const double debyeFactor_ = kB_ * eps0_ / qe_ / qe_;
+  const double debyeFactor_ = BOLTZMANNCONSTANT * VACUUMPERMITTIVITY / ELECTRONCHARGE / ELECTRONCHARGE;
   // const double PI_ = 4.0 * atan(1.0);
   const double PI_ = PI;
   const double qeOverkB_ = qe_ / kB_;
@@ -124,8 +125,11 @@ class ArgonMinimalTransport : public MolecularTransport {
   MFEM_HOST_DEVICE double computeThirdOrderElectronThermalConductivity(const double *X_sp, const double debyeLength,
                                                                        const double Te, const double nondimTe);
 
-  virtual void computeMixtureAverageDiffusivity(const Vector &state, Vector &diffusivity);
-  MFEM_HOST_DEVICE void computeMixtureAverageDiffusivity(const double *state, double *diffusivity);
+  virtual void computeMixtureAverageDiffusivity(const Vector &state, const Vector &Efield, Vector &diffusivity,
+                                                bool unused);
+  MFEM_HOST_DEVICE virtual void computeMixtureAverageDiffusivity(const double *state, const double *Efield,
+                                                                 double *diffusivity, bool unused);
+  MFEM_HOST_DEVICE void GetThermalConductivities(const double *conserved, const double *primitive, double *kappa);
 
   // These are used to compute third-order electron thermal conductivity based on standard Chapman--Enskog method.
   MFEM_HOST_DEVICE double L11ee(const double *Q2) { return Q2[0]; }
@@ -188,8 +192,13 @@ class ArgonMixtureTransport : public ArgonMinimalTransport {
 
   MFEM_HOST_DEVICE double computeThirdOrderElectronThermalConductivity(const double *X_sp,
                                                                        const collisionInputs &collInputs);
-  //
-  // virtual void computeMixtureAverageDiffusivity(const Vector &state, Vector &diffusivity);
+
+  virtual void computeMixtureAverageDiffusivity(const Vector &state, const Vector &Efield, Vector &diffusivity,
+                                                bool unused);
+  MFEM_HOST_DEVICE virtual void computeMixtureAverageDiffusivity(const double *state, const double *Efield,
+                                                                 double *diffusivity, bool unused);
+
+  MFEM_HOST_DEVICE void GetThermalConductivities(const double *conserved, const double *primitive, double *kappa);
 };
 
 #endif  // ARGON_TRANSPORT_HPP_

src/calorically_perfect.cpp

@@ -118,6 +118,7 @@ CaloricallyPerfectThermoChem::CaloricallyPerfectThermoChem(mfem::ParMesh *pmesh,
   tpsP_->getInput("loMach/calperfect/linear-solver-rtol", rtol_, 1e-12);
   tpsP_->getInput("loMach/calperfect/linear-solver-max-iter", max_iter_, 1000);
   tpsP_->getInput("loMach/calperfect/linear-solver-verbosity", pl_solve_, 0);
+  tpsP_->getInput("loMach/calperfect/numerical-integ", numerical_integ_, true);
 }
 
 CaloricallyPerfectThermoChem::~CaloricallyPerfectThermoChem() {
@@ -358,7 +359,7 @@ void CaloricallyPerfectThermoChem::initializeSelf() {
 
   // Wall BCs
   {
-    std::cout << "There are " << pmesh_->bdr_attributes.Max() << " boundary attributes!" << std::endl;
+    if (rank0_) std::cout << "There are " << pmesh_->bdr_attributes.Max() << " boundary attributes" << std::endl;
     Array<int> attr_wall(pmesh_->bdr_attributes.Max());
     attr_wall = 0;
 
@@ -371,7 +372,7 @@ void CaloricallyPerfectThermoChem::initializeSelf() {
       tpsP_->getRequiredInput((basepath + "/type").c_str(), type);
 
       if (type == "viscous_isothermal") {
-        std::cout << "Adding patch = " << patch << " to isothermal wall list!" << std::endl;
+        if (rank0_) std::cout << "Adding patch = " << patch << " to isothermal wall list" << std::endl;
 
         attr_wall = 0;
         attr_wall[patch - 1] = 1;
@@ -381,7 +382,6 @@ void CaloricallyPerfectThermoChem::initializeSelf() {
 
         ConstantCoefficient *Twall_coeff = new ConstantCoefficient();
         Twall_coeff->constant = Twall;
-
         AddTempDirichletBC(Twall_coeff, attr_wall);
 
         ConstantCoefficient *Qt_bc_coeff = new ConstantCoefficient();
@@ -557,7 +557,7 @@ void CaloricallyPerfectThermoChem::initializeOperators() {
   MqInv_->SetMaxIter(max_iter_);
 
   LQ_form_ = new ParBilinearForm(sfes_);
-  auto *lqd_blfi = new DiffusionIntegrator(*thermal_diff_coeff_);
+  auto *lqd_blfi = new DiffusionIntegrator(*thermal_diff_total_coeff_);
   if (numerical_integ_) {
     lqd_blfi->SetIntRule(&ir_di);
   }