From 6b68c92e9f9061fdc54309216b56a77b7ec73cd1 Mon Sep 17 00:00:00 2001 From: "Todd A. Oliver" Date: Tue, 21 Feb 2023 16:31:49 -0600 Subject: [PATCH] Some logic to be able to use BC info in GradFaceIntegrator (#198) The idea here is to take advantage of all the data that has already been set up in class BCIntegrator by making class GradFaceIntegrator its friend. Then we can use the internal std::unordered_map objects in BCIntegrator, which map the boundary attribute to the appropriate BoundaryCondition class, to get access to any required boundary condition information. In principle that is all that is necessary, but the data we need it is not in a convenient form. To address this, the method computeBdrPrimitiveStateForGradient has been added to the BoundaryCondition class. In the base class, this method just sets the boundary state to the internal state, which results in the same (incorrect) behavior we have now. In future commits, this method will be implemented in the children of BoundaryCondition in order to provide the right boundary state. --- src/BCintegrator.hpp | 2 ++ src/BoundaryCondition.cpp | 4 ++++ src/BoundaryCondition.hpp | 10 ++++++++++ src/M2ulPhyS.cpp | 13 +++---------- src/faceGradientIntegration.cpp | 32 ++++++++++++++++++++++++++++---- src/faceGradientIntegration.hpp | 4 +++- 6 files changed, 50 insertions(+), 15 deletions(-) diff --git a/src/BCintegrator.hpp b/src/BCintegrator.hpp index 18d3ececc..0c95c114f 100644 --- a/src/BCintegrator.hpp +++ b/src/BCintegrator.hpp @@ -47,6 +47,8 @@ using namespace mfem; // Boundary face term: class BCintegrator : public NonlinearFormIntegrator { + friend class GradFaceIntegrator; + protected: MPI_Groups *groupsMPI; diff --git a/src/BoundaryCondition.cpp b/src/BoundaryCondition.cpp index d77530aa9..56f9c2b82 100644 --- a/src/BoundaryCondition.cpp +++ b/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; diff --git a/src/BoundaryCondition.hpp b/src/BoundaryCondition.hpp index da9961cd6..4ca7f2c99 100644 --- a/src/BoundaryCondition.hpp +++ b/src/BoundaryCondition.hpp @@ -76,6 +76,16 @@ class BoundaryCondition { virtual void computeBdrFlux(Vector &normal, Vector &stateIn, DenseMatrix &gradState, double radius, Vector transip, double delta, 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; diff --git a/src/M2ulPhyS.cpp b/src/M2ulPhyS.cpp index 2d9cb3103..fdf976d34 100644 --- a/src/M2ulPhyS.cpp +++ b/src/M2ulPhyS.cpp @@ -598,17 +598,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)); rhsOperator = new RHSoperator(iter, dim, num_equation, order, eqSystem, max_char_speed, intRules, intRuleType, d_fluxClass, diff --git a/src/faceGradientIntegration.cpp b/src/faceGradientIntegration.cpp index dbc873106..790a43660 100644 --- a/src/faceGradientIntegration.cpp +++ b/src/faceGradientIntegration.cpp @@ -33,8 +33,9 @@ #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) + : dim(_dim), num_equation(_num_equation), intRules(_intRules), bc_(bc) {} void GradFaceIntegrator::AssembleFaceVector(const FiniteElement &el1, const FiniteElement &el2, FaceElementTransformations &Tr, const Vector &elfun, Vector &elvect) { @@ -76,6 +77,7 @@ void GradFaceIntegrator::AssembleFaceVector(const FiniteElement &el1, const Fini } 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); @@ -83,11 +85,33 @@ void GradFaceIntegrator::AssembleFaceVector(const FiniteElement &el1, const Fini // 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 elfun1_mat.MultTranspose(shape1, iUp1); - elfun2_mat.MultTranspose(shape2, iUp2); + if (Tr.Elem2No < 0) { + assert(bc_ != NULL); + + const int attr = Tr.Attribute; + std::unordered_map::const_iterator ibc = bc_->inletBCmap.find(attr); + std::unordered_map::const_iterator obc = bc_->outletBCmap.find(attr); + std::unordered_map::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 { + elfun2_mat.MultTranspose(shape2, iUp2); + } // Compute average // NB: Code below is mathematically equivalent to mean = 0.5*(iUp1 diff --git a/src/faceGradientIntegration.hpp b/src/faceGradientIntegration.hpp index 134839bdd..61ebfdef4 100644 --- a/src/faceGradientIntegration.hpp +++ b/src/faceGradientIntegration.hpp @@ -34,6 +34,7 @@ #include +#include "BCintegrator.hpp" #include "tps_mfem_wrap.hpp" using namespace mfem; @@ -44,9 +45,10 @@ class GradFaceIntegrator : public NonlinearFormIntegrator { const int dim; const int num_equation; IntegrationRules *intRules; + BCintegrator *bc_; // NB: GradFaceIntegrator is a friend of BCintegrator public: - GradFaceIntegrator(IntegrationRules *_intRules, const int _dim, const int _num_equation); + GradFaceIntegrator(IntegrationRules *_intRules, const int _dim, const int _num_equation, BCintegrator *bc = NULL); virtual void AssembleFaceVector(const FiniteElement &el1, const FiniteElement &el2, FaceElementTransformations &Tr, const Vector &elfun, Vector &elvect);