Add option specifying whether to include BC terms in gradient calc (#198

) Eventually useBCinGrad = true will be the default (and only) option. But, during development, it is useful to be able to turn it off so we can ensure we don't introduce any unintended effects by running the regression tests (which will have to have their 'truth' solutions regenerated once this change is ready.
pecos · Jul 18, 2023 · 1db8986 · 1db8986
1 parent 0acbc69
commit 1db8986
Show file tree

Hide file tree

Showing 4 changed files with 28 additions and 16 deletions.
diff --git a/src/M2ulPhyS.cpp b/src/M2ulPhyS.cpp
@@ -606,7 +606,7 @@ void M2ulPhyS::initVariables() {
 
   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));
+  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,
@@ -2980,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);

diff --git a/src/faceGradientIntegration.cpp b/src/faceGradientIntegration.cpp
@@ -34,8 +34,8 @@
 
 // Implementation of class FaceIntegrator
 GradFaceIntegrator::GradFaceIntegrator(IntegrationRules *_intRules, const int _dim, const int _num_equation,
-                                       BCintegrator *bc)
-    : dim(_dim), num_equation(_num_equation), intRules(_intRules), bc_(bc) {}
+                                       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) {
@@ -96,18 +96,22 @@ void GradFaceIntegrator::AssembleFaceVector(const FiniteElement &el1, const Fini
     if (Tr.Elem2No < 0) {
       assert(bc_ != NULL);
 
-      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);
+      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;
       }
     } else {
       elfun2_mat.MultTranspose(shape2, iUp2);

diff --git a/src/faceGradientIntegration.hpp b/src/faceGradientIntegration.hpp
@@ -45,10 +45,13 @@ class GradFaceIntegrator : public NonlinearFormIntegrator {
   const int dim;
   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, BCintegrator *bc = NULL);
+  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);

diff --git a/src/run_configuration.hpp b/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;