Skip to content

Dev v2f improvements #336

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
shaering wants to merge 14 commits into
base: main
Choose a base branch
from
Open

Dev v2f improvements #336

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
14 commits
Select commit Hold shift + click to select a range
7deda8f
adding spark, wip
shaering May 15, 2026
f741b34
added sparking to reacting flow. builds, not tested
shaering May 15, 2026
3f9e4b7
fix to spark coordinates
shaering May 23, 2026
1a19236
Series of improvements to v2f and other small changes: 1) improved sm…
shaering Jun 21, 2026
cafa9ad
merging in main
shaering Jun 21, 2026
da78c20
style
shaering Jun 21, 2026
7ef10b0
style
shaering Jun 21, 2026
eb37a4e
style
shaering Jun 21, 2026
d6fb13b
style
shaering Jun 21, 2026
f20d741
adjusted turbulent pipe test for zetaModel changes
Jun 23, 2026
e017d24
version seems to work with higher-order v2f and give decent results. …
shaering Jun 26, 2026
e57008b
adding robust form of f to allow higher-order v2f runs. fix to inlet…
shaering Jul 1, 2026
204e9bb
fix to joule heating scaling to prevent heating outside of torch
shaering Jul 1, 2026
4112b53
git being strange...
shaering Jul 1, 2026
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
68 changes: 34 additions & 34 deletions src/BCintegrator.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -37,13 +37,13 @@
#include "outletBC.hpp"
#include "wallBC.hpp"

BCintegrator::BCintegrator(bool _mpiRoot, MPI_Groups *_groupsMPI, ParMesh *_mesh, ParFiniteElementSpace *_vfes,
IntegrationRules *_intRules, RiemannSolverTPS *rsolver_, double &_dt, double *_time,
GasMixture *_mixture, GasMixture *d_mixture, Fluxes *_fluxClass, ParGridFunction *_Up,
ParGridFunction *_gradUp, const boundaryFaceIntegrationData &boundary_face_data,
const int _dim, const int _num_equation, double &_max_char_speed, RunConfiguration &_runFile,
Array<int> &local_attr, const int &_maxIntPoints, const int &_maxDofs,
ParGridFunction *distance)
BCintegrator::BCintegrator(bool _mpiRoot, MPI_Groups* _groupsMPI, ParMesh* _mesh, ParFiniteElementSpace* _vfes,
IntegrationRules* _intRules, RiemannSolverTPS* rsolver_, double& _dt, double* _time,
GasMixture* _mixture, GasMixture* d_mixture, Fluxes* _fluxClass, ParGridFunction* _Up,
ParGridFunction* _gradUp, const boundaryFaceIntegrationData& boundary_face_data,
const int _dim, const int _num_equation, double& _max_char_speed, RunConfiguration& _runFile,
Array<int>& local_attr, const int& _maxIntPoints, const int& _maxDofs,
ParGridFunction* distance)
: groupsMPI(_groupsMPI),
config(_runFile),
rsolver(rsolver_),
Expand Down Expand Up @@ -129,8 +129,8 @@ BCintegrator::BCintegrator(bool _mpiRoot, MPI_Groups *_groupsMPI, ParMesh *_mesh

// Inlets
if (NumBCelems > 0 && inletBCmap.size() > 0) {
Mesh *mesh_bc = vfes->GetMesh();
FaceElementTransformations *tr;
Mesh* mesh_bc = vfes->GetMesh();
FaceElementTransformations* tr;

for (int i = 0; i < local_attr.Size(); i++) {
int attr = local_attr[i];
Expand All @@ -145,15 +145,15 @@ BCintegrator::BCintegrator(bool _mpiRoot, MPI_Groups *_groupsMPI, ParMesh *_mesh
}
}

std::unordered_map<int, BoundaryCondition *>::const_iterator ibc = inletBCmap.find(attr);
std::unordered_map<int, BoundaryCondition*>::const_iterator ibc = inletBCmap.find(attr);
if (ibc != inletBCmap.end()) ibc->second->setElementList(list);
}
}

