Skip to content

Commit

Permalink
baby steps on infrastructure to get grad BC fix into gpu path (#198)
Browse files Browse the repository at this point in the history
  • Loading branch information
@trevilo
trevilo committed May 11, 2023
1 parent 66773e4 commit 255899e
Show file tree
Hide file tree
Showing 13 changed files with 64 additions and 271 deletions.
14 changes: 10 additions & 4 deletions src/BCintegrator.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -111,10 +111,10 @@ BCintegrator::BCintegrator(MPI_Groups *_groupsMPI, ParMesh *_mesh, ParFiniteElem
if (patchInMesh) {
WallData wallData = config.GetWallData(w);

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(),
config.useBCinGrad);
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(), config.useBCinGrad);
}
}

Expand Down Expand Up @@ -264,6 +264,12 @@ void BCintegrator::integrateBCs(Vector &y, const Vector &x, const elementIndexin
}
}

void BCintegrator::integrateGradientBCs(Vector &y, const Vector &x, const elementIndexingData &elem_index_data) {
for (auto bc = wallBCmap.begin(); bc != wallBCmap.end(); bc++) {
bc->second->integrateGradientBC(y, x, elem_index_data, boundary_face_data_, maxIntPoints, maxDofs);
}
}

void BCintegrator::retrieveGradientsData_gpu(ParGridFunction *gradUp, DenseTensor &elGradUp, Array<int> &vdofs,
const int &num_equation, const int &dim, const int &totalDofs,
const int &elDofs) {
Expand Down
1 change: 1 addition & 0 deletions src/BCintegrator.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -101,6 +101,7 @@ class BCintegrator : public NonlinearFormIntegrator {

void updateBCMean(ParGridFunction *Up);
void integrateBCs(Vector &y, const Vector &x, const elementIndexingData &elem_index_data);
void integrateGradientBCs(Vector &y, const Vector &x, const elementIndexingData &elem_index_data);

// GPU functions
static void retrieveGradientsData_gpu(ParGridFunction *gradUp, DenseTensor &elGradUp, Array<int> &vdofs,
Expand Down
4 changes: 4 additions & 0 deletions src/BoundaryCondition.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -106,6 +106,10 @@ class BoundaryCondition {
const boundaryFaceIntegrationData &boundary_face_data, const int &maxIntPoints,
const int &maxDofs) = 0;

virtual void integrateGradientBC(Vector &y, const Vector &x, const elementIndexingData &elem_index_data,
const boundaryFaceIntegrationData &boundary_face_data, const int &maxIntPoints,
const int &maxDofs) {}

static void copyValues(const Vector &orig, Vector &target, const double &mult);
};

Expand Down
3 changes: 1 addition & 2 deletions src/M2ulPhyS.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -597,8 +597,7 @@ void M2ulPhyS::initVariables() {
cout << "Unknown ODE solver type: " << config.GetTimeIntegratorType() << '\n';
}


gradUp_A = new GradNonLinearForm(gradUpfes, intRules, dim, num_equation);
gradUp_A = new GradNonLinearForm(gradUpfes, intRules, dim, num_equation, bcIntegrator);
gradUp_A->AddInteriorFaceIntegrator(new GradFaceIntegrator(intRules, dim, num_equation));
gradUp_A->AddBdrFaceIntegrator(new GradFaceIntegrator(intRules, dim, num_equation, bcIntegrator, config.useBCinGrad));

Expand Down
14 changes: 6 additions & 8 deletions src/faceGradientIntegration.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -38,9 +38,7 @@ GradFaceIntegrator::GradFaceIntegrator(IntegrationRules *_intRules, const int _d
: 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) {

FaceElementTransformations &Tr, const Vector &elfun, Vector &elvect) {
// Compute the term <nU,[w]> on the interior faces.
Vector nor(dim);
Vector mean(num_equation);
Expand Down Expand Up @@ -105,15 +103,15 @@ void GradFaceIntegrator::AssembleFaceVector(const FiniteElement &el1, const Fini
std::unordered_map<int, BoundaryCondition *>::const_iterator wbc = bc_->wallBCmap.find(attr);
if (ibc != bc_->inletBCmap.end()) {
ibc->second->computeBdrPrimitiveStateForGradient(nbdrInlet, iUp1, iUp2);
nbdrInlet++;
nbdrInlet++;
}
if (obc != bc_->outletBCmap.end()) {
obc->second->computeBdrPrimitiveStateForGradient(nbdrOutlet,iUp1, iUp2);
nbdrOutlet++;
obc->second->computeBdrPrimitiveStateForGradient(nbdrOutlet, iUp1, iUp2);
nbdrOutlet++;
}
if (wbc != bc_->wallBCmap.end()) {
wbc->second->computeBdrPrimitiveStateForGradient(nbdrWall,iUp1, iUp2);
nbdrWall++;
wbc->second->computeBdrPrimitiveStateForGradient(nbdrWall, iUp1, iUp2);
nbdrWall++;
}
} else {
iUp2 = iUp1;
Expand Down
14 changes: 7 additions & 7 deletions src/faceGradientIntegration.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -33,9 +33,10 @@
#define FACEGRADIENTINTEGRATION_HPP_

#include <tps_config.h>

#include "BCintegrator.hpp"
#include "tps_mfem_wrap.hpp"
#include "gradNonLinearForm.hpp"
#include "tps_mfem_wrap.hpp"

using namespace mfem;

Expand All @@ -49,19 +50,18 @@ class GradFaceIntegrator : public NonlinearFormIntegrator {
BCintegrator *bc_; // NB: GradFaceIntegrator is a friend of BCintegrator
const bool useBCinGrad_;

int nbdrInlet, nbdrOutlet, nbdrWall;
int nbdrInlet, nbdrOutlet, nbdrWall;

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

void initNbrInlet() {nbdrInlet = 0;}
void initNbrOutlet() {nbdrOutlet = 0;}
void initNbrWall() {nbdrWall = 0;}

void initNbrInlet() { nbdrInlet = 0; }
void initNbrOutlet() { nbdrOutlet = 0; }
void initNbrWall() { nbdrWall = 0; }
};

#endif // FACEGRADIENTINTEGRATION_HPP_
224 changes: 3 additions & 221 deletions src/gradNonLinearForm.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -35,8 +35,8 @@
//#include "doftrans.hpp"

GradNonLinearForm::GradNonLinearForm(ParFiniteElementSpace *_vfes, IntegrationRules *_intRules, const int _dim,
const int _num_equation)
: ParNonlinearForm(_vfes), vfes(_vfes), intRules(_intRules), dim(_dim), num_equation(_num_equation) {}
const int _num_equation, BCintegrator *bc)
: ParNonlinearForm(_vfes), vfes(_vfes), intRules(_intRules), dim(_dim), num_equation(_num_equation), bc_(bc) {}

void GradNonLinearForm::Mult(const ParGridFunction *Up, Vector &y) {
Vector x;
Expand All @@ -53,225 +53,7 @@ void GradNonLinearForm::Mult(const ParGridFunction *Up, Vector &y) {
gfi->initNbrInlet();
gfi->initNbrOutlet();
gfi->initNbrWall();

// After setting x as above, base class Mult does the right thing
ParNonlinearForm::Mult(x, y);

/*
// need to take over base class Mult for correct calcualtion
// of DG gradients. Modified from fem/pnonlinearform.cpp:Mult
// ParNonlinearForm::Mult(x, y) calls NonlinearForm::Mult(x, y),
// so both are here
// NonlinearForm::Mult...................................................
const Vector &px = Prolongate(x);
if (P) { aux2.SetSize(P->Height()); }
// If we are in parallel, ParNonLinearForm::Mult uses the aux2 vector. In
// serial, place the result directly in y (when there is no P).
Vector &py = P ? aux2 : y;
if (ext)
{
ext->Mult(px, py);
if (Serial())
{
if (cP) { cP->MultTranspose(py, y); }
const int N = ess_tdof_list.Size();
const auto tdof = ess_tdof_list.Read();
auto Y = y.ReadWrite();
MFEM_FORALL(i, N, Y[tdof[i]] = 0.0; );
}
// In parallel, the result is in 'py' which is an alias for 'aux2'.
return;
}
Array<int> vdofs;
Vector el_x, el_y;
const FiniteElement *fe;
ElementTransformation *T;
//DofTransformation *doftrans;
Mesh *mesh = fes->GetMesh();
py = 0.0;
if (dnfi.Size())
{
for (int i = 0; i < fes->GetNE(); i++)
{
fe = fes->GetFE(i);
//doftrans = fes->GetElementVDofs(i, vdofs);
T = fes->GetElementTransformation(i);
px.GetSubVector(vdofs, el_x);
//if (doftrans) {doftrans->InvTransformPrimal(el_x); }
for (int k = 0; k < dnfi.Size(); k++)
{
dnfi[k]->AssembleElementVector(*fe, *T, el_x, el_y);
//if (doftrans) {doftrans->TransformDual(el_y); }
py.AddElementVector(vdofs, el_y);
}
}
}
if (fnfi.Size())
{
FaceElementTransformations *tr;
const FiniteElement *fe1, *fe2;
Array<int> vdofs2;
for (int i = 0; i < mesh->GetNumFaces(); i++)
{
tr = mesh->GetInteriorFaceTransformations(i);
if (tr != NULL)
{
fes->GetElementVDofs(tr->Elem1No, vdofs);
fes->GetElementVDofs(tr->Elem2No, vdofs2);
vdofs.Append (vdofs2);
px.GetSubVector(vdofs, el_x);
fe1 = fes->GetFE(tr->Elem1No);
fe2 = fes->GetFE(tr->Elem2No);
for (int k = 0; k < fnfi.Size(); k++)
{
fnfi[k]->AssembleFaceVector(*fe1, *fe2, *tr, el_x, el_y);
py.AddElementVector(vdofs, el_y);
}
}
}
}
if (bfnfi.Size())
{
FaceElementTransformations *tr;
const FiniteElement *fe1, *fe2;
// Which boundary attributes need to be processed?
Array<int> bdr_attr_marker(mesh->bdr_attributes.Size() ?
mesh->bdr_attributes.Max() : 0);
bdr_attr_marker = 0;
for (int k = 0; k < bfnfi.Size(); k++)
{
if (bfnfi_marker[k] == NULL)
{
bdr_attr_marker = 1;
break;
}
Array<int> &bdr_marker = *bfnfi_marker[k];
MFEM_ASSERT(bdr_marker.Size() == bdr_attr_marker.Size(),
"invalid boundary marker for boundary face integrator #"
<< k << ", counting from zero");
for (int i = 0; i < bdr_attr_marker.Size(); i++)
{
bdr_attr_marker[i] |= bdr_marker[i];
}
}
// incrementers for each boundary type
//Nbdr_inlet = 0;
//Nbdr_outlet = 0;
//Nbdr_wall = 0;
//setNumInlet(0);
//setNumOutlet(0);
//setNumWall(0);
//int Nbdr_inlet = 0;
//int Nbdr_outlet = 0;
//int Nbdr_wall = 0;
//bfnfi.initNbrInlet();
//bfnfi.initNbrOutlet();
//bfnfi.initNbrWall();
// physical boundary loop <warp>
for (int i = 0; i < fes -> GetNBE(); i++)
{
const int bdr_attr = mesh->GetBdrAttribute(i);
if (bdr_attr_marker[bdr_attr-1] == 0) { continue; }
tr = mesh->GetBdrFaceTransformations (i);
if (tr != NULL)
{
fes->GetElementVDofs(tr->Elem1No, vdofs);
px.GetSubVector(vdofs, el_x);
fe1 = fes->GetFE(tr->Elem1No);
// The fe2 object is really a dummy and not used on the boundaries,
// but we can't dereference a NULL pointer, and we don't want to
// actually make a fake element.
fe2 = fe1;
for (int k = 0; k < bfnfi.Size(); k++) // remove loop, add assert that size is 1
{
if (bfnfi_marker[k] &&
(*bfnfi_marker[k])[bdr_attr-1] == 0) { continue; }
// IntegrationRule loop is in AsembleFaceVector
bfnfi[k]->AssembleFaceVector(*fe1, *fe2, *tr, el_x, el_y);
py.AddElementVector(vdofs, el_y);
}
}
}
}
if (Serial())
{
if (cP) { cP->MultTranspose(py, y); }
for (int i = 0; i < ess_tdof_list.Size(); i++)
{
y(ess_tdof_list[i]) = 0.0;
}
// y(ess_tdof_list[i]) = x(ess_tdof_list[i]);
}
// In parallel, the result is in 'py' which is an alias for 'aux2'.
// ParNonlinearForm::Mult...................................................
if (fnfi.Size())
{
MFEM_VERIFY(!NonlinearForm::ext, "Not implemented (extensions + faces");
// 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;
aux1.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());
for (int k = 0; k < fnfi.Size(); k++)
{
fnfi[k]->AssembleFaceVector(*fe1, *fe2, *tr, el_x, el_y);
aux2.AddElementVector(vdofs1, el_y.GetData());
}
}
}
P->MultTranspose(aux2, y);
const int N = ess_tdof_list.Size();
const auto idx = ess_tdof_list.Read();
auto Y_RW = y.ReadWrite();
MFEM_FORALL(i, N, Y_RW[idx[i]] = 0.0; );
*/

}
Loading

0 comments on commit 255899e

Please sign in to comment.