Merge pull request #80 from pecos/sponge_zone

Sponge zone added

.gitattributes

@@ -1,3 +1,5 @@
 test/ref_solns/cyl3d.mflow.2iters.bulkVisc.dtconst.h5 filter=lfs diff=lfs merge=lfs -text
 test/ref_solns/cyl.r0.p2.iter4.h5 filter=lfs diff=lfs merge=lfs -text
 test/ref_solns/cyl.r1.p2.iter4.h5 filter=lfs diff=lfs merge=lfs -text
+test/ref_solns/forcing-periodic-square.h5 filter=lfs diff=lfs merge=lfs -text
+test/meshes/periodic-square.mesh filter=lfs diff=lfs merge=lfs -text

src/dataStructures.hpp

@@ -42,13 +42,44 @@ using namespace mfem;
 // normal equal for both planes. The visc ratio will vary linearly
 // from 0 to viscRatio from the plane defined by pointInit and
 // the plane defined by point0.
+//  |-----------------------------------|
+//  |------------------------------<====|
+//  |-----------------------------------|
+//  p_init                     normal  p0
 struct linearlyVaryingVisc {
   Vector normal;
   Vector point0;
   Vector pointInit;
   double viscRatio;
 };
 
+enum SpongeZoneSolution {
+  USERDEF,   // User defined target solution
+  MIXEDOUT,  // Mixed-out values will be computed at the pInit plane and
+             // used as target solution in the sponge zone
+  NONE = -1  // Invalid value
+};
+
+struct SpongeZoneData {
+  Vector normal;  // These planes are defined in the same manner as
+  Vector point0;  // in the linearlyVaryingVisc struct case.
+  Vector pointInit;
+
+  SpongeZoneSolution szType;
+
+  double tol;  // Tolerance for finding nodes at pInit plane
+               // These points are used to compute mixed-out values.
+               // Required if szType==MIXEDOUT
+
+  Vector targetUp;  // Target user-defined solution defined with primitive
+                    // variables: rho, V, p. Required if szType==USERDEF
+
+  double multFactor;  // Factor multiplying the product of sigma*(U-U*) where
+                      // U* is the target solution. Currently sigam=a/l where
+                      // a is the mixed-out speed of sound and l is the lenght
+                      // of the spnge zone. [OPTIONAL]
+};
+
 struct volumeFaceIntegrationArrays {
   // nodes IDs and indirection array
   Array<int> nodesIDs;

src/em_options.hpp

@@ -59,54 +59,48 @@ class ElectromagneticOptions {
   double yinterp_max;  /**< End value for uniformly spaced points */
 
   ElectromagneticOptions()
-    :
-    mesh_file("hello.msh"), order(1), ref_levels(0),
-    max_iter(100), rtol(1e-6), atol(1e-10),
-    top_only(false), bot_only(false),
-    By_file("By.h5"), nBy(0),
-    yinterp_min(0.0), yinterp_max(1.0)
-  { }
+      : mesh_file("hello.msh"),
+        order(1),
+        ref_levels(0),
+        max_iter(100),
+        rtol(1e-6),
+        atol(1e-10),
+        top_only(false),
+        bot_only(false),
+        By_file("By.h5"),
+        nBy(0),
+        yinterp_min(0.0),
+        yinterp_max(1.0) {}
 
   void AddElectromagneticOptions(mfem::OptionsParser &args) {
     args.AddOption(&mesh_file, "-m", "--mesh", "Mesh file (for EM-only simulation)");
-    args.AddOption(&order, "-o", "--order",
-                   "Finite element order (polynomial degree) (for EM-only).");
-    args.AddOption(&ref_levels, "-r", "--ref",
-                   "Number of uniform refinements (for EM-only).");
-    args.AddOption(&max_iter, "-i", "--maxiter",
-                   "Maximum number of iterations (for EM-only).");
-    args.AddOption(&rtol, "-t", "--rtol",
-                   "Solver relative tolerance (for EM-only).");
-    args.AddOption(&atol, "-a", "--atol",
-                   "Solver absolute tolerance (for EM-only).");
-    args.AddOption(&top_only, "-top", "--top-only", "-ntop", "--no-top-only",
-                   "Run current through top branch only");
-    args.AddOption(&bot_only, "-bot", "--bot-only", "-nbot", "--no-bot-only",
-                   "Run current through bottom branch only");
-    args.AddOption(&By_file, "-by", "--byfile",
-                   "File for By interpolant output.");
-    args.AddOption(&nBy, "-ny", "--nyinterp",
-                   "Number of interpolation points.");
-    args.AddOption(&yinterp_min, "-y0", "--yinterpMin",
-                   "Minimum y interpolation value");
-    args.AddOption(&yinterp_max, "-y1", "--yinterpMax",
-                   "Maximum y interpolation value");
+    args.AddOption(&order, "-o", "--order", "Finite element order (polynomial degree) (for EM-only).");
+    args.AddOption(&ref_levels, "-r", "--ref", "Number of uniform refinements (for EM-only).");
+    args.AddOption(&max_iter, "-i", "--maxiter", "Maximum number of iterations (for EM-only).");
+    args.AddOption(&rtol, "-t", "--rtol", "Solver relative tolerance (for EM-only).");
+    args.AddOption(&atol, "-a", "--atol", "Solver absolute tolerance (for EM-only).");
+    args.AddOption(&top_only, "-top", "--top-only", "-ntop", "--no-top-only", "Run current through top branch only");
+    args.AddOption(&bot_only, "-bot", "--bot-only", "-nbot", "--no-bot-only", "Run current through bottom branch only");
+    args.AddOption(&By_file, "-by", "--byfile", "File for By interpolant output.");
+    args.AddOption(&nBy, "-ny", "--nyinterp", "Number of interpolation points.");
+    args.AddOption(&yinterp_min, "-y0", "--yinterpMin", "Minimum y interpolation value");
+    args.AddOption(&yinterp_max, "-y1", "--yinterpMax", "Maximum y interpolation value");
   }
 
   void print(std::ostream &out) {
     out << std::endl;
     out << "Electromagnetics options:" << std::endl;
-    out << "  mesh_file   = " << mesh_file   << std::endl;
-    out << "  order       = " << order       << std::endl;
-    out << "  ref_levels  = " << ref_levels  << std::endl;
-    out << "  ref_levels  = " << ref_levels  << std::endl;
-    out << "  max_iter    = " << max_iter    << std::endl;
-    out << "  rtol        = " << rtol        << std::endl;
-    out << "  atol        = " << atol        << std::endl;
-    out << "  top_only    = " << top_only    << std::endl;
-    out << "  bot_only    = " << top_only    << std::endl;
-    out << "  By_file     = " << By_file     << std::endl;
-    out << "  nBy         = " << nBy         << std::endl;
+    out << "  mesh_file   = " << mesh_file << std::endl;
+    out << "  order       = " << order << std::endl;
+    out << "  ref_levels  = " << ref_levels << std::endl;
+    out << "  ref_levels  = " << ref_levels << std::endl;
+    out << "  max_iter    = " << max_iter << std::endl;
+    out << "  rtol        = " << rtol << std::endl;
+    out << "  atol        = " << atol << std::endl;
+    out << "  top_only    = " << top_only << std::endl;
+    out << "  bot_only    = " << top_only << std::endl;
+    out << "  By_file     = " << By_file << std::endl;
+    out << "  nBy         = " << nBy << std::endl;
     out << "  yinterp_min = " << yinterp_min << std::endl;
     out << "  yinterp_max = " << yinterp_max << std::endl;
     out << std::endl;

src/equation_of_state.cpp

@@ -109,6 +109,28 @@ double EquationOfState::ComputeMaxCharSpeed(const Vector& state, const int dim)
   return vel + sound;
 }
 
+void EquationOfState::GetConservativesFromPrimitives(const Vector& primit, Vector& conserv, const int& dim,
+                                                     const int& num_equations) {
+  conserv = primit;
+
+  double v2 = 0.;
+  for (int d = 0; d < dim; d++) {
+    v2 += primit[1 + d];
+    conserv[1 + d] *= primit[0];
+  }
+  conserv[num_equations - 1] = primit[num_equations - 1] / (specific_heat_ratio - 1.) + 0.5 * primit[0] * v2;
+}
+
+void EquationOfState::GetPrimitivesFromConservatives(const Vector& conserv, Vector& primit, const int& dim,
+                                                     const int& num_equations) {
+  double p = ComputePressure(conserv, dim);
+  primit = conserv;
+
+  for (int d = 0; d < dim; d++) primit[1 + d] /= conserv[0];
+
+  primit[num_equations - 1] = p;
+}
+
 // GPU FUNCTIONS
 /*#ifdef _GPU_
 double EquationOfState::pressure( double *state,

src/equation_of_state.hpp

@@ -71,6 +71,9 @@ class EquationOfState {
 
   double ComputePressure(const Vector &state, int dim);
 
+  void GetPrimitivesFromConservatives(const Vector &conserv, Vector &primit, const int &dim, const int &num_equations);
+  void GetConservativesFromPrimitives(const Vector &primit, Vector &conserv, const int &dim, const int &num_equations);
+
   // Compute the maximum characteristic speed.
   double ComputeMaxCharSpeed(const Vector &state, const int dim);
 

src/forcing_terms.cpp

@@ -32,6 +32,7 @@
 #include "forcing_terms.hpp"
 
 #include <general/forall.hpp>
+#include <vector>
 
 ForcingTerms::ForcingTerms(const int &_dim, const int &_num_equation, const int &_order, const int &_intRuleType,
                            IntegrationRules *_intRules, ParFiniteElementSpace *_vfes, ParGridFunction *_Up,
@@ -213,6 +214,159 @@ void ConstantPressureGradient::updateTerms_gpu(const int numElems, const int off
 
 #endif
 
+SpongeZone::SpongeZone(const int &_dim, const int &_num_equation, const int &_order, const int &_intRuleType,
+                       Fluxes *_fluxClass, EquationOfState *_eqState, IntegrationRules *_intRules,
+                       ParFiniteElementSpace *_vfes, ParGridFunction *_Up, ParGridFunction *_gradUp,
+                       const volumeFaceIntegrationArrays &gpuArrays, RunConfiguration &_config)
+    : ForcingTerms(_dim, _num_equation, _order, _intRuleType, _intRules, _vfes, _Up, _gradUp, gpuArrays),
+      fluxes(_fluxClass),
+      eqState(_eqState),
+      szData(_config.GetSpongeZoneData()) {
+  targetU.SetSize(num_equation);
+  if (szData.szType == SpongeZoneSolution::USERDEF) {
+    Vector Up(num_equation);
+    Up[0] = szData.targetUp[0];
+    for (int d = 0; d < dim; d++) Up[1 + d] = szData.targetUp[1 + d];
+    Up[1 + dim] = szData.targetUp[4];
+    eqState->GetConservativesFromPrimitives(Up, targetU, dim, num_equation);
+  }
+
+  meanNormalFluxes.SetSize(num_equation + 1);
+
+  ParMesh *mesh = vfes->GetParMesh();
+  const FiniteElementCollection *fec = vfes->FEColl();
+  ParFiniteElementSpace dfes(mesh, fec, dim, Ordering::byNODES);
+  ParFiniteElementSpace fes(mesh, fec);
+
+  sigma = new ParGridFunction(&fes);
+  *sigma = 0.;
+  double *hSigma = sigma->HostWrite();
+
+  ParGridFunction coords(&dfes);
+  mesh->GetNodes(coords);
+  int ndofs = vfes->GetNDofs();
+
+  vector<int> nodesVec;
+  nodesVec.clear();
+  for (int n = 0; n < ndofs; n++) {
+    Vector Xn(dim);
+    for (int d = 0; d < dim; d++) Xn[d] = coords[n + d * ndofs];
+
+    // distance to the mix-out plane
+    double dist = 0.;
+    for (int d = 0; d < dim; d++) dist += szData.normal[d] * (Xn[d] - szData.pointInit[d]);
+
+    if (fabs(dist) < szData.tol) nodesVec.push_back(n);
+
+    // dist end plane
+    double distF = 0.;
+    for (int d = 0; d < dim; d++) distF += szData.normal[d] * (Xn[d] - szData.point0[d]);
+
+    if (dist < 0. && distF > 0.) {
+      double planeDistance = distF - dist;
+      hSigma[n] = -dist / planeDistance / planeDistance;
+    }
+  }
+
+  // find plane nodes
+  if (szData.szType == SpongeZoneSolution::MIXEDOUT) {
+    nodesInMixedOutPlane.SetSize(nodesVec.size());
+    for (int n = 0; n < nodesInMixedOutPlane.Size(); n++) nodesInMixedOutPlane[n] = nodesVec[n];
+
+    cout << "__________________________________________________" << endl;
+    cout << "Running with Mixed-Out Sponge Zone" << endl;
+    cout << "       Nodes searched with tolerance: " << szData.tol << endl;
+    cout << "       Num. of nodes in mix-out plane: " << nodesInMixedOutPlane.Size() << endl;
+    cout << "^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^" << endl;
+  }
+}
+
+SpongeZone::~SpongeZone() { delete sigma; }
+
+void SpongeZone::updateTerms(Vector &in) {
+  if (szData.szType == SpongeZoneSolution::MIXEDOUT) computeMixedOutValues();
+
+  addSpongeZoneForcing(in);
+}
+
+void SpongeZone::addSpongeZoneForcing(Vector &in) {
+  const double *ds = sigma->HostRead();
+  double *dataIn = in.HostReadWrite();
+  const double *dataUp = Up->HostRead();
+
+  int nnodes = vfes->GetNDofs();
+
+  Vector Un(num_equation), Up(num_equation);
+
+  // compute speed of sound
+  double gamma = eqState->GetSpecificHeatRatio();
+  eqState->GetPrimitivesFromConservatives(targetU, Up, dim, num_equation);
+  double speedSound = sqrt(gamma * Up[num_equation - 1] / Up[0]);
+
+  // add forcing to RHS, i.e., @in
+  for (int n = 0; n < nnodes; n++) {
+    double s = ds[n];
+    if (s > 0.) {
+      s *= szData.multFactor;
+      for (int eq = 0; eq < num_equation; eq++) Up[eq] = dataUp[n + eq * nnodes];
+      eqState->GetConservativesFromPrimitives(Up, Un, dim, num_equation);
+
+      for (int eq = 0; eq < num_equation; eq++) dataIn[n + eq * nnodes] -= speedSound * s * (Un[eq] - targetU[eq]);
+    }
+  }
+}
+
+void SpongeZone::computeMixedOutValues() {
+  int nnodes = vfes->GetNDofs();
+  const double *dataUp = Up->HostRead();
+  double gamma = eqState->GetSpecificHeatRatio();
+
+  // compute mean normal fluxes
+  meanNormalFluxes = 0.;
+
+  Vector Un(num_equation), Up(num_equation);
+  DenseMatrix f(num_equation, dim);
+  for (int n = 0; n < nodesInMixedOutPlane.Size(); n++) {
+    int node = nodesInMixedOutPlane[n];
+
+    for (int eq = 0; eq < num_equation; eq++) Up[eq] = dataUp[node + eq * nnodes];
+    eqState->GetConservativesFromPrimitives(Up, Un, dim, num_equation);
+    fluxes->ComputeConvectiveFluxes(Un, f);
+
+    for (int eq = 0; eq < num_equation; eq++)
+      for (int d = 0; d < dim; d++) meanNormalFluxes[eq] += szData.normal[d] * f(eq, d);
+  }
+  meanNormalFluxes[num_equation] = double(nodesInMixedOutPlane.Size());
+
+  // communicate different partitions
+  Vector sum(num_equation + 1);
+  MPI_Allreduce(meanNormalFluxes.HostRead(), sum.HostWrite(), num_equation + 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
+  meanNormalFluxes = sum;
+
+  // normalize
+  for (int eq = 0; eq < num_equation; eq++) meanNormalFluxes[eq] /= double(meanNormalFluxes[num_equation]);
+
+  // compute mixed-out variables considering the outlet pressure
+  double temp = 0.;
+  for (int d = 0; d < dim; d++) temp += meanNormalFluxes[1 + d] * szData.normal[d];
+  double A = 1. - 2. * gamma / (gamma - 1.);
+  double B = 2 * temp / (gamma - 1.);
+  double C = -2. * meanNormalFluxes[0] * meanNormalFluxes[num_equation - 1];
+  for (int d = 0; d < dim; d++) C += meanNormalFluxes[1 + d] * meanNormalFluxes[1 + d];
+  //   double p = (-B+sqrt(B*B-4.*A*C))/(2.*A);
+  double p = (-B - sqrt(B * B - 4. * A * C)) / (2. * A);  // real solution
+
+  Up[0] = meanNormalFluxes[0] * meanNormalFluxes[0] / (temp - p);
+  Up[num_equation - 1] = p;
+
+  for (int d = 0; d < dim; d++) Up[1 + d] = (meanNormalFluxes[1 + d] - p * szData.normal[d]) / meanNormalFluxes[0];
+
+  double v2 = 0.;
+  for (int d = 0; d < dim; d++) v2 += Up[1 + d] * Up[1 + d];
+
+  eqState->GetConservativesFromPrimitives(Up, targetU, dim, num_equation);
+}
+
 #ifdef _MASA_
 MASA_forcings::MASA_forcings(const int &_dim, const int &_num_equation, const int &_order, const int &_intRuleType,
                              IntegrationRules *_intRules, ParFiniteElementSpace *_vfes, ParGridFunction *_Up,

src/forcing_terms.hpp

@@ -37,6 +37,8 @@
 #include <mfem.hpp>
 
 #include "dataStructures.hpp"
+#include "equation_of_state.hpp"
+#include "fluxes.hpp"
 #include "run_configuration.hpp"
 
 #ifdef _MASA_
@@ -100,6 +102,32 @@ class ConstantPressureGradient : public ForcingTerms {
 #endif
 };
 
+class SpongeZone : public ForcingTerms {
+ private:
+  Fluxes *fluxes;
+  EquationOfState *eqState;
+
+  SpongeZoneData &szData;
+  Vector targetU;
+
+  Array<int> nodesInMixedOutPlane;
+
+  ParGridFunction *sigma;  // linearly varying factor
+
+  Vector meanNormalFluxes;
+
+  void computeMixedOutValues();
+  void addSpongeZoneForcing(Vector &in);
+
+ public:
+  SpongeZone(const int &_dim, const int &_num_equation, const int &_order, const int &_intRuleType, Fluxes *_fluxClass,
+             EquationOfState *_eqState, IntegrationRules *_intRules, ParFiniteElementSpace *_vfes, ParGridFunction *_Up,
+             ParGridFunction *_gradUp, const volumeFaceIntegrationArrays &gpuArrays, RunConfiguration &_config);
+  ~SpongeZone();
+
+  virtual void updateTerms(Vector &in);
+};
+
 #ifdef _MASA_
 // Manufactured Solution using MASA
 class MASA_forcings : public ForcingTerms {

src/gradients.hpp

@@ -35,6 +35,7 @@
 // Class to manage gradients of primitives
 
 #include <tps_config.h>
+
 #include <general/forall.hpp>
 #include <mfem.hpp>