Skip to content

Commit

Permalink
compressible bc with face-based uvw (#284)
Browse files Browse the repository at this point in the history 
* adding dirichlet inlet bc with velocity specified relative to inlet patch coordinate system

* added 3 options for second tangent to be aligned with x, y, or z global axis.  Other added files are due to style being not correct

* style check

* some tests use old labels for visc sponge, fixing

* continued

---------

Co-authored-by: Sigfried Haering <shaering@oden.utexas.edu>
  • Loading branch information
@shaering
shaering and Sigfried Haering authored Jun 25, 2024
1 parent c1a8af5 commit f87bef1
Show file tree
Hide file tree
Showing 14 changed files with 227 additions and 47 deletions.
24 changes: 14 additions & 10 deletions src/BCintegrator.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -37,11 +37,11 @@
#include "outletBC.hpp"
#include "wallBC.hpp"

BCintegrator::BCintegrator(MPI_Groups *_groupsMPI, ParMesh *_mesh, ParFiniteElementSpace *_vfes,
IntegrationRules *_intRules, RiemannSolver *rsolver_, double &_dt, 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,
BCintegrator::BCintegrator(bool _mpiRoot, MPI_Groups *_groupsMPI, ParMesh *_mesh, ParFiniteElementSpace *_vfes,
IntegrationRules *_intRules, RiemannSolver *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),
Expand All @@ -65,6 +65,10 @@ BCintegrator::BCintegrator(MPI_Groups *_groupsMPI, ParMesh *_mesh, ParFiniteElem
outletBCmap.clear();
wallBCmap.clear();

mpiRoot = _mpiRoot;
time = *_time;
pTime = _time;

// Init inlet BCs
for (size_t in = 0; in < config.GetInletPatchType()->size(); in++) {
std::pair<int, InletType> patchANDtype = (*config.GetInletPatchType())[in];
Expand Down Expand Up @@ -222,17 +226,17 @@ void BCintegrator::initBCs() {
}

void BCintegrator::computeBdrFlux(const int attr, Vector &normal, Vector &stateIn, DenseMatrix &gradState,
Vector transip, double delta, double distance, Vector &bdrFlux) {
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, distance, bdrFlux);
ibc->second->computeBdrFlux(normal, stateIn, gradState, transip, delta, time, distance, bdrFlux);
if (obc != outletBCmap.end())
obc->second->computeBdrFlux(normal, stateIn, gradState, transip, delta, distance, bdrFlux);
obc->second->computeBdrFlux(normal, stateIn, gradState, transip, delta, time, distance, bdrFlux);
if (wbc != wallBCmap.end())
wbc->second->computeBdrFlux(normal, stateIn, gradState, transip, delta, distance, bdrFlux);
wbc->second->computeBdrFlux(normal, stateIn, gradState, transip, delta, time, distance, bdrFlux);

// BCmap[attr]->computeBdrFlux(normal, stateIn, gradState, radius, bdrFlux);
}
Expand Down Expand Up @@ -417,7 +421,7 @@ void BCintegrator::AssembleFaceVector(const FiniteElement &el1, const FiniteElem
radius = transip[0];
}

computeBdrFlux(Tr.Attribute, nor, funval1, iGradUp, transip, delta, d1, fluxN);
computeBdrFlux(Tr.Attribute, nor, funval1, iGradUp, transip, delta, *pTime, d1, fluxN);
fluxN *= ip.weight;

if (config.isAxisymmetric()) {
Expand Down
17 changes: 11 additions & 6 deletions src/BCintegrator.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -86,16 +86,21 @@ class BCintegrator : public NonlinearFormIntegrator {
std::unordered_map<int, BoundaryCondition *> outletBCmap;
std::unordered_map<int, BoundaryCondition *> wallBCmap;

bool mpiRoot;
double time;
double *pTime;

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

public:
BCintegrator(MPI_Groups *_groupsMPI, ParMesh *_mesh, ParFiniteElementSpace *_vfes, IntegrationRules *_intRules,
RiemannSolver *rsolver_, double &_dt, 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, RiemannSolver *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,
Expand Down
2 changes: 1 addition & 1 deletion src/BoundaryCondition.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -74,7 +74,7 @@ class BoundaryCondition {
virtual ~BoundaryCondition();

virtual void computeBdrFlux(Vector &normal, Vector &stateIn, DenseMatrix &gradState, Vector transip, double delta,
double distance, Vector &bdrFlux) = 0;
double time, double distance, Vector &bdrFlux) = 0;

/** \brief Set the boundary state used in the gradient evaluation
*
Expand Down
43 changes: 30 additions & 13 deletions src/M2ulPhyS.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -685,11 +685,14 @@ void M2ulPhyS::initVariables() {
Array<int> local_attr;
getAttributesInPartition(local_attr);

double *pTime;
pTime = &time;

bcIntegrator = NULL;
if (local_attr.Size() > 0) {
bcIntegrator = new BCintegrator(groupsMPI, mesh, vfes, intRules, rsolver, dt, mixture, d_mixture, d_fluxClass, Up,
gradUp, gpu_precomputed_data_.boundary_face_data, dim, num_equation, max_char_speed,
config, local_attr, maxIntPoints, maxDofs, distance_);
bcIntegrator = new BCintegrator(rank0_, groupsMPI, mesh, vfes, intRules, rsolver, dt, pTime, mixture, d_mixture,
d_fluxClass, Up, gradUp, gpu_precomputed_data_.boundary_face_data, dim,
num_equation, max_char_speed, config, local_attr, maxIntPoints, maxDofs, distance_);
}

// A->SetAssemblyLevel(AssemblyLevel::PARTIAL);
Expand Down Expand Up @@ -2783,19 +2786,19 @@ void M2ulPhyS::parseViscosityOptions() {
tpsP->getInput("viscosityMultiplierFunction/isEnabled", config.linViscData.isEnabled, false);
if (config.linViscData.isEnabled) {
auto normal = config.linViscData.normal.HostWrite();
tpsP->getRequiredVecElem("viscosityMultiplierFunction/norm", normal[0], 0);
tpsP->getRequiredVecElem("viscosityMultiplierFunction/norm", normal[1], 1);
tpsP->getRequiredVecElem("viscosityMultiplierFunction/norm", normal[2], 2);
tpsP->getRequiredVecElem("viscosityMultiplierFunction/normal", normal[0], 0);
tpsP->getRequiredVecElem("viscosityMultiplierFunction/normal", normal[1], 1);
tpsP->getRequiredVecElem("viscosityMultiplierFunction/normal", normal[2], 2);

auto point0 = config.linViscData.point0.HostWrite();
tpsP->getRequiredVecElem("viscosityMultiplierFunction/p0", point0[0], 0);
tpsP->getRequiredVecElem("viscosityMultiplierFunction/p0", point0[1], 1);
tpsP->getRequiredVecElem("viscosityMultiplierFunction/p0", point0[2], 2);
tpsP->getRequiredVecElem("viscosityMultiplierFunction/point", point0[0], 0);
tpsP->getRequiredVecElem("viscosityMultiplierFunction/point", point0[1], 1);
tpsP->getRequiredVecElem("viscosityMultiplierFunction/point", point0[2], 2);

auto pointInit = config.linViscData.pointInit.HostWrite();
tpsP->getRequiredVecElem("viscosityMultiplierFunction/pInit", pointInit[0], 0);
tpsP->getRequiredVecElem("viscosityMultiplierFunction/pInit", pointInit[1], 1);
tpsP->getRequiredVecElem("viscosityMultiplierFunction/pInit", pointInit[2], 2);
// auto pointInit = config.linViscData.pointInit.HostWrite();
// tpsP->getRequiredVecElem("viscosityMultiplierFunction/pInit", pointInit[0], 0);
// tpsP->getRequiredVecElem("viscosityMultiplierFunction/pInit", pointInit[1], 1);
// tpsP->getRequiredVecElem("viscosityMultiplierFunction/pInit", pointInit[2], 2);

tpsP->getRequiredInput("viscosityMultiplierFunction/width", config.linViscData.width);
tpsP->getRequiredInput("viscosityMultiplierFunction/viscosityRatio", config.linViscData.viscRatio);
Expand Down Expand Up @@ -3563,6 +3566,9 @@ void M2ulPhyS::parseBCInputs() {
// Inlet Bcs
std::map<std::string, InletType> inletMapping;
inletMapping["subsonic"] = SUB_DENS_VEL;
inletMapping["subsonicFaceBasedX"] = SUB_DENS_VEL_FACE_X;
inletMapping["subsonicFaceBasedY"] = SUB_DENS_VEL_FACE_Y;
inletMapping["subsonicFaceBasedZ"] = SUB_DENS_VEL_FACE_Z;
inletMapping["nonreflecting"] = SUB_DENS_VEL_NR;
inletMapping["nonreflectingConstEntropy"] = SUB_VEL_CONST_ENT;

Expand All @@ -3582,6 +3588,17 @@ void M2ulPhyS::parseBCInputs() {
config.inletBC.Append(density);
config.inletBC.Append(uvw, 3);
}

// NOTE: if there is a need, make face-based general by specifying tangent, leave
// as either tan = x, y, or z for now as fixing requires changing BC interface
// across the board
// FaceBased requires an additional tangent definition cooresponding to the w-component in uvw
// if (type == "subsonicFaceBased") {
// Array<double> tangent2;
// tpsP->getRequiredVec((basepath + "/w-tangent").c_str(), tangent2, 3);
// config.inletBC.Append(tangent2, 3);
//}

// For multi-component gas, require (numActiveSpecies)-more inputs.
if ((config.workFluid != DRY_AIR) && (config.numSpecies > 1)) {
grvy_printf(GRVY_INFO, "\nInlet mass fraction of background species will not be used. \n");
Expand Down
16 changes: 11 additions & 5 deletions src/dataStructures.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -128,11 +128,17 @@ enum SpeciesPrimitiveType { MASS_FRACTION, MOLE_FRACTION, NUMBER_DENSITY, NUM_SP
enum boundaryCategory { INLET, OUTLET, WALL, NUM_BC_CATEGORIES };

enum InletType {
UNI_DENS_VEL, // uniform inlet with density and vel specified
INTERPOLATE, // from an external data file
SUB_DENS_VEL, // Subsonic inlet specified by the density and velocity components
SUB_DENS_VEL_NR, // Non-reflecting subsonic inlet specified by the density and velocity components
SUB_VEL_CONST_ENT // Subsonic non-reflecting. Specified vel, keeps entropy constant
UNI_DENS_VEL, // uniform inlet with density and vel specified
INTERPOLATE, // from an external data file
SUB_DENS_VEL, // Subsonic inlet specified by the density and velocity components
SUB_DENS_VEL_FACE_X, // Subsonic inlet specified by the density and velocity component relative to the inlet face
// with u = normal, w = global x
SUB_DENS_VEL_FACE_Y, // Subsonic inlet specified by the density and velocity component relative to the inlet face
// with u = normal, w = global y
SUB_DENS_VEL_FACE_Z, // Subsonic inlet specified by the density and velocity component relative to the inlet face
// with u = normal, w = global z
SUB_DENS_VEL_NR, // Non-reflecting subsonic inlet specified by the density and velocity components
SUB_VEL_CONST_ENT // Subsonic non-reflecting. Specified vel, keeps entropy constant
};

enum OutletType {
Expand Down
146 changes: 145 additions & 1 deletion src/inletBC.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -139,6 +139,11 @@ InletBC::InletBC(MPI_Groups *_groupsMPI, Equations _eqSystem, RiemannSolver *_rs
boundaryU.UseDevice(true);
boundaryU.SetSize(bdrN * num_equation_);
boundaryU = 0.;

boundaryUp.UseDevice(true);
boundaryUp.SetSize(bdrN * num_equation_);
boundaryUp = 0.;

Vector iState, iUp;
iState.UseDevice(false);
iUp.UseDevice(false);
Expand Down Expand Up @@ -210,6 +215,7 @@ InletBC::InletBC(MPI_Groups *_groupsMPI, Equations _eqSystem, RiemannSolver *_rs

meanUp.Read();
boundaryU.Read();
boundaryUp.Read();
tangent1.Read();
tangent2.Read();
}
Expand Down Expand Up @@ -433,14 +439,29 @@ void InletBC::initBoundaryU(ParGridFunction *Up) {
}

boundaryU.Read();
boundaryUp.Read();
}

void InletBC::computeBdrFlux(Vector &normal, Vector &stateIn, DenseMatrix &gradState, Vector transip, double delta,
double distance, Vector &bdrFlux) {
double time, double distance, Vector &bdrFlux) {
Vector tangentW(dim_);
for (int d = 0; d < dim_; d++) tangentW[d] = 0.0;
switch (inletType_) {
case SUB_DENS_VEL:
subsonicReflectingDensityVelocity(normal, stateIn, bdrFlux);
break;
case SUB_DENS_VEL_FACE_X:
tangentW[0] = 1.0;
subsonicReflectingDensityVelocityFace(normal, tangentW, stateIn, transip, time, bdrFlux);
break;
case SUB_DENS_VEL_FACE_Y:
tangentW[1] = 1.0;
subsonicReflectingDensityVelocityFace(normal, tangentW, stateIn, transip, time, bdrFlux);
break;
case SUB_DENS_VEL_FACE_Z:
tangentW[2] = 1.0;
subsonicReflectingDensityVelocityFace(normal, tangentW, stateIn, transip, time, bdrFlux);
break;
case SUB_DENS_VEL_NR:
subsonicNonReflectingDensityVelocity(normal, stateIn, gradState, bdrFlux);
break;
Expand Down Expand Up @@ -734,6 +755,114 @@ void InletBC::subsonicReflectingDensityVelocity(Vector &normal, Vector &stateIn,
rsolver->Eval(stateIn, state2, normal, bdrFlux, true);
}

/// Specifying subsonic vel and rho wher v is relative to the FACE-coordinate system specifyied by the face normal and
/// tangentW
void InletBC::subsonicReflectingDensityVelocityFace(Vector &normal, Vector tangentW, Vector &stateIn, Vector transip,
double time, Vector &bdrFlux) {
const double p = mixture->ComputePressure(stateIn);

Vector state2(num_equation_);
state2 = stateIn;

// double Rgas = mixture->GetGasConstant();
// double pi = 3.14159265359;
Vector unitNorm;

double tRamp, wt;
tRamp = 0.1; // make this readable
wt = time / tRamp;
wt = min(wt, 1.0);
wt = 1.0;

// injectsion relative to face
double Un = wt * inputState[1];
double Ut = wt * inputState[2];

// unit normals pointing INTO of domain
unitNorm = normal;
double mod = 0.;
for (int d = 0; d < dim_; d++) mod += normal[d] * normal[d];
unitNorm *= -1.0 / sqrt(mod); // inward-facing normal

// t2 aligned with tangent-w
for (int d = 0; d < dim_; d++) tangent2[d] = tangentW[d];

// ensure normal is orthogonal to tangent-w
{
Vector projt2(dim_);
double tmag, tn;
tmag = 0.0;
tn = 0.0;
for (int d = 0; d < dim_; d++) tmag += tangent2[d] * tangent2[d];
for (int d = 0; d < dim_; d++) tn += tangent2[d] * unitNorm[d];
for (int d = 0; d < dim_; d++) projt2[d] = (tn / tmag) * tangent2[d];
for (int d = 0; d < dim_; d++) unitNorm[d] -= projt2[d];
}

// t1 is then orthogonal to both normal and y-axis (clockwise, looking down +)
tangent1[0] = +(unitNorm[1] * tangent2[2] - unitNorm[2] * tangent2[1]);
tangent1[1] = -(unitNorm[0] * tangent2[2] - unitNorm[2] * tangent2[0]);
tangent1[2] = +(unitNorm[0] * tangent2[1] - unitNorm[1] * tangent2[0]);

// aligned with face coords
state2[0] = inputState[0];
// if using temp as the input input, then... p / (Rgas * inputState[0]);
state2[1] = state2[0] * Un;
state2[2] = state2[0] * Ut;
if (nvel_ == 3) state2[3] = state2[0] * 0.0;

if (eqSystem == NS_PASSIVE) {
state2[num_equation_ - 1] = 0.;
} else if (numActiveSpecies_ > 0) {
for (int sp = 0; sp < numActiveSpecies_; sp++) {
// NOTE: inlet BC does not specify total energy. therefore skips one index.
// NOTE: regardless of dim_ension, inletBC save the first 4 elements for density and velocity.
state2[nvel_ + 2 + sp] = inputState[4 + sp];
}
}

// transform from face coords to global
{
DenseMatrix M(dim_, dim_);
for (int d = 0; d < dim_; d++) {
M(0, d) = unitNorm[d];
M(1, d) = tangent1[d];
if (dim_ == 3) M(2, d) = tangent2[d];
}
DenseMatrix invM(dim_, dim_);
mfem::CalcInverse(M, invM);
Vector momN(dim_), momX(dim_);
for (int d = 0; d < dim_; d++) momN[d] = state2[1 + d];
invM.Mult(momN, momX);
for (int d = 0; d < dim_; d++) state2[1 + d] = momX[d];
}

if (eqSystem == NS_PASSIVE) {
state2[num_equation_ - 1] = 0.;
} else if (numActiveSpecies_ > 0) {
for (int sp = 0; sp < numActiveSpecies_; sp++) {
// NOTE: inlet BC does not specify total energy. therefore skips one index.
// NOTE: regardless of dim_ension, inletBC save the first 4 elements for density and velocity.
state2[nvel_ + 2 + sp] = inputState[4 + sp];
}
}

// NOTE: If two-temperature, BC for electron temperature is T_e = T_h, where the total pressure is p.
Vector tmpU(num_equation_);
tmpU = state2;
for (int eq = 1; eq <= dim_; eq++) tmpU[eq] = 2.0 * state2[eq] - stateIn[eq];
mixture->modifyEnergyForPressure(tmpU, tmpU, p, true);
rsolver->Eval(stateIn, tmpU, normal, bdrFlux, true);

// store primitive in boundaryU for gradient calcs
Vector iUp(num_equation_);
mixture->GetPrimitivesFromConservatives(tmpU, iUp);
for (int eq = 0; eq < num_equation_; eq++) {
boundaryUp[eq + bdrN * num_equation_] = iUp[eq];
}
bdrN++;
}

void InletBC::integrateInlets_gpu(Vector &y, const Vector &x, const elementIndexingData &elem_index_data,
const boundaryFaceIntegrationData &boundary_face_data, Array<int> &listElems,
Array<int> &offsetsBoundaryU) {
Expand Down Expand Up @@ -887,6 +1016,21 @@ void InletBC::interpInlet_gpu(const mfem::Vector &x, const elementIndexingData &
pluginInputState(d_inputState, u2, nvel, numActiveSpecies);
d_mix->modifyEnergyForPressure(u2, u2, p, true);
break;
case InletType::SUB_DENS_VEL_FACE_X:
p = d_mix->ComputePressure(u1);
pluginInputState(d_inputState, u2, nvel, numActiveSpecies);
d_mix->modifyEnergyForPressure(u2, u2, p, true);
break;
case InletType::SUB_DENS_VEL_FACE_Y:
p = d_mix->ComputePressure(u1);
pluginInputState(d_inputState, u2, nvel, numActiveSpecies);
d_mix->modifyEnergyForPressure(u2, u2, p, true);
break;
case InletType::SUB_DENS_VEL_FACE_Z:
p = d_mix->ComputePressure(u1);
pluginInputState(d_inputState, u2, nvel, numActiveSpecies);
d_mix->modifyEnergyForPressure(u2, u2, p, true);
break;
case InletType::SUB_DENS_VEL_NR:
printf("INLET BC NOT IMPLEMENTED");
// exit(1);
Expand Down
Loading

0 comments on commit f87bef1

Please sign in to comment.