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);