From a720e333e9006adc70225b406ecefc4c6249aa4b Mon Sep 17 00:00:00 2001 From: "Todd A. Oliver" Date: Tue, 21 Feb 2023 13:00:59 -0600 Subject: [PATCH 1/6] Eliminate some dead `_GPU_` code. Specifically, parts of GradNonLinearForm and GradFaceIntegrator had _GPU_ specific code that was superseded by device-specific code in the Gradients class and thus was never called. This code has been eliminated to improve hygiene in preparation for changes to address #198. --- src/faceGradientIntegration.cpp | 31 +++---------- src/gradNonLinearForm.cpp | 78 +++------------------------------ src/gradNonLinearForm.hpp | 6 --- 3 files changed, 12 insertions(+), 103 deletions(-) diff --git a/src/faceGradientIntegration.cpp b/src/faceGradientIntegration.cpp index 6d67cce56..2b19e2f67 100644 --- a/src/faceGradientIntegration.cpp +++ b/src/faceGradientIntegration.cpp @@ -60,18 +60,17 @@ void GradFaceIntegrator::AssembleFaceVector(const FiniteElement &el1, const Fini 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; @@ -104,20 +103,12 @@ void GradFaceIntegrator::AssembleFaceVector(const FiniteElement &el1, const Fini 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 } iUp2(eq) = sum; mean(eq) = (iUp1(eq) + iUp2(eq)) / 2.; @@ -133,18 +124,10 @@ void GradFaceIntegrator::AssembleFaceVector(const FiniteElement &el1, const Fini 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 } } } diff --git a/src/gradNonLinearForm.cpp b/src/gradNonLinearForm.cpp index 87ad42ada..243a645a1 100644 --- a/src/gradNonLinearForm.cpp +++ b/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 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 &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_ diff --git a/src/gradNonLinearForm.hpp b/src/gradNonLinearForm.hpp index 6be197070..7ec050d3a 100644 --- a/src/gradNonLinearForm.hpp +++ b/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 &vdofs, Vector &y, const Vector &el_y, const int &dof, - const int &elDof, const int &num_equation, const int &dim); -#endif }; #endif // GRADNONLINEARFORM_HPP_ From 032c82af8e1201e5f438d897ec469378ce58f2a2 Mon Sep 17 00:00:00 2001 From: "Todd A. Oliver" Date: Tue, 21 Feb 2023 13:30:55 -0600 Subject: [PATCH 2/6] Some clean up in class GradFaceIntegrator In preparation for further refactoring as part of addressing #198 --- src/faceGradientIntegration.cpp | 61 ++++++++++----------------------- 1 file changed, 19 insertions(+), 42 deletions(-) diff --git a/src/faceGradientIntegration.cpp b/src/faceGradientIntegration.cpp index 2b19e2f67..dbc873106 100644 --- a/src/faceGradientIntegration.cpp +++ b/src/faceGradientIntegration.cpp @@ -39,24 +39,17 @@ GradFaceIntegrator::GradFaceIntegrator(IntegrationRules *_intRules, const int _d void GradFaceIntegrator::AssembleFaceVector(const FiniteElement &el1, const FiniteElement &el2, FaceElementTransformations &Tr, const Vector &elfun, Vector &elvect) { // Compute the term 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; @@ -66,7 +59,6 @@ void GradFaceIntegrator::AssembleFaceVector(const FiniteElement &el1, const Fini 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.; @@ -82,7 +74,6 @@ 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); for (int i = 0; i < ir->GetNPoints(); i++) { @@ -95,41 +86,27 @@ void GradFaceIntegrator::AssembleFaceVector(const FiniteElement &el1, const Fini 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++) { - sum += shape1[k] * elfun1_mat(k, eq); - } - iUp1(eq) = sum; - sum = 0.; - for (int k = 0; k < dof2; k++) { - sum += shape2[k] * elfun2_mat(k, eq); - } - iUp2(eq) = sum; - mean(eq) = (iUp1(eq) + iUp2(eq)) / 2.; - } + elfun1_mat.MultTranspose(shape1, iUp1); + 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++) { - elvect1_mat(k, eq + d * num_equation) += du1n * shape1(k); - } - for (int k = 0; k < dof2; k++) { - elvect2_mat(k, eq + d * num_equation) -= du2n * shape2(k); - } + 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); } } From 0acbc69464fe6e421b331e1be9ccbd6653ce7ade Mon Sep 17 00:00:00 2001 From: "Todd A. Oliver" Date: Tue, 21 Feb 2023 16:31:49 -0600 Subject: [PATCH 3/6] 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 aea131573..0a001142e 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 f95bc89e9..5ff2364c1 100644 --- a/src/BoundaryCondition.hpp +++ b/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; diff --git a/src/M2ulPhyS.cpp b/src/M2ulPhyS.cpp index 9082c388b..e30248d54 100644 --- a/src/M2ulPhyS.cpp +++ b/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)); 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); From 1db8986520bf6d80252978cf1667aaca5b3d2141 Mon Sep 17 00:00:00 2001 From: "Todd A. Oliver" Date: Tue, 21 Feb 2023 18:59:29 -0600 Subject: [PATCH 4/6] 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. --- src/M2ulPhyS.cpp | 4 +++- src/faceGradientIntegration.cpp | 32 ++++++++++++++++++-------------- src/faceGradientIntegration.hpp | 5 ++++- src/run_configuration.hpp | 3 +++ 4 files changed, 28 insertions(+), 16 deletions(-) diff --git a/src/M2ulPhyS.cpp b/src/M2ulPhyS.cpp index e30248d54..97d6fed0d 100644 --- 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 index 790a43660..0f0061edd 100644 --- 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::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); + if (useBCinGrad_) { + 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 { + iUp2 = iUp1; } } else { elfun2_mat.MultTranspose(shape2, iUp2); diff --git a/src/faceGradientIntegration.hpp b/src/faceGradientIntegration.hpp index 61ebfdef4..f2a16f0a3 100644 --- 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 index cb8246071..3764b65ce 100644 --- a/src/run_configuration.hpp +++ b/src/run_configuration.hpp @@ -185,6 +185,9 @@ class RunConfiguration { std::vector wallBC; std::vector> wallPatchType; + // Apply BCs in gradient calculation + bool useBCinGrad; + // Passive scalar data Array arrayPassiveScalar; From fbce85126726924959a321414a3f803d1f05961a Mon Sep 17 00:00:00 2001 From: "Todd A. Oliver" Date: Tue, 21 Feb 2023 19:02:41 -0600 Subject: [PATCH 5/6] Towards #198 --- Add primitive BC for isothermal wall --- src/wallBC.cpp | 24 ++++++++++++++++++++++++ src/wallBC.hpp | 1 + 2 files changed, 25 insertions(+) diff --git a/src/wallBC.cpp b/src/wallBC.cpp index 90c605d75..d07adad1d 100644 --- a/src/wallBC.cpp +++ b/src/wallBC.cpp @@ -236,6 +236,30 @@ void WallBC::computeBdrFlux(Vector &normal, Vector &stateIn, DenseMatrix &gradSt } } +void WallBC::computeBdrPrimitiveStateForGradient(const Vector &primIn, Vector &primBC) const { + primBC = primIn; + + switch (wallType_) { + case INV: + // TODO(trevilo): fix + break; + case VISC_ADIAB: + // TODO(trevilo): fix + break; + case VISC_ISOTH: + // no-slip + for (int i = 0; i < nvel_; i++) { + primBC[1 + i] = 0.0; + } + // isothermal + primBC[nvel_ + 1] = wallTemp_; + break; + case VISC_GNRL: + // TODO(trevilo): fix + break; + } +} + void WallBC::integrationBC(Vector &y, const Vector &x, const elementIndexingData &elem_index_data, ParGridFunction *Up, ParGridFunction *gradUp, const boundaryFaceIntegrationData &boundary_face_data, const int &maxIntPoints, const int &maxDofs) { diff --git a/src/wallBC.hpp b/src/wallBC.hpp index a277b81aa..c6cc10d00 100644 --- a/src/wallBC.hpp +++ b/src/wallBC.hpp @@ -83,6 +83,7 @@ class WallBC : public BoundaryCondition { void computeBdrFlux(Vector &normal, Vector &stateIn, DenseMatrix &gradState, Vector transip, double delta, double distance, Vector &bdrFlux); + void computeBdrPrimitiveStateForGradient(const Vector &primIn, Vector &primBC) const override; virtual void initBCs(); From fa158f20190f542672a72a869fe689c32efef097 Mon Sep 17 00:00:00 2001 From: "Todd A. Oliver" Date: Wed, 22 Feb 2023 08:50:05 -0600 Subject: [PATCH 6/6] Switch isothermal wall BC to mass-conserving variant (but only when using recent mods to account for BC in gradient) --- src/BCintegrator.cpp | 3 ++- src/wallBC.cpp | 19 ++++++++++++++----- src/wallBC.hpp | 5 ++++- 3 files changed, 20 insertions(+), 7 deletions(-) diff --git a/src/BCintegrator.cpp b/src/BCintegrator.cpp index 8607e6eab..d06a0eca1 100644 --- a/src/BCintegrator.cpp +++ b/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); } } diff --git a/src/wallBC.cpp b/src/wallBC.cpp index d07adad1d..85271b433 100644 --- a/src/wallBC.cpp +++ b/src/wallBC.cpp @@ -37,7 +37,8 @@ WallBC::WallBC(RiemannSolver *_rsolver, GasMixture *_mixture, GasMixture *d_mixture, Equations _eqSystem, Fluxes *_fluxClass, ParFiniteElementSpace *_vfes, IntegrationRules *_intRules, double &_dt, const int _dim, const int _num_equation, int _patchNumber, WallType _bcType, const WallData _inputData, - const boundaryFaceIntegrationData &boundary_face_data, const int &_maxIntPoints, bool axisym) + const boundaryFaceIntegrationData &boundary_face_data, const int &_maxIntPoints, bool axisym, + bool useBCinGrad) : BoundaryCondition(_rsolver, _mixture, _eqSystem, _vfes, _intRules, _dt, _dim, _num_equation, _patchNumber, 1, axisym), // so far walls do not require ref. length. Left at 1 wallType_(_bcType), @@ -45,7 +46,8 @@ WallBC::WallBC(RiemannSolver *_rsolver, GasMixture *_mixture, GasMixture *d_mixt d_mixture_(d_mixture), fluxClass(_fluxClass), boundary_face_data_(boundary_face_data), - maxIntPoints_(_maxIntPoints) { + maxIntPoints_(_maxIntPoints), + useBCinGrad_(useBCinGrad) { // Initialize bc state. // bcState_.prim.SetSize(num_equation_); // bcState_.primIdxs.SetSize(num_equation_); @@ -470,10 +472,15 @@ void WallBC::computeAdiabaticWallFlux(Vector &normal, Vector &stateIn, DenseMatr void WallBC::computeIsothermalWallFlux(Vector &normal, Vector &stateIn, DenseMatrix &gradState, Vector transip, double delta, Vector &bdrFlux) { Vector wallState(num_equation_); - mixture->computeStagnantStateWithTemp(stateIn, wallTemp_, wallState); - // TODO(kevin): set stangant state with two separate temperature. + wallState = stateIn; - if (eqSystem == NS_PASSIVE) wallState[num_equation_ - 1] = stateIn[num_equation_ - 1]; + if (useBCinGrad_) { + for (int i = 0; i < nvel_; i++) { + wallState[i + 1] *= -1.0; + } + } else { + mixture->computeStagnantStateWithTemp(stateIn, wallTemp_, wallState); + } // Normal convective flux rsolver->Eval(stateIn, wallState, normal, bdrFlux, true); @@ -486,6 +493,8 @@ void WallBC::computeIsothermalWallFlux(Vector &normal, Vector &stateIn, DenseMat for (int d = 0; d < dim_; d++) bcFlux_.normal[d] = unitNorm[d]; // evaluate viscous fluxes at the wall + wallState = stateIn; + mixture->computeStagnantStateWithTemp(stateIn, wallTemp_, wallState); Vector wallViscF(num_equation_); fluxClass->ComputeBdrViscousFluxes(wallState, gradState, transip, delta, 0.0, bcFlux_, wallViscF); wallViscF *= sqrt(normN); // in case normal is not a unit vector.. diff --git a/src/wallBC.hpp b/src/wallBC.hpp index c6cc10d00..d970544bc 100644 --- a/src/wallBC.hpp +++ b/src/wallBC.hpp @@ -47,6 +47,8 @@ class WallBC : public BoundaryCondition { const WallType wallType_; const WallData wallData_; + const bool useBCinGrad_; + GasMixture *d_mixture_; // only used in the device. Fluxes *fluxClass; @@ -78,7 +80,8 @@ class WallBC : public BoundaryCondition { WallBC(RiemannSolver *rsolver_, GasMixture *_mixture, GasMixture *d_mixture, Equations _eqSystem, Fluxes *_fluxClass, ParFiniteElementSpace *_vfes, IntegrationRules *_intRules, double &_dt, const int _dim, const int _num_equation, int _patchNumber, WallType _bcType, const WallData _inputData, - const boundaryFaceIntegrationData &boundary_face_data, const int &maxIntPoints, bool axisym); + const boundaryFaceIntegrationData &boundary_face_data, const int &maxIntPoints, bool axisym, + bool useBCinGrad = false); ~WallBC(); void computeBdrFlux(Vector &normal, Vector &stateIn, DenseMatrix &gradState, Vector transip, double delta,