Toward boundary gradient fix: part 1

src/BCintegrator.cpp

@@ -115,7 +115,8 @@ BCintegrator::BCintegrator(MPI_Groups *_groupsMPI, ParMesh *_mesh, ParFiniteElem
 
       wallBCmap[patchType.first] = new WallBC(rsolver, mixture, d_mixture, _runFile.GetEquationSystem(), fluxClass,
                                               vfes, intRules, _dt, dim, num_equation, patchType.first, patchType.second,
-                                              wallData, boundary_face_data_, _maxIntPoints, config.isAxisymmetric());
+                                              wallData, boundary_face_data_, _maxIntPoints, config.isAxisymmetric(),
+                                              config.useBCinGrad);
     }
   }
 

src/BCintegrator.hpp

@@ -47,6 +47,8 @@ using namespace mfem;
 
 // Boundary face term: <F.n(u),[w]>
 class BCintegrator : public NonlinearFormIntegrator {
+  friend class GradFaceIntegrator;
+
  protected:
   MPI_Groups *groupsMPI;
 

src/BoundaryCondition.cpp

@@ -52,6 +52,10 @@ BoundaryCondition::BoundaryCondition(RiemannSolver *_rsolver, GasMixture *_mixtu
 
 BoundaryCondition::~BoundaryCondition() {}
 
+void BoundaryCondition::computeBdrPrimitiveStateForGradient(const Vector &stateIn, Vector &stateBC) const {
+  stateBC = stateIn;
+}
+
 double BoundaryCondition::aggregateArea(int bndry_patchnum, MPI_Comm bc_comm) {
   double area = 0;
 

src/BoundaryCondition.hpp

@@ -76,6 +76,16 @@ class BoundaryCondition {
   virtual void computeBdrFlux(Vector &normal, Vector &stateIn, DenseMatrix &gradState, Vector transip, double delta,
                               double distance, Vector &bdrFlux) = 0;
 
+  /** \brief Set the boundary state used in the gradient evaluation
+   *
+   *  The jump in the state at the boundary appears in the
+   *  right-hand-side of the gradient solve.  This method sets this
+   *  boundary state to be equal to the interior state such that this
+   *  term is zero.  If that is not what you want, you must override
+   *  this method in a derived class.
+   */
+  virtual void computeBdrPrimitiveStateForGradient(const Vector &stateIn, Vector &stateBC) const;
+
   // holding function for any miscellaneous items needed to initialize BCs
   // prior to use (and require MPI)
   virtual void initBCs() = 0;

src/M2ulPhyS.cpp

@@ -603,17 +603,10 @@ void M2ulPhyS::initVariables() {
       cout << "Unknown ODE solver type: " << config.GetTimeIntegratorType() << '\n';
   }
 
-  //   gradUp_A = new ParNonlinearForm(gradUpfes);
-  //   gradUp_A->AddInteriorFaceIntegrator(
-  //       new GradFaceIntegrator(intRules, dim, num_equation) );
-  gradUp_A = new GradNonLinearForm(
-#ifdef _GPU_
-      vfes,
-#else
-      gradUpfes,
-#endif
-      intRules, dim, num_equation);
+
+  gradUp_A = new GradNonLinearForm(gradUpfes, intRules, dim, num_equation);
   gradUp_A->AddInteriorFaceIntegrator(new GradFaceIntegrator(intRules, dim, num_equation));
+  gradUp_A->AddBdrFaceIntegrator(new GradFaceIntegrator(intRules, dim, num_equation, bcIntegrator, config.useBCinGrad));
 
   rhsOperator =
       new RHSoperator(iter, dim, num_equation, order, eqSystem, max_char_speed, intRules, intRuleType, d_fluxClass,
@@ -2987,6 +2980,8 @@ void M2ulPhyS::parseReactionInputs() {
 }
 
 void M2ulPhyS::parseBCInputs() {
+  tpsP->getInput("boundaryConditions/useBCinGrad", config.useBCinGrad, false);
+
   // number of BC regions defined
   int numWalls, numInlets, numOutlets;
   tpsP->getInput("boundaryConditions/numWalls", numWalls, 0);

src/faceGradientIntegration.cpp

@@ -33,45 +33,37 @@
 #include "faceGradientIntegration.hpp"
 
 // Implementation of class FaceIntegrator
-GradFaceIntegrator::GradFaceIntegrator(IntegrationRules *_intRules, const int _dim, const int _num_equation)
-    : dim(_dim), num_equation(_num_equation), intRules(_intRules) {}
+GradFaceIntegrator::GradFaceIntegrator(IntegrationRules *_intRules, const int _dim, const int _num_equation,
+                                       BCintegrator *bc, bool useBCinGrad)
+    : dim(_dim), num_equation(_num_equation), intRules(_intRules), bc_(bc), useBCinGrad_(useBCinGrad) {}
 
 void GradFaceIntegrator::AssembleFaceVector(const FiniteElement &el1, const FiniteElement &el2,
                                             FaceElementTransformations &Tr, const Vector &elfun, Vector &elvect) {
   // Compute the term <nU,[w]> on the interior faces.
-  Vector shape1;
-  shape1.UseDevice(false);
-  Vector shape2;
-  shape2.UseDevice(false);
   Vector nor(dim);
-  nor.UseDevice(false);
-  Vector mean;
-  mean.UseDevice(false);
-  mean.SetSize(num_equation);
+  Vector mean(num_equation);
+  Vector iUp1(num_equation), iUp2(num_equation);
+  Vector du1n(dim * num_equation), du2n(dim * num_equation);
 
   const int dof1 = el1.GetDof();
   const int dof2 = el2.GetDof();
 
-  shape1.SetSize(dof1);
-  shape2.SetSize(dof2);
+  Vector shape1(dof1);
+  Vector shape2(dof2);
 
-  elfun.Read();
-  elvect.UseDevice(true);
   elvect.SetSize((dof1 + dof2) * num_equation * dim);
   elvect = 0.0;
 
-// #ifdef _GPU_
-//   DenseMatrix elfun1_mat(elfun.GetData(), dof1, num_equation);
-//   DenseMatrix elfun2_mat(elfun.GetData()+dof1*num_equation,
-//                          dof2, num_equation);
-#ifndef _GPU_
-  DenseMatrix elfun1_mat(elfun.GetData(), dof1, num_equation * dim);
-  DenseMatrix elfun2_mat(elfun.GetData() + dof1 * num_equation * dim, dof2, num_equation * dim);
+  // NB: Incoming Vector &elfun is in gradient space.  The first
+  // component is the state, and the rest is zero.  See
+  // GradNonLinearForm::Mult in gradNonLinearForm.cpp for details.
+  DenseMatrix elfun1_mat(elfun.GetData(), dof1, num_equation);
+  DenseMatrix elfun2_mat(elfun.GetData() + dof1 * num_equation * dim, dof2, num_equation);
+
   DenseMatrix elvect1_mat(elvect.GetData(), dof1, num_equation * dim);
   DenseMatrix elvect2_mat(elvect.GetData() + dof1 * num_equation * dim, dof2, num_equation * dim);
   elvect1_mat = 0.;
   elvect2_mat = 0.;
-#endif
 
   // Integration order calculation from DGTraceIntegrator
   int intorder;
@@ -83,70 +75,66 @@ void GradFaceIntegrator::AssembleFaceVector(const FiniteElement &el1, const Fini
   if (el1.Space() == FunctionSpace::Pk) {
     intorder++;
   }
-  // IntegrationRules IntRules2(0, Quadrature1D::GaussLobatto);
   const IntegrationRule *ir = &intRules->Get(Tr.GetGeometryType(), intorder);
 
+  // Quadrature point loop
   for (int i = 0; i < ir->GetNPoints(); i++) {
     const IntegrationPoint &ip = ir->IntPoint(i);
 
     Tr.SetAllIntPoints(&ip);  // set face and element int. points
 
     // Calculate basis functions on both elements at the face
     el1.CalcShape(Tr.GetElement1IntPoint(), shape1);
-    el2.CalcShape(Tr.GetElement2IntPoint(), shape2);
+    if (Tr.Elem2No < 0) {
+      shape2 = shape1;
+    } else {
+      el2.CalcShape(Tr.GetElement2IntPoint(), shape2);
+    }
 
     // Interpolate U values at the point
-    Vector iUp1, iUp2;
-    iUp1.UseDevice(false);
-    iUp2.UseDevice(false);
-    iUp1.SetSize(num_equation);
-    iUp2.SetSize(num_equation);
-    for (int eq = 0; eq < num_equation; eq++) {
-      double sum = 0.;
-      for (int k = 0; k < dof1; k++) {
-#ifdef _GPU_
-        sum += shape1(k) * elfun(k + eq * dof1);
-#else
-        sum += shape1[k] * elfun1_mat(k, eq);
-#endif
-      }
-      iUp1(eq) = sum;
-      sum = 0.;
-      for (int k = 0; k < dof2; k++) {
-#ifdef _GPU_
-        sum += shape2(k) * elfun(k + eq * dof2 + dof1 * num_equation);
-#else
-        sum += shape2[k] * elfun2_mat(k, eq);
-#endif
+    elfun1_mat.MultTranspose(shape1, iUp1);
+    if (Tr.Elem2No < 0) {
+      assert(bc_ != NULL);
+
+      if (useBCinGrad_) {
+        const int attr = Tr.Attribute;
+        std::unordered_map<int, BoundaryCondition *>::const_iterator ibc = bc_->inletBCmap.find(attr);
+        std::unordered_map<int, BoundaryCondition *>::const_iterator obc = bc_->outletBCmap.find(attr);
+        std::unordered_map<int, BoundaryCondition *>::const_iterator wbc = bc_->wallBCmap.find(attr);
+        if (ibc != bc_->inletBCmap.end()) {
+          ibc->second->computeBdrPrimitiveStateForGradient(iUp1, iUp2);
+        }
+        if (obc != bc_->outletBCmap.end()) {
+          obc->second->computeBdrPrimitiveStateForGradient(iUp1, iUp2);
+        }
+        if (wbc != bc_->wallBCmap.end()) {
+          wbc->second->computeBdrPrimitiveStateForGradient(iUp1, iUp2);
+        }
+      } else {
+        iUp2 = iUp1;
       }
-      iUp2(eq) = sum;
-      mean(eq) = (iUp1(eq) + iUp2(eq)) / 2.;
+    } else {
+      elfun2_mat.MultTranspose(shape2, iUp2);
     }
 
-    // Get the normal vector and the flux on the face
-    CalcOrtho(Tr.Jacobian(), nor);
+    // Compute average
+    // NB: Code below is mathematically equivalent to mean = 0.5*(iUp1
+    // + iUp2), but mfem::Vector does not provide operator+
+    mean.Set(0.5, iUp1);
+    mean.Add(0.5, iUp2);
 
+    // Get the normal vector
+    CalcOrtho(Tr.Jacobian(), nor);
     nor *= ip.weight;
 
     for (int d = 0; d < dim; d++) {
       for (int eq = 0; eq < num_equation; eq++) {
-        const double du1n = (mean(eq) - iUp1(eq)) * nor(d);
-        const double du2n = (mean(eq) - iUp2(eq)) * nor(d);
-        for (int k = 0; k < dof1; k++) {
-#ifdef _GPU_
-          elvect(k + eq + d * num_equation) += du1n * shape1(k);
-#else
-          elvect1_mat(k, eq + d * num_equation) += du1n * shape1(k);
-#endif
-        }
-        for (int k = 0; k < dof2; k++) {
-#ifdef _GPU_
-          elvect(k + eq + d * num_equation + dof1 * num_equation * dim) -= du2n * shape2(k);
-#else
-          elvect2_mat(k, eq + d * num_equation) -= du2n * shape2(k);
-#endif
-        }
+        du1n[eq + d * num_equation] = (mean(eq) - iUp1(eq)) * nor(d);
+        du2n[eq + d * num_equation] = (iUp2(eq) - mean(eq)) * nor(d);
       }
     }
+
+    AddMultVWt(shape1, du1n, elvect1_mat);
+    AddMultVWt(shape2, du2n, elvect2_mat);
   }
 }

src/faceGradientIntegration.hpp

@@ -34,6 +34,7 @@
 
 #include <tps_config.h>
 
+#include "BCintegrator.hpp"
 #include "tps_mfem_wrap.hpp"
 
 using namespace mfem;
@@ -45,8 +46,12 @@ class GradFaceIntegrator : public NonlinearFormIntegrator {
   const int num_equation;
   IntegrationRules *intRules;
 
+  BCintegrator *bc_;  // NB: GradFaceIntegrator is a friend of BCintegrator
+  const bool useBCinGrad_;
+
  public:
-  GradFaceIntegrator(IntegrationRules *_intRules, const int _dim, const int _num_equation);
+  GradFaceIntegrator(IntegrationRules *_intRules, const int _dim, const int _num_equation, BCintegrator *bc = NULL,
+                     bool useBCinGrad = false);
 
   virtual void AssembleFaceVector(const FiniteElement &el1, const FiniteElement &el2, FaceElementTransformations &Tr,
                                   const Vector &elfun, Vector &elvect);

src/gradNonLinearForm.cpp

@@ -39,82 +39,14 @@ GradNonLinearForm::GradNonLinearForm(ParFiniteElementSpace *_vfes, IntegrationRu
 
 void GradNonLinearForm::Mult(const ParGridFunction *Up, Vector &y) {
   Vector x;
-
-#ifdef _GPU_
-  //   TestNonLinearForm::setToZero_gpu(y,y.Size()); already done in gradient.cpd
-  x.UseDevice(true);
-  if (fnfi.Size()) {
-    x.NewMemoryAndSize(Up->GetMemory(), vfes->GetNDofs() * num_equation, false);
-    Mult_gpu(x, y);
-  }
-#else
   const double *data = Up->GetData();
+
+  // NB: Setting x here accounts for the fact that the space Up lives
+  // in is different from the space of the gradient.
   x.SetSize(vfes->GetNDofs() * num_equation * dim);
   x = 0.;
   for (int i = 0; i < vfes->GetNDofs() * num_equation; i++) x(i) = data[i];
-  ParNonlinearForm::Mult(x, y);
-#endif
-}
-
-#ifdef _GPU_
-void GradNonLinearForm::Mult_gpu(const Vector &x, Vector &y) {
-  // INTEGRATION SHARED FACES
-  const Vector &px = Prolongate(x);
-
-  const int dofs = vfes->GetNDofs();
-
-  // Terms over shared interior faces in parallel.
-  ParFiniteElementSpace *pfes = ParFESpace();
-  ParMesh *pmesh = pfes->GetParMesh();
-  FaceElementTransformations *tr;
-  const FiniteElement *fe1, *fe2;
-  Array<int> vdofs1, vdofs2;
-  Vector el_x, el_y;
-  el_x.UseDevice(true);
-  el_y.UseDevice(true);
-
-  // if (bfnfi.Size()==0) y.HostReadWrite();
-  X.MakeRef(aux1, 0);  // aux1 contains P.x
-  X.ExchangeFaceNbrData();
-  const int n_shared_faces = pmesh->GetNSharedFaces();
-  for (int i = 0; i < n_shared_faces; i++) {
-    tr = pmesh->GetSharedFaceTransformations(i, true);
-    int Elem2NbrNo = tr->Elem2No - pmesh->GetNE();
 
-    fe1 = pfes->GetFE(tr->Elem1No);
-    fe2 = pfes->GetFaceNbrFE(Elem2NbrNo);
-
-    pfes->GetElementVDofs(tr->Elem1No, vdofs1);
-    pfes->GetFaceNbrElementVDofs(Elem2NbrNo, vdofs2);
-
-    el_x.SetSize(vdofs1.Size() + vdofs2.Size());
-    X.GetSubVector(vdofs1, el_x.GetData());
-    X.FaceNbrData().GetSubVector(vdofs2, el_x.GetData() + vdofs1.Size());
-
-    const int elDof1 = fe1->GetDof();
-    const int elDof2 = fe2->GetDof();
-    for (int k = 0; k < fnfi.Size(); k++) {
-      fnfi[k]->AssembleFaceVector(*fe1, *fe2, *tr, el_x, el_y);
-      // y.AddElementVector(vdofs1, el_y.GetData());
-      GradNonLinearForm::IndexedAddToGlobalMemory(vdofs1, y, el_y, dofs, elDof1, num_equation, dim);
-    }
-  }
-}
-
-void GradNonLinearForm::IndexedAddToGlobalMemory(const Array<int> &vdofs, Vector &y, const Vector &el_y, const int &dof,
-                                                 const int &elDof, const int &num_equation, const int &dim) {
-  double *dy = y.ReadWrite();
-  const double *d_ely = el_y.Read();
-  auto d_vdofs = vdofs.Read();
-
-  MFEM_FORALL(i, elDof, {
-    const int index = d_vdofs[i];
-    for (int eq = 0; eq < num_equation; eq++) {
-      for (int d = 0; d < dim; d++) {
-        dy[index + eq * dof + d * num_equation * dof] += d_ely[i + eq * elDof + d * num_equation * elDof];
-      }
-    }
-  });
+  // After setting x as above, base class Mult does the right thing
+  ParNonlinearForm::Mult(x, y);
 }
-
-#endif  // _GPU_

src/gradNonLinearForm.hpp

@@ -51,12 +51,6 @@ class GradNonLinearForm : public ParNonlinearForm {
   GradNonLinearForm(ParFiniteElementSpace *f, IntegrationRules *intRules, const int dim, const int num_equation);
 
   void Mult(const ParGridFunction *Up, Vector &y);
-
-#ifdef _GPU_
-  void Mult_gpu(const Vector &x, Vector &y);
-  static void IndexedAddToGlobalMemory(const Array<int> &vdofs, Vector &y, const Vector &el_y, const int &dof,
-                                       const int &elDof, const int &num_equation, const int &dim);
-#endif
 };
 
 #endif  // GRADNONLINEARFORM_HPP_

src/run_configuration.hpp

@@ -185,6 +185,9 @@ class RunConfiguration {
   std::vector<WallData> wallBC;
   std::vector<pair<int, WallType>> wallPatchType;
 
+  // Apply BCs in gradient calculation
+  bool useBCinGrad;
+
   // Passive scalar data
   Array<passiveScalarData*> arrayPassiveScalar;