Commit c8718cf4 authored by Gianmarco Mengaldo's avatar Gianmarco Mengaldo
Browse files

Merge branch 'fix/cfs-dgopt' into 'master'

Fix boundary conditions after DG optimisation

This MR fixes issues from the DG optimisation introduced in MR !555 inside the compressible flow solver, whereby MPI deadlocks under the use of most boundary conditions due to the use of `ExtractTracePhys` inside the boundary condition evaluation. This is fixed by pulling these calls out to the parent function and passing them through to the boundary condition functions, which is in fact more efficient and avoids the use of these functions repeatedly for each RHS evaluation.

See merge request !657
parents 480c850d 0f998d37
......@@ -32,6 +32,7 @@ v4.3.3
**CompressibleFlowSolver**:
- Fix issue with residual output (!647)
- Issues with 1D Euler solver fixed (!565)
- Fix deadlocking issue with boundary conditions (!657)
**Packaging**:
- Fix NekMesh dependencies for DEB package (!650)
......
......@@ -249,11 +249,20 @@ namespace Nektar
fields[i]->GetFwdBwdTracePhys(inarray[i], Fwd[i], Bwd[i]);
fields[0]->GetTrace()->Upwind(Vn, Fwd[i], Bwd[i], numflux[i]);
}
// Extract internal values of the scalar variables for Neumann bcs
Array< OneD, Array<OneD, NekDouble > > uplus(nScalars);
for (i = 0; i < nScalars; ++i)
{
uplus[i] = Array<OneD, NekDouble>(nTracePts, 0.0);
fields[i]->ExtractTracePhys(inarray[i], uplus[i]);
}
// Modify the values in case of boundary interfaces
if (fields[0]->GetBndCondExpansions().num_elements())
{
v_WeakPenaltyO1(fields, inarray, numflux);
v_WeakPenaltyO1(fields, inarray, uplus, numflux);
}
// Splitting the numerical flux into the dimensions
......@@ -273,7 +282,8 @@ namespace Nektar
*/
void DiffusionLDGNS::v_WeakPenaltyO1(
const Array<OneD, MultiRegions::ExpListSharedPtr> &fields,
const Array<OneD, Array<OneD, NekDouble> > &inarray,
const Array<OneD, Array<OneD, NekDouble> > &inarray,
const Array<OneD, Array<OneD, NekDouble> > &uplus,
Array<OneD, Array<OneD, NekDouble> > &penaltyfluxO1)
{
int cnt;
......@@ -290,17 +300,13 @@ namespace Nektar
Array<OneD, NekDouble> Tw(nTracePts, m_Twall);
Array< OneD, Array<OneD, NekDouble > > scalarVariables(nScalars);
Array< OneD, Array<OneD, NekDouble > > uplus(nScalars);
// Extract internal values of the scalar variables for Neumann bcs
for (i = 0; i < nScalars; ++i)
{
scalarVariables[i] = Array<OneD, NekDouble>(nTracePts, 0.0);
uplus[i] = Array<OneD, NekDouble>(nTracePts, 0.0);
fields[i]->ExtractTracePhys(inarray[i], uplus[i]);
}
// Compute boundary conditions for velocities
for (i = 0; i < nScalars-1; ++i)
{
......@@ -502,7 +508,9 @@ namespace Nektar
Vmath::Svtvp(nTracePts, 1.0, m_traceNormals[i], 1,
Vn, 1, Vn, 1);
}
Array<OneD, NekDouble > qtemp(nTracePts);
// Evaulate Riemann flux
// qflux = \hat{q} \cdot u = q \cdot n
// Notice: i = 1 (first row of the viscous tensor is zero)
......@@ -520,11 +528,14 @@ namespace Nektar
// Multiply the Riemann flux by the trace normals
Vmath::Vmul(nTracePts, m_traceNormals[j], 1, qfluxtemp, 1,
qfluxtemp, 1);
// Extract the physical values of the solution at the boundaries
fields[i]->ExtractTracePhys(qfield[j][i], qtemp);
// Impose weak boundary condition with flux
if (fields[0]->GetBndCondExpansions().num_elements())
{
v_WeakPenaltyO2(fields, i, j, qfield[j][i], qfluxtemp);
v_WeakPenaltyO2(fields, i, j, qfield[j][i], qtemp, qfluxtemp);
}
// Store the final flux into qflux
......@@ -544,6 +555,7 @@ namespace Nektar
const int var,
const int dir,
const Array<OneD, const NekDouble> &qfield,
const Array<OneD, const NekDouble> &qtemp,
Array<OneD, NekDouble> &penaltyflux)
{
int cnt = 0;
......@@ -555,11 +567,7 @@ namespace Nektar
int nBndRegions = fields[var]->GetBndCondExpansions().num_elements();
Array<OneD, NekDouble > uterm(nTracePts);
Array<OneD, NekDouble > qtemp(nTracePts);
// Extract the physical values of the solution at the boundaries
fields[var]->ExtractTracePhys(qfield, qtemp);
// Loop on the boundary regions to apply appropriate bcs
for (i = 0; i < nBndRegions; ++i)
{
......
......@@ -93,6 +93,7 @@ namespace Nektar
virtual void v_WeakPenaltyO1(
const Array<OneD, MultiRegions::ExpListSharedPtr> &fields,
const Array<OneD, Array<OneD, NekDouble> > &inarray,
const Array<OneD, Array<OneD, NekDouble> > &uplus,
Array<OneD, Array<OneD, NekDouble> > &penaltyfluxO1);
virtual void v_NumericalFluxO2(
......@@ -106,6 +107,7 @@ namespace Nektar
const int var,
const int dir,
const Array<OneD, const NekDouble> &qfield,
const Array<OneD, const NekDouble> &qtemp,
Array<OneD, NekDouble> &penaltyflux);
virtual void v_SetHomoDerivs(
......
......@@ -326,6 +326,7 @@ namespace Nektar
const int n,
const NekDouble time,
int &cnt,
Array<OneD, Array<OneD, NekDouble> > &Fwd,
Array<OneD, Array<OneD, NekDouble> > &inarray)
{
std::string varName;
......@@ -335,48 +336,48 @@ namespace Nektar
{
if(boost::iequals(userDefStr,"Wall"))
{
WallBC(n, cnt, inarray);
WallBC(n, cnt, Fwd, inarray);
}
else if(boost::iequals(userDefStr,"WallViscous") ||
boost::iequals(userDefStr,"WallAdiabatic"))
{
// Wall Boundary Condition
WallViscousBC(n, cnt, inarray);
WallViscousBC(n, cnt, Fwd, inarray);
}
else if(boost::iequals(userDefStr,"Symmetry"))
{
// Symmetric Boundary Condition
SymmetryBC(n, cnt, inarray);
SymmetryBC(n, cnt, Fwd, inarray);
}
else if(boost::iequals(userDefStr,"RiemannInvariant"))
{
// Riemann invariant characteristic Boundary Condition
RiemannInvariantBC(n, cnt, inarray);
RiemannInvariantBC(n, cnt, Fwd, inarray);
}
else if(boost::iequals(userDefStr,"PressureOutflowNonReflective"))
{
// Pressure outflow non-reflective Boundary Condition
PressureOutflowNonReflectiveBC(n, cnt, inarray);
PressureOutflowNonReflectiveBC(n, cnt, Fwd, inarray);
}
else if(boost::iequals(userDefStr,"PressureOutflow"))
{
// Pressure outflow Boundary Condition
PressureOutflowBC(n, cnt, inarray);
PressureOutflowBC(n, cnt, Fwd, inarray);
}
else if(boost::iequals(userDefStr,"PressureOutflowFile"))
{
// Pressure outflow Boundary Condition from file
PressureOutflowFileBC(n, cnt, inarray);
PressureOutflowFileBC(n, cnt, Fwd, inarray);
}
else if(boost::iequals(userDefStr,"PressureInflowFile"))
{
// Pressure inflow Boundary Condition from file
PressureInflowFileBC(n, cnt, inarray);
PressureInflowFileBC(n, cnt, Fwd, inarray);
}
else if(boost::iequals(userDefStr,"ExtrapOrder0"))
{
// Extrapolation of the data at the boundaries
ExtrapOrder0BC(n, cnt, inarray);
ExtrapOrder0BC(n, cnt, Fwd, inarray);
}
else if(boost::iequals(userDefStr,"TimeDependent"))
{
......@@ -401,23 +402,15 @@ namespace Nektar
void CompressibleFlowSystem::WallBC(
int bcRegion,
int cnt,
Array<OneD, Array<OneD, NekDouble> > &Fwd,
Array<OneD, Array<OneD, NekDouble> > &physarray)
{
int i;
int nTracePts = GetTraceTotPoints();
int nVariables = physarray.num_elements();
const Array<OneD, const int> &traceBndMap
= m_fields[0]->GetTraceBndMap();
// Get physical values of the forward trace
Array<OneD, Array<OneD, NekDouble> > Fwd(nVariables);
for (i = 0; i < nVariables; ++i)
{
Fwd[i] = Array<OneD, NekDouble>(nTracePts);
m_fields[i]->ExtractTracePhys(physarray[i], Fwd[i]);
}
// Adjust the physical values of the trace to take
// user defined boundaries into account
int e, id1, id2, nBCEdgePts, eMax;
......@@ -498,23 +491,15 @@ namespace Nektar
void CompressibleFlowSystem::WallViscousBC(
int bcRegion,
int cnt,
Array<OneD, Array<OneD, NekDouble> > &Fwd,
Array<OneD, Array<OneD, NekDouble> > &physarray)
{
int i;
int nTracePts = GetTraceTotPoints();
int nVariables = physarray.num_elements();
const Array<OneD, const int> &traceBndMap
= m_fields[0]->GetTraceBndMap();
// Get physical values of the forward trace
Array<OneD, Array<OneD, NekDouble> > Fwd(nVariables);
for (i = 0; i < nVariables; ++i)
{
Fwd[i] = Array<OneD, NekDouble>(nTracePts);
m_fields[i]->ExtractTracePhys(physarray[i], Fwd[i]);
}
// Take into account that for PDE based shock capturing, eps = 0 at the
// wall. Adjust the physical values of the trace to take user defined
// boundaries into account
......@@ -573,23 +558,15 @@ namespace Nektar
void CompressibleFlowSystem::SymmetryBC(
int bcRegion,
int cnt,
Array<OneD, Array<OneD, NekDouble> > &Fwd,
Array<OneD, Array<OneD, NekDouble> > &physarray)
{
int i;
int nTracePts = GetTraceTotPoints();
int nVariables = physarray.num_elements();
const Array<OneD, const int> &traceBndMap
= m_fields[0]->GetTraceBndMap();
// Get physical values of the forward trace (from exp to phys)
Array<OneD, Array<OneD, NekDouble> > Fwd(nVariables);
for (i = 0; i < nVariables; ++i)
{
Fwd[i] = Array<OneD, NekDouble>(nTracePts);
m_fields[i]->ExtractTracePhys(physarray[i], Fwd[i]);
}
// Take into account that for PDE based shock capturing, eps = 0 at the
// wall.
int e, id1, id2, nBCEdgePts, eMax;
......@@ -669,11 +646,11 @@ namespace Nektar
void CompressibleFlowSystem::RiemannInvariantBC(
int bcRegion,
int cnt,
Array<OneD, Array<OneD, NekDouble> > &Fwd,
Array<OneD, Array<OneD, NekDouble> > &physarray)
{
int i, j;
int nTracePts = GetTraceTotPoints();
int nVariables = physarray.num_elements();
int nDimensions = m_spacedim;
const Array<OneD, const int> &traceBndMap
......@@ -706,14 +683,6 @@ namespace Nektar
Vmath::Vadd(nTracePts, VnInf, 1, tmp2, 1, VnInf, 1);
}
// Get physical values of the forward trace
Array<OneD, Array<OneD, NekDouble> > Fwd(nVariables);
for (i = 0; i < nVariables; ++i)
{
Fwd[i] = Array<OneD, NekDouble>(nTracePts);
m_fields[i]->ExtractTracePhys(physarray[i], Fwd[i]);
}
// Computing the normal velocity for characteristics coming
// from inside the computational domain
Array<OneD, NekDouble > Vn (nTracePts, 0.0);
......@@ -926,6 +895,7 @@ namespace Nektar
void CompressibleFlowSystem::PressureOutflowNonReflectiveBC(
int bcRegion,
int cnt,
Array<OneD, Array<OneD, NekDouble> > &Fwd,
Array<OneD, Array<OneD, NekDouble> > &physarray)
{
int i, j;
......@@ -962,14 +932,6 @@ namespace Nektar
Vmath::Vadd(nTracePts, VnInf, 1, tmp2, 1, VnInf, 1);
}
// Get physical values of the forward trace
Array<OneD, Array<OneD, NekDouble> > Fwd(nVariables);
for (i = 0; i < nVariables; ++i)
{
Fwd[i] = Array<OneD, NekDouble>(nTracePts);
m_fields[i]->ExtractTracePhys(physarray[i], Fwd[i]);
}
// Computing the normal velocity for characteristics coming
// from inside the computational domain
Array<OneD, NekDouble > Vn (nTracePts, 0.0);
......@@ -1086,6 +1048,7 @@ namespace Nektar
void CompressibleFlowSystem::PressureOutflowBC(
int bcRegion,
int cnt,
Array<OneD, Array<OneD, NekDouble> > &Fwd,
Array<OneD, Array<OneD, NekDouble> > &physarray)
{
int i, j;
......@@ -1122,14 +1085,6 @@ namespace Nektar
Vmath::Vadd(nTracePts, VnInf, 1, tmp2, 1, VnInf, 1);
}
// Get physical values of the forward trace
Array<OneD, Array<OneD, NekDouble> > Fwd(nVariables);
for (i = 0; i < nVariables; ++i)
{
Fwd[i] = Array<OneD, NekDouble>(nTracePts);
m_fields[i]->ExtractTracePhys(physarray[i], Fwd[i]);
}
// Computing the normal velocity for characteristics coming
// from inside the computational domain
Array<OneD, NekDouble > Vn (nTracePts, 0.0);
......@@ -1247,6 +1202,7 @@ namespace Nektar
void CompressibleFlowSystem::PressureOutflowFileBC(
int bcRegion,
int cnt,
Array<OneD, Array<OneD, NekDouble> > &Fwd,
Array<OneD, Array<OneD, NekDouble> > &physarray)
{
int i, j;
......@@ -1283,15 +1239,6 @@ namespace Nektar
Vmath::Vadd(nTracePts, VnInf, 1, tmp2, 1, VnInf, 1);
}
// Get physical values of the forward trace
Array<OneD, Array<OneD, NekDouble> > Fwd(nVariables);
for (i = 0; i < nVariables; ++i)
{
Fwd[i] = Array<OneD, NekDouble>(nTracePts, 0.0);
m_fields[i]->ExtractTracePhys(physarray[i], Fwd[i]);
}
// Computing the normal velocity for characteristics coming
// from inside the computational domain
Array<OneD, NekDouble > Vn (nTracePts, 0.0);
......@@ -1411,6 +1358,7 @@ namespace Nektar
void CompressibleFlowSystem::PressureInflowFileBC(
int bcRegion,
int cnt,
Array<OneD, Array<OneD, NekDouble> > &Fwd,
Array<OneD, Array<OneD, NekDouble> > &physarray)
{
int i, j;
......@@ -1447,15 +1395,6 @@ namespace Nektar
Vmath::Vadd(nTracePts, VnInf, 1, tmp2, 1, VnInf, 1);
}
// Get physical values of the forward trace
Array<OneD, Array<OneD, NekDouble> > Fwd(nVariables);
for (i = 0; i < nVariables; ++i)
{
Fwd[i] = Array<OneD, NekDouble>(nTracePts, 0.0);
m_fields[i]->ExtractTracePhys(physarray[i], Fwd[i]);
}
// Computing the normal velocity for characteristics coming
// from inside the computational domain
Array<OneD, NekDouble > Vn (nTracePts, 0.0);
......@@ -1574,26 +1513,18 @@ namespace Nektar
void CompressibleFlowSystem::ExtrapOrder0BC(
int bcRegion,
int cnt,
Array<OneD, Array<OneD, NekDouble> > &Fwd,
Array<OneD, Array<OneD, NekDouble> > &physarray)
{
int i, j;
int e, pnt;
int id1, id2, nBCEdgePts;
int nTracePts = GetTraceTotPoints();
int nVariables = physarray.num_elements();
int nDimensions = m_spacedim;
const Array<OneD, const int> &traceBndMap
= m_fields[0]->GetTraceBndMap();
// Get physical values of the forward trace
Array<OneD, Array<OneD, NekDouble> > Fwd(nVariables);
for (i = 0; i < nVariables; ++i)
{
Fwd[i] = Array<OneD, NekDouble>(nTracePts);
m_fields[i]->ExtractTracePhys(physarray[i], Fwd[i]);
}
int eMax;
eMax = m_fields[0]->GetBndCondExpansions()[bcRegion]->GetExpSize();
......
......@@ -164,43 +164,52 @@ namespace Nektar
const int n,
const NekDouble time,
int &cnt,
Array<OneD, Array<OneD, NekDouble> > &Fwd,
Array<OneD, Array<OneD, NekDouble> > &inarray);
void WallBC(
int bcRegion,
int cnt,
Array<OneD, Array<OneD, NekDouble> > &Fwd,
Array<OneD, Array<OneD, NekDouble> > &physarray);
void WallViscousBC(
int bcRegion,
int cnt,
Array<OneD, Array<OneD, NekDouble> > &Fwd,
Array<OneD, Array<OneD, NekDouble> > &physarray);
void SymmetryBC(
int bcRegion,
int cnt,
Array<OneD, Array<OneD, NekDouble> > &Fwd,
Array<OneD, Array<OneD, NekDouble> > &physarray);
void RiemannInvariantBC(
int bcRegion,
int cnt,
Array<OneD, Array<OneD, NekDouble> > &Fwd,
Array<OneD, Array<OneD, NekDouble> > &physarray);
void PressureOutflowNonReflectiveBC(
int bcRegion,
int cnt,
Array<OneD, Array<OneD, NekDouble> > &Fwd,
Array<OneD, Array<OneD, NekDouble> > &physarray);
void PressureOutflowBC(
int bcRegion,
int cnt,
Array<OneD, Array<OneD, NekDouble> > &Fwd,
Array<OneD, Array<OneD, NekDouble> > &physarray);
void PressureOutflowFileBC(
int bcRegion,
int cnt,
Array<OneD, Array<OneD, NekDouble> > &Fwd,
Array<OneD, Array<OneD, NekDouble> > &physarray);
void PressureInflowFileBC(
int bcRegion,
int cnt,
Array<OneD, Array<OneD, NekDouble> > &Fwd,
Array<OneD, Array<OneD, NekDouble> > &physarray);
void ExtrapOrder0BC(
int bcRegion,
int cnt,
Array<OneD, Array<OneD, NekDouble> > &Fwd,
Array<OneD, Array<OneD, NekDouble> > &physarray);
void GetVelocityVector(
const Array<OneD, Array<OneD, NekDouble> > &physfield,
......
......@@ -274,6 +274,15 @@ namespace Nektar
{
std::string varName;
int cnt = 0;
int nTracePts = GetTraceTotPoints();
int nvariables = inarray.num_elements();
Array<OneD, Array<OneD, NekDouble> > Fwd(nvariables);
for (int i = 0; i < nvariables; ++i)
{
Fwd[i] = Array<OneD, NekDouble>(nTracePts);
m_fields[i]->ExtractTracePhys(inarray[i], Fwd[i]);
}
// loop over Boundary Regions
for (int n = 0; n < m_fields[0]->GetBndConditions().num_elements(); ++n)
......@@ -289,12 +298,12 @@ namespace Nektar
}
else
{
SetCommonBC(type,n,time, cnt,inarray);
SetCommonBC(type, n, time, cnt, Fwd, inarray);
}
// no User Defined conditions provided so skip cnt
// this line is left in case solver specific condition is added.
cnt += m_fields[0]->GetBndCondExpansions()[n]->GetExpSize();
// no User Defined conditions provided so skip cnt
// this line is left in case solver specific condition is added.
cnt += m_fields[0]->GetBndCondExpansions()[n]->GetExpSize();
}
}
}
......@@ -223,7 +223,16 @@ namespace Nektar
{
std::string varName;
int cnt = 0;
int nvariables = inarray.num_elements();
int nTracePts = GetTraceTotPoints();
Array<OneD, Array<OneD, NekDouble> > Fwd(nvariables);
for (int i = 0; i < nvariables; ++i)
{
Fwd[i] = Array<OneD, NekDouble>(nTracePts);
m_fields[i]->ExtractTracePhys(inarray[i], Fwd[i]);
}
std::string userDefStr;
int nreg = m_fields[0]->GetBndConditions().num_elements();
// Loop over Boundary Regions
......@@ -240,16 +249,16 @@ namespace Nektar
else if(boost::iequals(userDefStr, "IsentropicVortex"))
{
// Isentropic Vortex Boundary Condition
SetBoundaryIsentropicVortex(n, time, cnt, inarray);
SetBoundaryIsentropicVortex(n, time, cnt, Fwd, inarray);
}
else if (boost::iequals(userDefStr,"RinglebFlow"))
{
SetBoundaryRinglebFlow(n, time, cnt, inarray);
SetBoundaryRinglebFlow(n, time, cnt, Fwd, inarray);
}
else
{
// set up userdefined BC common to all solvers
SetCommonBC(userDefStr,n,time, cnt,inarray);
SetCommonBC(userDefStr, n, time, cnt, Fwd, inarray);
}
}
......@@ -390,21 +399,13 @@ namespace Nektar
void EulerCFE::SetBoundaryIsentropicVortex(
int bcRegion,
NekDouble time,
int cnt,
int cnt,
Array<OneD, Array<OneD, NekDouble> > &Fwd,
Array<OneD, Array<OneD, NekDouble> > &physarray)
{
int nvariables = physarray.num_elements();
int nTraceNumPoints = GetTraceTotPoints();
Array<OneD, Array<OneD, NekDouble> > Fwd(nvariables);
const Array<OneD, const int> &bndTraceMap = m_fields[0]->GetTraceBndMap();
// Get physical values of the forward trace (from exp to phys)
for (int i = 0; i < nvariables; ++i)
{
Fwd[i] = Array<OneD, NekDouble>(nTraceNumPoints);
m_fields[i]->ExtractTracePhys(physarray[i], Fwd[i]);
}
int id2, e_max;
e_max = m_fields[0]->GetBndCondExpansions()[bcRegion]->GetExpSize();
......@@ -836,11 +837,10 @@ namespace Nektar
int bcRegion,
NekDouble time,
int cnt,
Array<OneD, Array<OneD, NekDouble> > &Fwd,
Array<OneD, Array<OneD, NekDouble> > &physarray)
{
int nvariables = physarray.num_elements();
int nTraceNumPoints = GetTraceTotPoints();
Array<OneD, Array<OneD, NekDouble> > Fwd(nvariables);
// For 3DHomogenoeus1D
int n_planes = 1;
......@@ -849,14 +849,6 @@ namespace Nektar
int nPointsTot = m_fields[0]->GetTotPoints();
int nPointsTot_plane = m_fields[0]->GetPlane(0)->GetTotPoints();