Merge dev-implicit-chem-7sp-supg-tested, reimplementing streamwise stability for lomach solvers

src/M2ulPhyS.cpp

@@ -237,7 +237,7 @@ void M2ulPhyS::initMixtureAndTransportModels() {
                                               TPSCommWorld, trans_data, trans_tables);
         if (!success) exit(ERROR);
 
-          // Instantiate LteTransport class
+        // Instantiate LteTransport class
 #if defined(_CUDA_) || defined(_HIP_)
         // Tables from above have host pointers.  Must get device
         // pointers here before instantiating device class
@@ -1784,7 +1784,7 @@ void M2ulPhyS::initSolutionAndVisualizationVectors() {
         // visualizationNames_.push_back(std::string("rxn_rate: " + config.reactionEquations[r]));
       }
     }  // if (config.workFluid != DRY_AIR)
-  }    // if tpsP->isVisualizationMode()
+  }  // if tpsP->isVisualizationMode()
 
   // If mms, add conserved and exact solution.
 #ifdef HAVE_MASA
@@ -2094,7 +2094,7 @@ void M2ulPhyS::solveStep() {
       exit(ERROR);
 #endif
     }  // plane dump
-  }    // step check
+  }  // step check
 
   average->addSample(iter, d_mixture);
 }
@@ -2465,7 +2465,9 @@ void M2ulPhyS::Check_NAN() {
   int dof = vfes->GetNDofs();
 
 #ifdef _GPU_
-  { local_print = M2ulPhyS::Check_NaN_GPU(U, dof * num_equation, loc_print); }
+  {
+    local_print = M2ulPhyS::Check_NaN_GPU(U, dof * num_equation, loc_print);
+  }
   if (local_print > 0) {
     cout << "Found a NaN!" << endl;
   }
@@ -4234,7 +4236,7 @@ void M2ulPhyS::updateVisualizationVariables() {
         dataVis[visualIdxs.rxn + r][n] = progressRates[r];
       }
     }  // if (!isDryAir)
