Low-Mach averaging

src/algebraicSubgridModels.cpp

@@ -114,7 +114,7 @@ void AlgebraicSubgridModels::initializeSelf() {
   // gridScaleZ_.SetSize(sfes_truevsize);
   // resolution_gf_.SetSpace(sfes_);
 
-  if (verbose) grvy_printf(ginfo, "AlgebraicSubgridModels vectors and gf initialized...\n");
+  if (rank0_) grvy_printf(ginfo, "AlgebraicSubgridModels vectors and gf initialized...\n");
 
   // exports
   toFlow_interface_.eddy_viscosity = &subgridVisc_gf_;

src/averaging.cpp

@@ -75,6 +75,9 @@ Averaging::Averaging(AveragingOptions &opts, std::string output_name) {
 
   mean_output_name_ = "mean_";
   mean_output_name_.append(output_name);
+
+  enable_mean_continuation_ = opts.enable_mean_continuation_;
+  zero_variances_ = opts.zero_variances_;
 }
 
 Averaging::~Averaging() {
@@ -119,6 +122,7 @@ void Averaging::registerField(std::string name, const ParGridFunction *field_to_
 
   // and store those fields in an AveragingFamily object that gets appended to the avg_families_ vector
   avg_families_.emplace_back(AveragingFamily(name, field_to_average, mean, vari, vari_start_index, vari_components));
+  std::cout << "Fam size: " << avg_families_.size() << " with addition of " << name << endl;
 }
 
 void Averaging::initializeViz() {
@@ -198,15 +202,37 @@ void Averaging::addSample(const int &iter, GasMixture *mixture) {
 
   assert(avg_families_.size() >= 1);
 
+  /*
+  if (!enable_mean_continuation_) {
+    ns_mean_ = 0;
+    enable_mean_continuation_ = true;
+  }
+
+  if (zero_variances_) {
+    ns_vari_ = 0;
+    zero_variances_ = false;
+  }
+  */
+
   if (iter % sample_interval_ == 0 && iter >= step_start_mean_) {
-    if (iter == step_start_mean_ && rank0_) cout << "Starting mean calculation." << endl;
+    if (iter == step_start_mean_ && rank0_) cout << "Starting mean calculation at iter " << iter << endl;
+
+    if (ns_mean_ == 0) {
+      for (size_t ifam = 0; ifam < avg_families_.size(); ifam++) {
+        if (avg_families_[ifam].mean_fcn_ != nullptr) {
+          *avg_families_[ifam].mean_fcn_ = 0.0;
+        }
+      }
+    }
 
+    // If we got here, then either this is our first time averaging
+    // or the variances have been restarted.  Either way, valid to
+    // set variances to zero.
     if (ns_vari_ == 0) {
-      // If we got here, then either this is our first time averaging
-      // or the variances have been restarted.  Either way, valid to
-      // set variances to zero.
       for (size_t ifam = 0; ifam < avg_families_.size(); ifam++) {
-        *avg_families_[ifam].vari_fcn_ = 0.0;
+        if (avg_families_[ifam].vari_fcn_ != nullptr) {
+          *avg_families_[ifam].vari_fcn_ = 0.0;
+        }
       }
     }
 
@@ -251,18 +277,28 @@ void Averaging::addSampleInternal() {
 
     const double *d_inst = inst->Read();
     double *d_mean = mean->ReadWrite();
-    double *d_vari = vari->ReadWrite();
+    double *d_vari = nullptr;
+    if (vari != nullptr) {
+      d_vari = vari->ReadWrite();
+    }
 
     // Extract size information for use on device
     const int dof = mean->ParFESpace()->GetNDofs();  // dofs per scalar field
     const int neq = mean->ParFESpace()->GetVDim();   // number of scalar variables in mean field
 
-    const int d_vari_start = fam.vari_start_index_;
-    const int d_vari_components = fam.vari_components_;
+    int d_vari_start = 0;
+    int d_vari_components = 0;
+    if (vari != nullptr) {
+      d_vari_start = fam.vari_start_index_;
+      d_vari_components = fam.vari_components_;
+    }
 
     // Extract sample size information for use on device
     double d_ns_mean = (double)ns_mean_;
-    double d_ns_vari = (double)ns_vari_;
+    double d_ns_vari = 0;
+    if (vari != nullptr) {
+      d_ns_vari = (double)ns_vari_;
+    }
 
     // "Loop" over all dofs and update statistics
     MFEM_FORALL(n, dof, {
@@ -296,7 +332,8 @@ void Averaging::addSampleInternal() {
           const double uinst_i = d_inst[n + i * dof];
           const double umean_i = d_mean[n + i * dof];
           delta_i = uinst_i - umean_i;
-          for (int j = d_vari_start + 1; j < d_vari_start + d_vari_components; j++) {
+
+          for (int j = i + 1; j < d_vari_start + d_vari_components; j++) {
             const double uinst_j = d_inst[n + j * dof];
             const double umean_j = d_mean[n + j * dof];
             delta_j = uinst_j - umean_j;

src/averaging.hpp

@@ -177,6 +177,12 @@ class Averaging {
   /// Visualization directory (i.e., paraview dumped to mean_output_name_)
   std::string mean_output_name_;
 
+  /// flag to continue mean from restart
+  bool enable_mean_continuation_;
+
+  /// flag to restart variances
+  bool zero_variances_;
+
   /// mfem paraview data collection, used to write viz files
   ParaViewDataCollection *pvdc_ = nullptr;
 
@@ -282,7 +288,10 @@ class Averaging {
   int GetSamplesMean() { return ns_mean_; }
   int GetSamplesRMS() { return ns_vari_; }
   int GetSamplesInterval() { return sample_interval_; }
+
   bool ComputeMean() { return compute_mean_; }
+  bool ContinueMean() { return enable_mean_continuation_; }
+  bool RestartRMS() { return zero_variances_; }
 
   void SetSamplesMean(int &samples) { ns_mean_ = samples; }
   void SetSamplesRMS(int &samples) { ns_vari_ = samples; }

src/calorically_perfect.cpp

@@ -47,6 +47,10 @@
 using namespace mfem;
 using namespace mfem::common;
 
+/// forward declarations
+double temp_rt3d(const Vector &x, double t);
+double temp_channel(const Vector &x, double t);
+
 MFEM_HOST_DEVICE double Sutherland(const double T, const double mu_star, const double T_star, const double S_star) {
   const double T_rat = T / T_star;
   const double T_rat_32 = T_rat * sqrt(T_rat);
@@ -239,6 +243,9 @@ void CaloricallyPerfectThermoChem::initializeSelf() {
     std::cout << "exports set..." << endl;
   }
 
+  tpsP_->getInput("loMach/calperfect/numerical-integ", numerical_integ_, true);
+  tpsP_->getInput("loMach/calperfect/over-integ", over_integrate_, false);
+
   //-----------------------------------------------------
   // 2) Set the initial condition
   //-----------------------------------------------------
@@ -256,9 +263,26 @@ void CaloricallyPerfectThermoChem::initializeSelf() {
   // 2) For restarts, this IC is overwritten by the restart field,
   // which is read later.
 
-  ConstantCoefficient t_ic_coef;
-  t_ic_coef.constant = T_ic_;
-  Tn_gf_.ProjectCoefficient(t_ic_coef);
+  tpsP_->getInput("loMach/calperfect/ic", ic_string_, std::string(""));
+
+  // set IC if we have one at this point
+  if (!ic_string_.empty()) {
+    if (ic_string_ == "rt3D") {
+      if (rank0_) std::cout << "Setting rt3D IC..." << std::endl;
+      FunctionCoefficient t_excoeff(temp_rt3d);
+      t_excoeff.SetTime(0.0);
+      Tn_gf_.ProjectCoefficient(t_excoeff);
+    } else if (ic_string_ == "channel") {
+      if (rank0_) std::cout << "Setting channel IC..." << std::endl;
+      FunctionCoefficient t_excoeff(temp_channel);
+      t_excoeff.SetTime(0.0);
+      Tn_gf_.ProjectCoefficient(t_excoeff);
+    }
+  } else {
+    ConstantCoefficient t_ic_coef;
+    t_ic_coef.constant = T_ic_;
+    Tn_gf_.ProjectCoefficient(t_ic_coef);
+  }
 
   Tn_gf_.GetTrueDofs(Tn_);
   Tnm1_gf_.SetFromTrueDofs(Tn_);
@@ -334,16 +358,13 @@ void CaloricallyPerfectThermoChem::initializeSelf() {
     Array<int> attr_outlet(pmesh_->bdr_attributes.Max());
     attr_outlet = 0;
 
-    // No code for this yet, so die if detected
-    // assert(numOutlets == 0);
-
-    // But... outlets will just get homogeneous Neumann on T, so
+    // Outlets will just get homogeneous Neumann on T, so
     // basically need to do nothing.
   }
 
   // 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;
 
@@ -356,7 +377,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;
@@ -391,9 +412,17 @@ void CaloricallyPerfectThermoChem::initializeOperators() {
   // unsteady: p+p [+p] = 2p [3p]
   // convection: p+p+(p-1) [+p] = 3p-1 [4p-1]
   // diffusion: (p-1)+(p-1) [+p] = 2p-2 [3p-2]
-  const IntegrationRule &ir_i = gll_rules_.Get(sfes_->GetFE(0)->GetGeomType(), 2 * order_ + 1);
-  const IntegrationRule &ir_nli = gll_rules_.Get(sfes_->GetFE(0)->GetGeomType(), 4 * order_);
-  const IntegrationRule &ir_di = gll_rules_.Get(sfes_->GetFE(0)->GetGeomType(), 3 * order_ - 1);
+  int nir_i = 2 * order_ - 1;
+  int nir_nli = 2 * order_ - 1;
+  int nir_di = 2 * order_ - 1;
+  if (over_integrate_) {
+    nir_i = 3 * order_ + 1;
+    nir_nli = 4 * order_;
+    nir_di = 3 * order_ - 1;
+  }
+  const IntegrationRule &ir_i = gll_rules_.Get(sfes_->GetFE(0)->GetGeomType(), nir_i);
+  const IntegrationRule &ir_nli = gll_rules_.Get(sfes_->GetFE(0)->GetGeomType(), nir_nli);
+  const IntegrationRule &ir_di = gll_rules_.Get(sfes_->GetFE(0)->GetGeomType(), nir_di);
   if (rank0_) std::cout << "Integration rules set" << endl;
 
   // coefficients for operators
@@ -733,6 +762,18 @@ void CaloricallyPerfectThermoChem::initializeViz(ParaViewDataCollection &pvdc) {
   pvdc.RegisterField("Qt", &Qt_gf_);
 }
 
+void CaloricallyPerfectThermoChem::initializeStats(Averaging &average, IODataOrganizer &io, bool continuation) {
+  if (average.ComputeMean()) {
+    // fields for averaging
+    average.registerField(std::string("temperature"), &Tn_gf_, false, 0, 1);
+
+    // io init
+    io.registerIOFamily("Time-averaged temperature", "/meanTemp", average.GetMeanField(std::string("temperature")),
+                        false, continuation, sfec_);
+    io.registerIOVar("/meanTemp", "<T>", 0, true);
+  }
+}
+
 /**
    Update boundary conditions for Temperature, useful for
    ramping a dirichlet condition over time
@@ -1426,3 +1467,37 @@ double temp_inlet(const Vector &coords, double t) {
   return temp;
 }
 #endif
+
+double temp_rt3d(const Vector &x, double t) {
+  double CC = 0.05;
+  double twoPi = 6.28318530718;
+  double yWidth = 0.1;
+  double yInt, dy, wt;
+  double temp, dT;
+  double Tlo = 100.0;
+  double Thi = 1500.0;
+
+  yInt = std::cos(twoPi * x[0]) + std::cos(twoPi * x[2]);
+  yInt *= CC;
+  yInt += 4.0;
+
+  dy = x[1] - yInt;
+  dT = Thi - Tlo;
+
+  wt = 0.5 * (tanh(-dy / yWidth) + 1.0);
+  temp = Tlo + wt * dT;
+
+  return temp;
+}
+
+/// Used to set the channel IC
+double temp_channel(const Vector &x, double t) {
+  double Thi = 400.0;
+  double Tlo = 200.0;
+  double temp = 300.0;
+
+  // expects channel height (-1,1)
+  temp = Tlo + (Thi - Tlo) * 0.5 * (x(1) + 1.0);
+
+  return temp;
+}

src/calorically_perfect.hpp

@@ -43,6 +43,7 @@ class Tps;
 
 #include <iostream>
 
+#include "averaging.hpp"
 #include "dirichlet_bc_helper.hpp"
 #include "io.hpp"
 #include "thermo_chem_base.hpp"
@@ -70,10 +71,13 @@ class CaloricallyPerfectThermoChem : public ThermoChemModelBase {
   double dt_;
   double time_;
 
+  std::string ic_string_;
+
   // Flags
   bool rank0_;                      /**< true if this is rank 0 */
   bool partial_assembly_ = false;   /**< Enable/disable partial assembly of forms. */
   bool numerical_integ_ = true;     /**< Enable/disable numerical integration rules of forms. */
+  bool over_integrate_ = false;     /**< Enable/disable numerical integration with higher order. */
   bool constant_viscosity_ = false; /**< Enable/disable constant viscosity */
   bool constant_density_ = false;   /**< Enable/disable constant density */
   bool domain_is_open_ = false;     /**< true if domain is open */
@@ -215,6 +219,7 @@ class CaloricallyPerfectThermoChem : public ThermoChemModelBase {
   void step() final;
   void initializeIO(IODataOrganizer &io) final;
   void initializeViz(ParaViewDataCollection &pvdc) final;
+  void initializeStats(Averaging &average, IODataOrganizer &io, bool continuation) final;
 
   void screenHeader(std::vector<std::string> &header) const final;
   void screenValues(std::vector<double> &values) final;