// Outlets
if (NumBCelems > 0 && outletBCmap.size() > 0) {
Mesh *mesh_bc = vfes->GetMesh();
FaceElementTransformations *tr;
Mesh* mesh_bc = vfes->GetMesh();
FaceElementTransformations* tr;

for (int i = 0; i < local_attr.Size(); i++) {
int attr = local_attr[i];
Expand All @@ -168,15 +168,15 @@ BCintegrator::BCintegrator(bool _mpiRoot, MPI_Groups *_groupsMPI, ParMesh *_mesh
}
}

std::unordered_map<int, BoundaryCondition *>::const_iterator obc = outletBCmap.find(attr);
std::unordered_map<int, BoundaryCondition*>::const_iterator obc = outletBCmap.find(attr);
if (obc != outletBCmap.end()) obc->second->setElementList(list);
}
}

// Walls
if (NumBCelems > 0 && wallBCmap.size() > 0) {
Mesh *mesh_bc = vfes->GetMesh();
FaceElementTransformations *tr;
Mesh* mesh_bc = vfes->GetMesh();
FaceElementTransformations* tr;

for (int i = 0; i < local_attr.Size(); i++) {
int attr = local_attr[i];
Expand All @@ -191,7 +191,7 @@ BCintegrator::BCintegrator(bool _mpiRoot, MPI_Groups *_groupsMPI, ParMesh *_mesh
}
}

std::unordered_map<int, BoundaryCondition *>::const_iterator wbc = wallBCmap.find(attr);
std::unordered_map<int, BoundaryCondition*>::const_iterator wbc = wallBCmap.find(attr);
if (wbc != wallBCmap.end()) wbc->second->setElementList(list);
}
}
Expand Down Expand Up @@ -225,11 +225,11 @@ void BCintegrator::initBCs() {
}
}