-  }    // for (int n = 0; n < ndofs; n++)
+  }  // for (int n = 0; n < ndofs; n++)
 }
 
 void M2ulPhyS::evaluatePlasmaConductivityGF() {

src/calorically_perfect.cpp

@@ -58,10 +58,12 @@ MFEM_HOST_DEVICE double Sutherland(const double T, const double mu_star, const d
 }
 
 CaloricallyPerfectThermoChem::CaloricallyPerfectThermoChem(mfem::ParMesh *pmesh, LoMachOptions *loMach_opts,
-                                                           temporalSchemeCoefficients &time_coeff, TPS::Tps *tps)
+                                                           temporalSchemeCoefficients &time_coeff,
+                                                           ParGridFunction *gridScale, TPS::Tps *tps)
     : tpsP_(tps), pmesh_(pmesh), time_coeff_(time_coeff) {
   rank0_ = (pmesh_->GetMyRank() == 0);
   order_ = loMach_opts->order;
+  gridScale_gf_ = gridScale;
 
   std::string visc_model;
   tpsP_->getInput("loMach/calperfect/viscosity-model", visc_model, std::string("sutherland"));
@@ -132,6 +134,9 @@ CaloricallyPerfectThermoChem::CaloricallyPerfectThermoChem(mfem::ParMesh *pmesh,
   tpsP_->getInput("loMach/calperfect/msolve-atol", mass_inverse_atol_, default_atol_);
   tpsP_->getInput("loMach/calperfect/msolve-max-iter", mass_inverse_max_iter_, max_iter_);
   tpsP_->getInput("loMach/calperfect/msolve-verbosity", mass_inverse_pl_, pl_solve_);
+
+  // artificial diffusion (SUPG)
+  tpsP_->getInput("loMach/calperfect/streamwise-stabilization", sw_stab_, false);
 }
 
 CaloricallyPerfectThermoChem::~CaloricallyPerfectThermoChem() {
@@ -164,6 +169,20 @@ CaloricallyPerfectThermoChem::~CaloricallyPerfectThermoChem() {
   delete rho_over_dt_coeff_;
   delete rho_coeff_;
 
+  delete umag_coeff_;
+  delete gscale_coeff_;
+  delete visc_coeff_;
+  delete visc_inv_coeff_;
+  delete reh1_coeff_;
+  delete reh2_coeff_;
+  delete Reh_coeff_;
+  delete csupg_coeff_;
+  delete uw1_coeff_;
+  delete uw2_coeff_;
+  delete upwind_coeff_;
+  delete swdiff_coeff_;
+  delete supg_coeff_;
+
   // allocated in initializeSelf
   delete sfes_;
   delete sfec_;
@@ -448,6 +467,32 @@ void CaloricallyPerfectThermoChem::initializeOperators() {
   rhon_next_coeff_ = new GridFunctionCoefficient(&rn_gf_);
   rhou_coeff_ = new ScalarVectorProductCoefficient(*rhon_next_coeff_, *un_next_coeff_);
 
+  // artifical diffusion coefficients
+  if (sw_stab_) {
+    umag_coeff_ = new VectorMagnitudeCoefficient(*un_next_coeff_);
+    gscale_coeff_ = new GridFunctionCoefficient(gridScale_gf_);
+    visc_coeff_ = new GridFunctionCoefficient(&visc_gf_);
+    visc_inv_coeff_ = new PowerCoefficient(*visc_coeff_, -1.0);
+
+    // compute Reh
+    reh1_coeff_ = new ProductCoefficient(*rho_coeff_, *visc_inv_coeff_);
+    reh2_coeff_ = new ProductCoefficient(*reh1_coeff_, *gscale_coeff_);
+    Reh_coeff_ = new ProductCoefficient(*reh1_coeff_, *umag_coeff_);
+
+    // Csupg
+    csupg_coeff_ = new TransformedCoefficient(Reh_coeff_, csupgFactor);
+
+    // compute upwind magnitude
+    uw1_coeff_ = new ProductCoefficient(*rho_coeff_, *csupg_coeff_);
+    uw2_coeff_ = new ProductCoefficient(*uw1_coeff_, *gscale_coeff_);
+    upwind_coeff_ = new ProductCoefficient(*uw2_coeff_, *umag_coeff_);
+
+    // streamwise diffusion direction
+    swdiff_coeff_ = new TransformedMatrixVectorCoefficient(un_next_coeff_, &streamwiseTensor);
+
+    supg_coeff_ = new ScalarMatrixProductCoefficient(*upwind_coeff_, *swdiff_coeff_);
+  }
+
   At_form_ = new ParBilinearForm(sfes_);
   auto *at_blfi = new ConvectionIntegrator(*rhou_coeff_);
   if (numerical_integ_) {
@@ -483,6 +528,7 @@ void CaloricallyPerfectThermoChem::initializeOperators() {
   MsRho_form_->FormSystemMatrix(empty, MsRho_);
   if (rank0_) std::cout << "CaloricallyPerfectThermoChem MsRho operator set" << endl;
 
+  // Helmholtz
   Ht_form_ = new ParBilinearForm(sfes_);
   auto *hmt_blfi = new MassIntegrator(*rho_over_dt_coeff_);
   auto *hdt_blfi = new DiffusionIntegrator(*thermal_diff_total_coeff_);
@@ -491,6 +537,15 @@ void CaloricallyPerfectThermoChem::initializeOperators() {
     hmt_blfi->SetIntRule(&ir_di);
     hdt_blfi->SetIntRule(&ir_di);
   }
+  // SUPG
+  if (sw_stab_) {
+    auto *sdt_blfi = new DiffusionIntegrator(*supg_coeff_);
+    // SUPG diffusion
+    // if (numerical_integ_) {
+    //   sdt_blfi->SetIntRule(&ir_di);
+    // }
+    Ht_form_->AddDomainIntegrator(sdt_blfi);
+  }
   Ht_form_->AddDomainIntegrator(hmt_blfi);
   Ht_form_->AddDomainIntegrator(hdt_blfi);
   Ht_form_->Assemble();
@@ -574,6 +629,18 @@ void CaloricallyPerfectThermoChem::initializeOperators() {
     lqd_blfi->SetIntRule(&ir_di);
   }
   LQ_form_->AddDomainIntegrator(lqd_blfi);
+
+  // DiffusionIntegrator *slqd_blfi;
+  if (sw_stab_) {
+    // slqd_blfi = new DiffusionIntegrator(*supg_coeff_);
+    auto *slqd_blfi = new DiffusionIntegrator(*supg_coeff_);
+    // SUPG diffusion
+    // if (numerical_integ_) {
+    //   slqd_blfi->SetIntRule(&ir_di);
+    // }
+    LQ_form_->AddDomainIntegrator(slqd_blfi);
+  }
+
   if (partial_assembly_) {
     LQ_form_->SetAssemblyLevel(AssemblyLevel::PARTIAL);
   }

src/calorically_perfect.hpp

@@ -48,6 +48,7 @@ class Tps;
 #include "io.hpp"
 #include "thermo_chem_base.hpp"
 #include "tps_mfem_wrap.hpp"
+#include "utils.hpp"
 
 using VecFuncT = void(const Vector &x, double t, Vector &u);
 using ScalarFuncT = double(const Vector &x, double t);
@@ -136,6 +137,9 @@ class CaloricallyPerfectThermoChem : public ThermoChemModelBase {
   // Initial temperature value (if constant IC)
   double T_ic_;
 
+  // streamwise-stabilization
+  bool sw_stab_;
+
   // FEM related fields and objects
 
   // Scalar \f$H^1\f$ finite element collection.
@@ -155,6 +159,8 @@ class CaloricallyPerfectThermoChem : public ThermoChemModelBase {
   ParGridFunction R0PM0_gf_;
   ParGridFunction Qt_gf_;
 
+  ParGridFunction *gridScale_gf_ = nullptr;
+
   // ParGridFunction *buffer_tInlet_ = nullptr;
   GridFunctionCoefficient *temperature_bc_field_ = nullptr;
 
@@ -171,6 +177,20 @@ class CaloricallyPerfectThermoChem : public ThermoChemModelBase {
   GridFunctionCoefficient *rho_over_dt_coeff_ = nullptr;
   GridFunctionCoefficient *rho_coeff_ = nullptr;
 
+  VectorMagnitudeCoefficient *umag_coeff_ = nullptr;
+  GridFunctionCoefficient *gscale_coeff_ = nullptr;
+  GridFunctionCoefficient *visc_coeff_ = nullptr;
+  PowerCoefficient *visc_inv_coeff_ = nullptr;
+  ProductCoefficient *reh1_coeff_ = nullptr;
+  ProductCoefficient *reh2_coeff_ = nullptr;
+  ProductCoefficient *Reh_coeff_ = nullptr;
+  TransformedCoefficient *csupg_coeff_ = nullptr;
+  ProductCoefficient *uw1_coeff_ = nullptr;
+  ProductCoefficient *uw2_coeff_ = nullptr;
+  ProductCoefficient *upwind_coeff_ = nullptr;
+  TransformedMatrixVectorCoefficient *swdiff_coeff_ = nullptr;
+  ScalarMatrixProductCoefficient *supg_coeff_ = nullptr;
+
   // operators and solvers
   ParBilinearForm *At_form_ = nullptr;
   ParBilinearForm *Ms_form_ = nullptr;
@@ -233,7 +253,7 @@ class CaloricallyPerfectThermoChem : public ThermoChemModelBase {
 
  public:
   CaloricallyPerfectThermoChem(mfem::ParMesh *pmesh, LoMachOptions *loMach_opts, temporalSchemeCoefficients &timeCoeff,
-                               TPS::Tps *tps);
+                               ParGridFunction *gridScale, TPS::Tps *tps);
   virtual ~CaloricallyPerfectThermoChem();
 
   // Functions overriden from base class

src/em_options.hpp

@@ -52,8 +52,8 @@ class ElectromagneticOptions {
   double atol;                            /**< Linear solver absolute tolerance */
   double preconditioner_background_sigma; /**< Uniform conductivity to use to preconditioner (ignored if <= 0) */
 
-  bool top_only; /**< Flag to specify current in top rings only */
-  bool bot_only; /**< Flag to specify current in bottom rings only */
+  bool top_only;         /**< Flag to specify current in top rings only */
+  bool bot_only;         /**< Flag to specify current in bottom rings only */
   bool variable_current; /**< Flag to specify variable current in each ring */
 
   bool evaluate_magnetic_field;

src/loMach.cpp

@@ -166,11 +166,12 @@ void LoMachSolver::initialize() {
   if (loMach_opts_.thermo_solver == "constant-property") {
     thermo_ = new ConstantPropertyThermoChem(pmesh_, loMach_opts_.order, tpsP_);
   } else if (loMach_opts_.thermo_solver == "calorically-perfect") {
-    thermo_ = new CaloricallyPerfectThermoChem(pmesh_, &loMach_opts_, temporal_coeff_, tpsP_);
+    thermo_ =
+        new CaloricallyPerfectThermoChem(pmesh_, &loMach_opts_, temporal_coeff_, (meshData_->getGridScale()), tpsP_);
   } else if (loMach_opts_.thermo_solver == "lte-thermo-chem") {
     thermo_ = new LteThermoChem(pmesh_, &loMach_opts_, temporal_coeff_, tpsP_);
   } else if (loMach_opts_.thermo_solver == "reacting-flow") {
-    thermo_ = new ReactingFlow(pmesh_, &loMach_opts_, temporal_coeff_, tpsP_);
+    thermo_ = new ReactingFlow(pmesh_, &loMach_opts_, temporal_coeff_, (meshData_->getGridScale()), tpsP_);
   } else {
     // Unknown choice... die
     if (rank0_) {
@@ -182,10 +183,11 @@ void LoMachSolver::initialize() {
   // Instantiate flow solver
   if (loMach_opts_.flow_solver == "zero-flow") {
     // No flow---set u = 0.  Primarily useful for testing thermochem models in isolation
-    flow_ = new ZeroFlow(pmesh_, 1);
+    flow_ = new ZeroFlow(pmesh_, 1, tpsP_);
   } else if (loMach_opts_.flow_solver == "tomboulides") {
     // Tomboulides flow solver
-    flow_ = new Tomboulides(pmesh_, loMach_opts_.order, loMach_opts_.order, temporal_coeff_, tpsP_);
+    flow_ = new Tomboulides(pmesh_, loMach_opts_.order, loMach_opts_.order, temporal_coeff_,
+                            (meshData_->getGridScale()), tpsP_);
   } else {
     // Unknown choice... die
     if (rank0_) {

src/lte_thermo_chem.cpp

@@ -453,7 +453,7 @@ void LteThermoChem::initializeSelf() {
 
   // outlet bc
   {
-      // Assumed homogeneous Nuemann on T, so nothing to do
+    // Assumed homogeneous Nuemann on T, so nothing to do
   }
 
   // Wall BCs

src/quasimagnetostatic.cpp

@@ -272,9 +272,7 @@ void QuasiMagnetostaticSolver3D::InitializeCurrent() {
   }
 
   if (rank0_) {
-    std::cout << "J0 = " << J0(0) << ", " << J0(1)
-              << ", " << J0(2) << ", " << J0(3)
-              << ", " << J0(4) << endl;
+    std::cout << "J0 = " << J0(0) << ", " << J0(1) << ", " << J0(2) << ", " << J0(3) << ", " << J0(4) << endl;
   }
 
   PWConstCoefficient J0coef(J0);
@@ -924,9 +922,7 @@ void QuasiMagnetostaticSolverAxiSym::InitializeCurrent() {
   }
 
   if (rank0_) {
-    std::cout << "J0 = " << J0(0) << ", " << J0(1)
-              << ", " << J0(2) << ", " << J0(3)
-              << ", " << J0(4) << endl;
+    std::cout << "J0 = " << J0(0) << ", " << J0(1) << ", " << J0(2) << ", " << J0(3) << ", " << J0(4) << endl;
   }
 
   FunctionCoefficient radius_coeff(radius);