From f058d3305d5c63a5f2f655fb3e8400945f69b627 Mon Sep 17 00:00:00 2001
From: Sigfried Haering <shaering@oden.utexas.edu>
Date: Wed, 26 Jun 2024 14:59:35 -0700
Subject: [PATCH] adding 4-step time-averaging in addition to p-filter for SGS
 mu_T

---
 src/algebraicSubgridModels.cpp | 38 +++++++++++++++++++++++
 src/algebraicSubgridModels.hpp |  6 ++++
 src/utils.cpp                  | 57 ++++++++++++++++++++++++++++++++++
 src/utils.hpp                  |  1 +
 4 files changed, 102 insertions(+)

diff --git a/src/algebraicSubgridModels.cpp b/src/algebraicSubgridModels.cpp
index 028fe1e76..b02d5b482 100644
--- a/src/algebraicSubgridModels.cpp
+++ b/src/algebraicSubgridModels.cpp
@@ -75,6 +75,8 @@ AlgebraicSubgridModels::AlgebraicSubgridModels(mfem::ParMesh *pmesh, LoMachOptio
     }
     sgs_model_nFilter_ = 0;
   }
+
+  tpsP_->getInput("loMach/sgsSmooth", sgs_model_smooth_, false);
 }
 
 AlgebraicSubgridModels::~AlgebraicSubgridModels() {
@@ -114,6 +116,16 @@ void AlgebraicSubgridModels::initializeSelf() {
   subgridVisc_gf_.SetSpace(sfes_);
   delta_.SetSize(sfes_truevsize);
 
+  muT_NM0_.SetSize(sfes_truevsize);
+  muT_NM1_.SetSize(sfes_truevsize);
+  muT_NM2_.SetSize(sfes_truevsize);
+  muT_NM3_.SetSize(sfes_truevsize);
+  muT_NM0_ = 0.0;
+  muT_NM1_ = 0.0;
+  muT_NM2_ = 0.0;
+  muT_NM3_ = 0.0;
+  aveSteps_ = 0;
+
   gradU_.SetSize(vfes_truevsize);
   gradV_.SetSize(vfes_truevsize);
   gradW_.SetSize(vfes_truevsize);
@@ -221,6 +233,32 @@ void AlgebraicSubgridModels::step() {
                 { d_muT_gf[i] = (1.0 - filter_alpha) * d_muT_gf[i] + filter_alpha * d_muT_filtered_gf[i]; });
     subgridVisc_gf_.GetTrueDofs(subgridVisc_);
   }
+
+  // clip any small negatives resulting from filtering
+  double *dmuT = subgridVisc_.HostReadWrite();
+  for (int i = 0; i < SdofInt_; i++) {
+    dmuT[i] = std::max(dmuT[i], 1.0e-12);
+  }
+
+  if (sgs_model_smooth_) {
+    aveSteps_++;
+    aveSteps_ = std::min(aveSteps_, 4);
+
+    // take average of recent steps
+    muT_NM0_ = subgridVisc_;
+    double Cave = 1.0 / (double)aveSteps_;
+    subgridVisc_ = 0.0;
+    subgridVisc_.Add(Cave, muT_NM0_);
+    subgridVisc_.Add(Cave, muT_NM1_);
+    subgridVisc_.Add(Cave, muT_NM2_);
+    subgridVisc_.Add(Cave, muT_NM3_);
+
+    // shift storage
+    muT_NM3_ = muT_NM2_;
+    muT_NM2_ = muT_NM1_;
+    muT_NM1_ = subgridVisc_;
+  }
+  subgridVisc_gf_.SetFromTrueDofs(subgridVisc_);
 }
 
 /**
diff --git a/src/algebraicSubgridModels.hpp b/src/algebraicSubgridModels.hpp
index 86f7653b2..32d1b2d82 100644
--- a/src/algebraicSubgridModels.hpp
+++ b/src/algebraicSubgridModels.hpp
@@ -130,6 +130,10 @@ class AlgebraicSubgridModels : public TurbModelBase {
 
   ParGridFunction subgridVisc_gf_;
   Vector subgridVisc_;
+  Vector muT_NM0_;
+  Vector muT_NM1_;
+  Vector muT_NM2_;
+  Vector muT_NM3_;
 
   // filter related
   ParGridFunction muT_NM1_gf_;
@@ -140,6 +144,8 @@ class AlgebraicSubgridModels : public TurbModelBase {
 
   double sgs_model_const_;
   int sgs_model_nFilter_;
+  bool sgs_model_smooth_;
+  int aveSteps_;
 
  public:
   AlgebraicSubgridModels(mfem::ParMesh *pmesh, LoMachOptions *loMach_opts, TPS::Tps *tps, ParGridFunction *gridScale,
diff --git a/src/utils.cpp b/src/utils.cpp
index 6af93458d..06477fadd 100644
--- a/src/utils.cpp
+++ b/src/utils.cpp
@@ -1156,6 +1156,63 @@ void Orthogonalize(Vector &v, MPI_Comm comm) {
   v -= global_sum / static_cast<double>(global_size);
 }
 
+void makeContinuous(ParGridFunction &u) {
+  FiniteElementSpace *fes = u.FESpace();
+
+  ParGridFunction au;
+  au.SetSpace(fes);
+  au = 0.0;
+
+  int elndofs;
+  Array<int> vdofs;
+  Vector vals;
+  Vector loc_data;
+  int vdim = fes->GetVDim();
+  Array<int> zones_per_vdof;
+  zones_per_vdof.SetSize(fes->GetVSize());
+  zones_per_vdof = 0;
+
+  for (int e = 0; e < fes->GetNE(); ++e) {
+    fes->GetElementVDofs(e, vdofs);
+    u.GetSubVector(vdofs, loc_data);
+    vals.SetSize(vdofs.Size());
+    const FiniteElement *el = fes->GetFE(e);
+    elndofs = el->GetDof();
+    for (int dof = 0; dof < elndofs; ++dof) {
+      vals(dof) = loc_data(dof);
+    }
+
+    // Accumulate values in all dofs, count the zones.
+    for (int j = 0; j < vdofs.Size(); j++) {
+      int ldof = vdofs[j];
+      au(ldof) += vals[j];
+      zones_per_vdof[ldof]++;
+    }
+  }
+
+  // Communication
+
+  // Count the zones globally.
+  GroupCommunicator &gcomm = u.ParFESpace()->GroupComm();
+  gcomm.Reduce<int>(zones_per_vdof, GroupCommunicator::Sum);
+  gcomm.Bcast(zones_per_vdof);
+
+  // Accumulate for all vdofs.
+  gcomm.Reduce<double>(au.GetData(), GroupCommunicator::Sum);
+  gcomm.Bcast<double>(au.GetData());
+
+  // Compute means.
+  for (int i = 0; i < au.Size(); i++) {
+    const int nz = zones_per_vdof[i];
+    if (nz) {
+      au(i) /= nz;
+    }
+  }
+
+  // copy back
+  u = au;
+}
+
 namespace mfem {
 GradientVectorGridFunctionCoefficient::GradientVectorGridFunctionCoefficient(const GridFunction *gf)
     : MatrixCoefficient((gf) ? gf->VectorDim() : 0) {
diff --git a/src/utils.hpp b/src/utils.hpp
index 4fb5ef324..c5c6dfd4a 100644
--- a/src/utils.hpp
+++ b/src/utils.hpp
@@ -173,6 +173,7 @@ void vectorGrad3D(ParGridFunction &u, ParGridFunction &gu, ParGridFunction &gv,
 void scalarGrad3D(ParGridFunction &u, ParGridFunction &gu);
 void vectorGrad3DV(FiniteElementSpace *fes, Vector u, Vector *gu, Vector *gv, Vector *gw);
 void scalarGrad3DV(FiniteElementSpace *fes, FiniteElementSpace *vfes, Vector u, Vector *gu);
+void makeContinuous(ParGridFunction &u);
 
 /// Eliminate essential BCs in an Operator and apply to RHS.
 /// rename this to something sensible "ApplyEssentialBC" or something