void BCintegrator::computeBdrFlux(const int attr, Vector &normal, Vector &stateIn, DenseMatrix &gradState,
Vector transip, double delta, double time, double distance, Vector &bdrFlux) {
std::unordered_map<int, BoundaryCondition *>::const_iterator ibc = inletBCmap.find(attr);
std::unordered_map<int, BoundaryCondition *>::const_iterator obc = outletBCmap.find(attr);
std::unordered_map<int, BoundaryCondition *>::const_iterator wbc = wallBCmap.find(attr);
void BCintegrator::computeBdrFlux(const int attr, Vector& normal, Vector& stateIn, DenseMatrix& gradState,
Vector transip, double delta, double time, double distance, Vector& bdrFlux) {
std::unordered_map<int, BoundaryCondition*>::const_iterator ibc = inletBCmap.find(attr);
std::unordered_map<int, BoundaryCondition*>::const_iterator obc = outletBCmap.find(attr);
std::unordered_map<int, BoundaryCondition*>::const_iterator wbc = wallBCmap.find(attr);

if (ibc != inletBCmap.end())
ibc->second->computeBdrFlux(normal, stateIn, gradState, transip, delta, time, distance, bdrFlux);
Expand All @@ -241,7 +241,7 @@ void BCintegrator::computeBdrFlux(const int attr, Vector &normal, Vector &stateI
// BCmap[attr]->computeBdrFlux(normal, stateIn, gradState, radius, bdrFlux);
}

void BCintegrator::updateBCMean(ParGridFunction *Up) {
void BCintegrator::updateBCMean(ParGridFunction* Up) {
for (auto bc = inletBCmap.begin(); bc != inletBCmap.end(); bc++) {
bc->second->updateMean(intRules, Up);
}
Expand All @@ -255,7 +255,7 @@ void BCintegrator::updateBCMean(ParGridFunction *Up) {
}
}

void BCintegrator::integrateBCs(Vector &y, const Vector &x, const elementIndexingData &elem_index_data) {
void BCintegrator::integrateBCs(Vector& y, const Vector& x, const elementIndexingData& elem_index_data) {
for (auto bc = inletBCmap.begin(); bc != inletBCmap.end(); bc++) {
bc->second->integrationBC(y, // output
x, elem_index_data, Up, gradUp, boundary_face_data_, maxIntPoints, maxDofs);
Expand All @@ -272,12 +272,12 @@ void BCintegrator::integrateBCs(Vector &y, const Vector &x, const elementIndexin
}
}

void BCintegrator::retrieveGradientsData_gpu(ParGridFunction *gradUp, DenseTensor &elGradUp, Array<int> &vdofs,
const int &num_equation, const int &dim, const int &totalDofs,
const int &elDofs) {
void BCintegrator::retrieveGradientsData_gpu(ParGridFunction* gradUp, DenseTensor& elGradUp, Array<int>& vdofs,
const int& num_equation, const int& dim, const int& totalDofs,
const int& elDofs) {
#ifdef _GPU_
const double *d_GradUp = gradUp->Read();
double *d_elGradUp = elGradUp.ReadWrite();
const double* d_GradUp = gradUp->Read();
double* d_elGradUp = elGradUp.ReadWrite();
auto d_vdofs = vdofs.Read();

MFEM_FORALL(i, elDofs, {
Expand All @@ -292,15 +292,15 @@ void BCintegrator::retrieveGradientsData_gpu(ParGridFunction *gradUp, DenseTenso
#endif
}

void BCintegrator::AssembleFaceVector(const FiniteElement &el1, const FiniteElement &el2,
FaceElementTransformations &Tr, const Vector &elfun, Vector &elvect) {
void BCintegrator::AssembleFaceVector(const FiniteElement& el1, const FiniteElement& el2,
FaceElementTransformations& Tr, const Vector& elfun, Vector& elvect) {
Vector shape1;
Vector funval1(num_equation);
Vector nor(dim);
Vector fluxN(num_equation);

#ifndef _GPU_
const double *dataGradUp = gradUp->HostRead();
const double* dataGradUp = gradUp->HostRead();
#endif

const int dof1 = el1.GetDof();
Expand All @@ -323,7 +323,7 @@ void BCintegrator::AssembleFaceVector(const FiniteElement &el1, const FiniteElem

// element size
double delta;
Mesh *mesh = vfes->GetMesh();
Mesh* mesh = vfes->GetMesh();
delta = mesh->GetElementSize(Tr.Elem1No, 1) / el1.GetOrder();

#ifdef _GPU_
Expand Down Expand Up @@ -361,12 +361,12 @@ void BCintegrator::AssembleFaceVector(const FiniteElement &el1, const FiniteElem
intorder++;
}

const IntegrationRule *ir = &intRules->Get(Tr.GetGeometryType(), intorder);
const IntegrationRule* ir = &intRules->Get(Tr.GetGeometryType(), intorder);
const int numActiveSpecies = mixture->GetNumActiveSpecies();
const int nvel = mixture->GetNumVels();

for (int i = 0; i < ir->GetNPoints(); i++) {
const IntegrationPoint &ip = ir->IntPoint(i);
const IntegrationPoint& ip = ir->IntPoint(i);

Tr.SetAllIntPoints(&ip); // set face and element int. points

Expand Down
74 changes: 37 additions & 37 deletions src/BCintegrator.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -51,74 +51,74 @@ class BCintegrator : public NonlinearFormIntegrator {
friend class M2ulPhyS;

protected:
MPI_Groups *groupsMPI;
MPI_Groups* groupsMPI;

RunConfiguration &config;
RunConfiguration& config;

RiemannSolverTPS *rsolver;
GasMixture *mixture;
Fluxes *fluxClass;
RiemannSolverTPS* rsolver;
GasMixture* mixture;
Fluxes* fluxClass;

double &max_char_speed;
IntegrationRules *intRules;
double& max_char_speed;
IntegrationRules* intRules;

ParMesh *mesh;
ParMesh* mesh;

// pointer to finite element space
ParFiniteElementSpace *vfes;
ParFiniteElementSpace* vfes;

// pointer to primitive varibales
ParGridFunction *Up;
ParGridFunction* Up;

ParGridFunction *gradUp;
ParGridFunction* gradUp;

ParGridFunction *distance_;
ParGridFunction* distance_;

const boundaryFaceIntegrationData &boundary_face_data_;
const boundaryFaceIntegrationData& boundary_face_data_;

const int dim;
const int num_equation;

const int &maxIntPoints;
const int &maxDofs;
const int& maxIntPoints;
const int& maxDofs;

std::unordered_map<int, BoundaryCondition *> inletBCmap;
std::unordered_map<int, BoundaryCondition *> outletBCmap;
std::unordered_map<int, BoundaryCondition *> wallBCmap;
std::unordered_map<int, BoundaryCondition*> inletBCmap;
std::unordered_map<int, BoundaryCondition*> outletBCmap;
std::unordered_map<int, BoundaryCondition*> wallBCmap;

bool mpiRoot;
double time;
double *pTime;
double* pTime;

// void calcMeanState();
void computeBdrFlux(const int attr, Vector &normal, Vector &stateIn, DenseMatrix &gradState, Vector transip,
double delta, double time, double distance, Vector &bdrFlux);
void computeBdrFlux(const int attr, Vector& normal, Vector& stateIn, DenseMatrix& gradState, Vector transip,
double delta, double time, double distance, Vector& bdrFlux);

public:
BCintegrator(bool _mpiRoot, MPI_Groups *_groupsMPI, ParMesh *_mesh, ParFiniteElementSpace *_vfes,
IntegrationRules *_intRules, RiemannSolverTPS *rsolver_, double &_dt, double *_time, GasMixture *mixture,
GasMixture *d_mixture, Fluxes *_fluxClass, ParGridFunction *_Up, ParGridFunction *_gradUp,
const boundaryFaceIntegrationData &boundary_face_data, const int _dim, const int _num_equation,
double &_max_char_speed, RunConfiguration &_runFile, Array<int> &local_bdr_attr,
const int &_maxIntPoints, const int &_maxDofs, ParGridFunction *distance_);
BCintegrator(bool _mpiRoot, MPI_Groups* _groupsMPI, ParMesh* _mesh, ParFiniteElementSpace* _vfes,
IntegrationRules* _intRules, RiemannSolverTPS* rsolver_, double& _dt, double* _time, GasMixture* mixture,
GasMixture* d_mixture, Fluxes* _fluxClass, ParGridFunction* _Up, ParGridFunction* _gradUp,
const boundaryFaceIntegrationData& boundary_face_data, const int _dim, const int _num_equation,
double& _max_char_speed, RunConfiguration& _runFile, Array<int>& local_bdr_attr,
const int& _maxIntPoints, const int& _maxDofs, ParGridFunction* distance_);
~BCintegrator();

virtual void AssembleFaceVector(const FiniteElement &el1, const FiniteElement &el2, FaceElementTransformations &Tr,
const Vector &elfun, Vector &elvect);
virtual void AssembleFaceVector(const FiniteElement& el1, const FiniteElement& el2, FaceElementTransformations& Tr,
const Vector& elfun, Vector& elvect);
void initBCs();

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

// GPU functions
static void retrieveGradientsData_gpu(ParGridFunction *gradUp, DenseTensor &elGradUp, Array<int> &vdofs,
const int &num_equation, const int &dim, const int &totalDofs,
const int &elDofs);
static void retrieveGradientsData_gpu(ParGridFunction* gradUp, DenseTensor& elGradUp, Array<int>& vdofs,
const int& num_equation, const int& dim, const int& totalDofs,
const int& elDofs);

boundaryCategory getAttributeCategory(int attr) const {
std::unordered_map<int, BoundaryCondition *>::const_iterator ibc = inletBCmap.find(attr);
std::unordered_map<int, BoundaryCondition *>::const_iterator obc = outletBCmap.find(attr);
std::unordered_map<int, BoundaryCondition *>::const_iterator wbc = wallBCmap.find(attr);
std::unordered_map<int, BoundaryCondition*>::const_iterator ibc = inletBCmap.find(attr);
std::unordered_map<int, BoundaryCondition*>::const_iterator obc = outletBCmap.find(attr);
std::unordered_map<int, BoundaryCondition*>::const_iterator wbc = wallBCmap.find(attr);
if (ibc != inletBCmap.end()) {
return INLET;
}
Expand Down
18 changes: 9 additions & 9 deletions src/BoundaryCondition.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -31,8 +31,8 @@
// -----------------------------------------------------------------------------------el-
#include "BoundaryCondition.hpp"

BoundaryCondition::BoundaryCondition(RiemannSolverTPS *_rsolver, GasMixture *_mixture, Equations _eqSystem,
ParFiniteElementSpace *_vfes, IntegrationRules *_intRules, double &_dt,
BoundaryCondition::BoundaryCondition(RiemannSolverTPS* _rsolver, GasMixture* _mixture, Equations _eqSystem,
ParFiniteElementSpace* _vfes, IntegrationRules* _intRules, double& _dt,
const int _dim, const int _num_equation, const int _patchNumber,
const double _refLength, bool axisym)
: rsolver(_rsolver),
Expand All @@ -52,7 +52,7 @@ BoundaryCondition::BoundaryCondition(RiemannSolverTPS *_rsolver, GasMixture *_mi

BoundaryCondition::~BoundaryCondition() {}

void BoundaryCondition::computeBdrPrimitiveStateForGradient(const Vector &stateIn, Vector &stateBC) const {
void BoundaryCondition::computeBdrPrimitiveStateForGradient(const Vector& stateIn, Vector& stateBC) const {
stateBC = stateIn;
}

Expand All @@ -63,8 +63,8 @@ double BoundaryCondition::aggregateArea(int bndry_patchnum, MPI_Comm bc_comm) {
for (int bel = 0; bel < vfes->GetNBE(); bel++) {
int attr = vfes->GetBdrAttribute(bel);
if (attr == bndry_patchnum) {
FaceElementTransformations *Tr = vfes->GetMesh()->GetBdrFaceTransformations(bel);
const IntegrationRule &ir = IntRules.Get(Tr->GetGeometryType(), Tr->OrderJ());
FaceElementTransformations* Tr = vfes->GetMesh()->GetBdrFaceTransformations(bel);
const IntegrationRule& ir = IntRules.Get(Tr->GetGeometryType(), Tr->OrderJ());

for (int p = 0; p < ir.GetNPoints(); p++) {
const IntegrationPoint ip = ir.IntPoint(p);
Expand All @@ -91,15 +91,15 @@ int BoundaryCondition::aggregateBndryFaces(int bndry_patchnum, MPI_Comm bc_comm)
return (nfacesTotal);
}

void BoundaryCondition::setElementList(Array<int> &_listElems) {
void BoundaryCondition::setElementList(Array<int>& _listElems) {
listElems.SetSize(_listElems.Size());
for (int i = 0; i < listElems.Size(); i++) listElems[i] = _listElems[i];
listElems.ReadWrite();
}

void BoundaryCondition::copyValues(const Vector &orig, Vector &target, const double &mult) {
const double *dOrig = orig.Read();
double *dTarget = target.Write();
void BoundaryCondition::copyValues(const Vector& orig, Vector& target, const double& mult) {
const double* dOrig = orig.Read();
double* dTarget = target.Write();

MFEM_FORALL(i, target.Size(), { dTarget[i] = dOrig[i] * mult; });
}
Loading
Loading