Explicit streamwise stabilization

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)
-    : tpsP_(tps), pmesh_(pmesh), time_coeff_(time_coeff) {
+                                                           temporalSchemeCoefficients &time_coeff,
+                                                           ParGridFunction *gridScale, TPS::Tps *tps)
+    : tpsP_(tps), pmesh_(pmesh), dim_(pmesh->Dimension()), 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,10 @@ 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_);
+
+  tps->getInput("loMach/calperfect/streamwise-stabilization", sw_stab_, false);
+  tps->getInput("loMach/calperfect/Reh_offset", re_offset_, 1.0);
+  tps->getInput("loMach/calperfect/Reh_factor", re_factor_, 0.1);
 }
 
 CaloricallyPerfectThermoChem::~CaloricallyPerfectThermoChem() {
@@ -147,10 +153,15 @@ CaloricallyPerfectThermoChem::~CaloricallyPerfectThermoChem() {
   delete HtInvPC_;
   delete MsInv_;
   delete MsInvPC_;
+  delete Mv_inv_;
+  delete Mv_inv_pc_;
   delete Ht_form_;
   delete MsRho_form_;
   delete Ms_form_;
+  delete Mv_form_;
   delete At_form_;
+  delete G_form_;
+  delete D_form_;
   delete rhou_coeff_;
   delete rhon_next_coeff_;
   delete un_next_coeff_;
@@ -167,6 +178,8 @@ CaloricallyPerfectThermoChem::~CaloricallyPerfectThermoChem() {
   // allocated in initializeSelf
   delete sfes_;
   delete sfec_;
+  delete vfes_;
+  delete vfec_;
 }
 
 void CaloricallyPerfectThermoChem::initializeSelf() {
@@ -178,6 +191,9 @@ void CaloricallyPerfectThermoChem::initializeSelf() {
   sfec_ = new H1_FECollection(order_);
   sfes_ = new ParFiniteElementSpace(pmesh_, sfec_);
 
+  vfec_ = new H1_FECollection(order_, dim_);
+  vfes_ = new ParFiniteElementSpace(pmesh_, vfec_, dim_);
+
   // Check if fully periodic mesh
   if (!(pmesh_->bdr_attributes.Size() == 0)) {
     temp_ess_attr_.SetSize(pmesh_->bdr_attributes.Max());
@@ -189,6 +205,7 @@ void CaloricallyPerfectThermoChem::initializeSelf() {
   if (rank0_) grvy_printf(ginfo, "CaloricallyPerfectThermoChem paces constructed...\n");
 
   int sfes_truevsize = sfes_->GetTrueVSize();
+  int vfes_truevsize = vfes_->GetTrueVSize();
 
   Qt_.SetSize(sfes_truevsize);
   Qt_ = 0.0;
@@ -238,8 +255,19 @@ void CaloricallyPerfectThermoChem::initializeSelf() {
   visc_gf_.SetSpace(sfes_);
   visc_gf_ = 0.0;
 
+  vel_gf_.SetSpace(vfes_);
+  tmpR1_gf_.SetSpace(vfes_);
+  tmpR0_gf_.SetSpace(sfes_);
+
+  swDiff_.SetSize(sfes_truevsize);
   tmpR0_.SetSize(sfes_truevsize);
+  tmpR0a_.SetSize(sfes_truevsize);
   tmpR0b_.SetSize(sfes_truevsize);
+  tmpR0c_.SetSize(sfes_truevsize);
+  tmpR1_.SetSize(vfes_truevsize);
+
+  gradT_.SetSize(vfes_truevsize);
+  gradT_ = 0.0;
 
   R0PM0_gf_.SetSpace(sfes_);
 
@@ -471,6 +499,61 @@ void CaloricallyPerfectThermoChem::initializeOperators() {
   Ms_form_->Assemble();
   Ms_form_->FormSystemMatrix(empty, Ms_);
 
+  // Divergence operator
+  D_form_ = new ParMixedBilinearForm(vfes_, sfes_);
+  VectorDivergenceIntegrator *vd_mblfi;
+  // if (axisym_) {
+  //  vd_mblfi = new VectorDivergenceIntegrator(radius_coeff);
+  // } else {
+  vd_mblfi = new VectorDivergenceIntegrator();
+  // }
+  if (numerical_integ_) {
+    vd_mblfi->SetIntRule(&ir_nli);
+  }
+  D_form_->AddDomainIntegrator(vd_mblfi);
+  if (partial_assembly_) {
+    D_form_->SetAssemblyLevel(AssemblyLevel::PARTIAL);
+  }
+  D_form_->Assemble();
+  D_form_->FormRectangularSystemMatrix(empty, empty, D_op_);
+
+  // Gradient
+  G_form_ = new ParMixedBilinearForm(sfes_, vfes_);
+  // auto *g_mblfi = new GradientIntegrator();
+  GradientIntegrator *g_mblfi;
+  // if (axisym_) {
+  //   g_mblfi = new GradientIntegrator(radius_coeff);
+  // } else {
+  g_mblfi = new GradientIntegrator();
+  //}
+  if (numerical_integ_) {
+    g_mblfi->SetIntRule(&ir_nli);
+  }
+  G_form_->AddDomainIntegrator(g_mblfi);
+  if (partial_assembly_) {
+    G_form_->SetAssemblyLevel(AssemblyLevel::PARTIAL);
+  }
+  G_form_->Assemble();
+  G_form_->FormRectangularSystemMatrix(empty, empty, G_op_);
+
+  // Mass matrix for the vector (gradT)
+  Mv_form_ = new ParBilinearForm(vfes_);
+  VectorMassIntegrator *mv_blfi;
+  // if (axisym_) {
+  //   mv_blfi = new VectorMassIntegrator(radius_coeff);
+  // } else {
+  mv_blfi = new VectorMassIntegrator();
+  //}
+  if (numerical_integ_) {
+    mv_blfi->SetIntRule(&ir_nli);
+  }
+  Mv_form_->AddDomainIntegrator(mv_blfi);
+  if (partial_assembly_) {
+    Mv_form_->SetAssemblyLevel(AssemblyLevel::PARTIAL);
+  }
+  Mv_form_->Assemble();
+  Mv_form_->FormSystemMatrix(temp_ess_tdof_, Mv_);
+
   // mass matrix with rho
   MsRho_form_ = new ParBilinearForm(sfes_);
   auto *msrho_blfi = new MassIntegrator(*rho_coeff_);
@@ -483,6 +566,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_);
@@ -517,6 +601,27 @@ void CaloricallyPerfectThermoChem::initializeOperators() {
   MsInv_->SetAbsTol(mass_inverse_atol_);
   MsInv_->SetMaxIter(mass_inverse_max_iter_);
 
+  // Inverse (unweighted) mass operator (velocity space)
+  if (partial_assembly_) {
+    Vector diag_pa(vfes_->GetTrueVSize());
+    Mv_form_->AssembleDiagonal(diag_pa);
+    Mv_inv_pc_ = new OperatorJacobiSmoother(diag_pa, empty);
+  } else {
+    Mv_inv_pc_ = new HypreSmoother(*Mv_.As<HypreParMatrix>());
+    dynamic_cast<HypreSmoother *>(Mv_inv_pc_)->SetType(HypreSmoother::Jacobi, smoother_passes_);
+    dynamic_cast<HypreSmoother *>(Mv_inv_pc_)->SetSOROptions(smoother_relax_weight_, smoother_relax_omega_);
+    dynamic_cast<HypreSmoother *>(Mv_inv_pc_)
+        ->SetPolyOptions(smoother_poly_order_, smoother_poly_fraction_, smoother_eig_est_);
+  }
+  Mv_inv_ = new CGSolver(vfes_->GetComm());
+  Mv_inv_->iterative_mode = false;
+  Mv_inv_->SetOperator(*Mv_);
+  Mv_inv_->SetPreconditioner(*Mv_inv_pc_);
+  Mv_inv_->SetPrintLevel(mass_inverse_pl_);
+  Mv_inv_->SetAbsTol(mass_inverse_atol_);
+  Mv_inv_->SetRelTol(mass_inverse_rtol_);
+  Mv_inv_->SetMaxIter(mass_inverse_max_iter_);
+
   HtInvPC_ = new HypreSmoother(*Ht_.As<HypreParMatrix>());
   dynamic_cast<HypreSmoother *>(HtInvPC_)->SetType(HypreSmoother::Jacobi, smoother_passes_);
   dynamic_cast<HypreSmoother *>(HtInvPC_)->SetSOROptions(hsmoother_relax_weight_, hsmoother_relax_omega_);
@@ -666,6 +771,17 @@ void CaloricallyPerfectThermoChem::step() {
   tmpR0_ = (dtP_ / Cp_);
   Ms_->AddMult(tmpR0_, resT_);
 
+  // Add streamwise stability to rhs
+  if (sw_stab_) {
+    // compute temp gradient (only really needed for sw-stab)
+    G_op_->Mult(Text_, tmpR1_);
+    Mv_inv_->Mult(tmpR1_, gradT_);
+
+    // streamwiseDiffusion(Tn_, swDiff_);
+    streamwiseDiffusion(gradT_, swDiff_);
+    resT_.Add(1.0, swDiff_);
+  }
+
   // Add natural boundary terms here later
   // NB: adiabatic natural BC is handled, but don't have ability to impose non-zero heat flux yet
 
@@ -1022,3 +1138,56 @@ void CaloricallyPerfectThermoChem::screenValues(std::vector<double> &values) {
     values[0] = thermo_pressure_ / ambient_pressure_;
   }
 }
+
+// void CaloricallyPerfectThermoChem::streamwiseDiffusion(Vector &phi, Vector &swDiff) {
+void CaloricallyPerfectThermoChem::streamwiseDiffusion(Vector &gradPhi, Vector &swDiff) {
+  (flow_interface_->velocity)->GetTrueDofs(tmpR0a_);
+  vel_gf_.SetFromTrueDofs(tmpR0a_);
+
+  // compute streamwise gradient of input field
+  // tmpR0_gf_.SetFromTrueDofs(phi);
+  // streamwiseGrad(dim_, tmpR0_gf_, vel_gf_, tmpR1_gf_);
+
+  tmpR1_gf_.SetFromTrueDofs(gradPhi);
+  streamwiseGrad(dim_, vel_gf_, tmpR1_gf_);
+
+  // divergence of sw-grad
+  tmpR1_gf_.GetTrueDofs(tmpR1_);
+  D_op_->Mult(tmpR1_, swDiff);
+
+  gridScale_gf_->GetTrueDofs(tmpR0b_);
+  // (turbModel_interface_->eddy_viscosity)->GetTrueDofs(tmpR0c_);
+  (flow_interface_->Reh)->GetTrueDofs(tmpR0c_);
+
+  upwindDiff(dim_, re_factor_, re_offset_, tmpR0a_, rn_, tmpR0b_, tmpR0c_, swDiff);
+
+  /*
+  const double *rho = rn_.HostRead();
+  const double *vel = tmpR0a_.HostRead();
+  const double *del = tmpR0b_.HostRead();
+  //const double *mu = visc_.HostRead();
+  //const double *muT = tmpR0c_.HostRead();
+  const double *Reh = tmpR0c_.HostRead();
+  double *data = swDiff.HostReadWrite();
+
+  int Sdof = rn_.Size();
+  for (int dof = 0; dof < Sdof; dof++) {
+    double Umag = 0.0;
+    for (int i = 0; i < dim_; i++) Umag += vel[i] * vel[i];
+    Umag = std::sqrt(Umag);
+
+    // element Re
+    //double Re = Umag * del[dof] * rho[dof] / (mu[dof] + muT[dof]);
+    double Re = Reh[dof];
+
+    // SUPG weight
+    double Csupg = 0.5 * (tanh(re_factor_ * Re - re_offset_) + 1.0);
+
+    // streamwise diffusion coeff
+    double CswDiff = Csupg * Umag * del[dof] * rho[dof];
+
+    // scaled streamwise Laplacian
+    data[dof] *= CswDiff;
+  }
+  */
+}

src/calorically_perfect.hpp

@@ -65,6 +65,7 @@ class CaloricallyPerfectThermoChem : public ThermoChemModelBase {
 
   // Mesh and discretization scheme info
   ParMesh *pmesh_ = nullptr;
+  int dim_;
   int order_;
   IntegrationRules gll_rules_;
   const temporalSchemeCoefficients &time_coeff_;
@@ -110,6 +111,11 @@ class CaloricallyPerfectThermoChem : public ThermoChemModelBase {
   double hsolve_rtol_;
   double hsolve_atol_;
 
+  // streamwise-stabilization
+  bool sw_stab_;
+  double re_offset_;
+  double re_factor_;
+
   // Boundary condition info
   Array<int> temp_ess_attr_; /**< List of patches with Dirichlet BC on temperature */
   Array<int> Qt_ess_attr_;   /**< List of patches with Dirichlet BC on Q (thermal divergence) */
@@ -144,6 +150,10 @@ class CaloricallyPerfectThermoChem : public ThermoChemModelBase {
   // Scalar \f$H^1\f$ finite element space.
   ParFiniteElementSpace *sfes_ = nullptr;
 
+  // Vector fe collection and space
+  FiniteElementCollection *vfec_ = nullptr;
+  ParFiniteElementSpace *vfes_ = nullptr;
+
   // Fields
   ParGridFunction Tnm1_gf_, Tnm2_gf_;
   ParGridFunction Tn_gf_, Tn_next_gf_, Text_gf_, resT_gf_;
@@ -155,6 +165,11 @@ class CaloricallyPerfectThermoChem : public ThermoChemModelBase {
   ParGridFunction R0PM0_gf_;
   ParGridFunction Qt_gf_;
 
+  ParGridFunction tmpR0_gf_;
+  ParGridFunction tmpR1_gf_;
+  ParGridFunction vel_gf_;
+  ParGridFunction *gridScale_gf_ = nullptr;
+
   // ParGridFunction *buffer_tInlet_ = nullptr;
   GridFunctionCoefficient *temperature_bc_field_ = nullptr;
 
@@ -174,10 +189,13 @@ class CaloricallyPerfectThermoChem : public ThermoChemModelBase {
   // operators and solvers
   ParBilinearForm *At_form_ = nullptr;
   ParBilinearForm *Ms_form_ = nullptr;
+  ParBilinearForm *Mv_form_ = nullptr;
   ParBilinearForm *MsRho_form_ = nullptr;
   ParBilinearForm *Ht_form_ = nullptr;
   ParBilinearForm *Mq_form_ = nullptr;
   ParBilinearForm *LQ_form_ = nullptr;
+  ParMixedBilinearForm *D_form_ = nullptr;
+  ParMixedBilinearForm *G_form_ = nullptr;
   ParLinearForm *LQ_bdry_ = nullptr;
 
   OperatorHandle LQ_;
@@ -186,20 +204,28 @@ class CaloricallyPerfectThermoChem : public ThermoChemModelBase {
   OperatorHandle Ms_;
   OperatorHandle MsRho_;
   OperatorHandle Mq_;
+  OperatorHandle Mv_;
+  OperatorHandle D_op_;
+  OperatorHandle G_op_;
 
   mfem::Solver *MsInvPC_ = nullptr;
   mfem::CGSolver *MsInv_ = nullptr;
   mfem::Solver *MqInvPC_ = nullptr;
   mfem::CGSolver *MqInv_ = nullptr;
   mfem::Solver *HtInvPC_ = nullptr;
   mfem::CGSolver *HtInv_ = nullptr;
+  mfem::Solver *Mv_inv_pc_ = nullptr;
+  mfem::CGSolver *Mv_inv_ = nullptr;
 
   // Vectors
   Vector Tn_, Tn_next_, Tnm1_, Tnm2_;
   Vector NTn_, NTnm1_, NTnm2_;
   Vector Text_;
   Vector resT_;
-  Vector tmpR0_, tmpR0b_;
+  Vector tmpR0_, tmpR0a_, tmpR0b_, tmpR0c_;
+  Vector tmpR1_;
+  Vector swDiff_;
+  Vector gradT_;
 
   Vector Qt_;
   Vector rn_;
@@ -233,7 +259,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
@@ -257,6 +283,8 @@ class CaloricallyPerfectThermoChem : public ThermoChemModelBase {
   void computeExplicitTempConvectionOP(bool extrap);
   void computeQt();
   void computeQtTO();
+  //  void streamwiseDiffusion(Vector &phi, Vector &swDiff);
+  void streamwiseDiffusion(Vector &gradPhi, Vector &swDiff);
 
   /// Return a pointer to the current temperature ParGridFunction.
   ParGridFunction *GetCurrentTemperature() { return &Tn_gf_; }

src/cases.cpp

@@ -204,7 +204,7 @@ double temp_channel(const Vector &coords, double t) {
   double Tlo = 200.0;
   double y = coords(1);
   double temp;
-  temp = Tlo + (y + 0.5) * (Thi - Tlo);
+  temp = Tlo + 0.5 * (y + 0.1) * (Thi - Tlo);
   return temp;
 }
 

src/loMach.cpp

@@ -166,9 +166,10 @@ 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_);
+    thermo_ = new LteThermoChem(pmesh_, &loMach_opts_, temporal_coeff_, (meshData_->getGridScale()), tpsP_);
   } else if (loMach_opts_.thermo_solver == "reacting-flow") {
     thermo_ = new ReactingFlow(pmesh_, &loMach_opts_, temporal_coeff_, tpsP_);
   } else {
@@ -185,7 +186,8 @@ void LoMachSolver::initialize() {
     flow_ = new ZeroFlow(pmesh_, 1);
   } 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_) {