Commit 149fbcb6 authored by Chris Cantwell's avatar Chris Cantwell
Browse files

Changed singular mesh vertex to be that of minimum ID for parallel consistency.

parent 3ce6a86d
...@@ -838,6 +838,7 @@ namespace Vmath ...@@ -838,6 +838,7 @@ namespace Vmath
} }
template LIB_UTILITIES_EXPORT int Imin( int n, const Nektar::NekDouble *x, const int incx); template LIB_UTILITIES_EXPORT int Imin( int n, const Nektar::NekDouble *x, const int incx);
template LIB_UTILITIES_EXPORT int Imin( int n, const int *x, const int incx);
/// \brief Return the minimum element in x - called vmin to avoid /// \brief Return the minimum element in x - called vmin to avoid
/// conflict with min /// conflict with min
......
...@@ -658,13 +658,13 @@ namespace Nektar ...@@ -658,13 +658,13 @@ namespace Nektar
vCommRow->AllReduce(s, LibUtilities::ReduceMin); vCommRow->AllReduce(s, LibUtilities::ReduceMin);
systemSingular = (s == 1 ? true : false); systemSingular = (s == 1 ? true : false);
// Count the number of boundary regions on each process // Find the minimum boundary vertex ID on each process
Array<OneD, int> bccounts(n, 0); Array<OneD, int> bcminvertid(n, 0);
bccounts[p] = bndCondExp.num_elements(); bcminvertid[p] = Dofs[0].begin()->first;
vCommRow->AllReduce(bccounts, LibUtilities::ReduceSum); vCommRow->AllReduce(bcminvertid, LibUtilities::ReduceSum);
// Find the process rank with the maximum number of boundary regions // Find the process rank with the minimum boundary vertex ID
int maxBCIdx = Vmath::Imax(n, bccounts, 1); int minIdx = Vmath::Imin(n, bcminvertid, 1);
// If the system is singular, the process with the maximum number of // If the system is singular, the process with the maximum number of
// BCs will set a Dirichlet vertex to make system non-singular. // BCs will set a Dirichlet vertex to make system non-singular.
...@@ -672,7 +672,7 @@ namespace Nektar ...@@ -672,7 +672,7 @@ namespace Nektar
// we do not try to set a Dirichlet vertex on a partition with no // we do not try to set a Dirichlet vertex on a partition with no
// intersection with the boundary. // intersection with the boundary.
meshVertId = 0; meshVertId = 0;
if(systemSingular == true && checkIfSystemSingular && maxBCIdx == p) if(systemSingular == true && checkIfSystemSingular && minIdx == p)
{ {
if (m_session->DefinesParameter("SingularElement")) if (m_session->DefinesParameter("SingularElement"))
{ {
...@@ -699,9 +699,9 @@ namespace Nektar ...@@ -699,9 +699,9 @@ namespace Nektar
} }
else else
{ {
bndSegExp = bndCondExp[bndCondExp.num_elements()-1] // Set pinned vertex to that with minimum vertex ID to
->GetExp(0)->as<LocalRegions::SegExp>(); // ensure consistency in parallel.
meshVertId = bndSegExp->GetGeom1D()->GetVid(0); meshVertId = bcminvertid[p];
} }
if (ReorderedGraphVertId[0].count(meshVertId) == 0) if (ReorderedGraphVertId[0].count(meshVertId) == 0)
...@@ -720,7 +720,7 @@ namespace Nektar ...@@ -720,7 +720,7 @@ namespace Nektar
// Firstly, we check that no other processors have this // Firstly, we check that no other processors have this
// vertex. If they do, then we mark the vertex as also being // vertex. If they do, then we mark the vertex as also being
// Dirichlet. // Dirichlet.
if (maxBCIdx != p) if (minIdx != p)
{ {
if (Dofs[0].count(meshVertId) > 0) if (Dofs[0].count(meshVertId) > 0)
{ {
...@@ -810,7 +810,7 @@ namespace Nektar ...@@ -810,7 +810,7 @@ namespace Nektar
Array<OneD, int> offsets(n, 0); Array<OneD, int> offsets(n, 0);
counts[p] = ReorderedGraphVertId[0].size(); counts[p] = ReorderedGraphVertId[0].size();
vCommRow->AllReduce(counts, LibUtilities::ReduceSum); vCommRow->AllReduce(counts, LibUtilities::ReduceSum);
for (i = 1; i < n; ++i) for (i = 1; i < n; ++i)
{ {
offsets[i] = offsets[i-1] + counts[i-1]; offsets[i] = offsets[i-1] + counts[i-1];
......
...@@ -662,14 +662,14 @@ namespace Nektar ...@@ -662,14 +662,14 @@ namespace Nektar
} }
} }
m_locToGloMap->UniversalAssembleBnd(tmp); m_locToGloMap->UniversalAssembleBnd(tmp);
// Now fill in all other Dirichlet coefficients. // Now fill in all other Dirichlet coefficients.
for(i = 0; i < m_bndCondExpansions.num_elements(); ++i) for(i = 0; i < m_bndCondExpansions.num_elements(); ++i)
{ {
if(m_bndConditions[i]->GetBoundaryConditionType() == if(m_bndConditions[i]->GetBoundaryConditionType() ==
SpatialDomains::eDirichlet) SpatialDomains::eDirichlet)
{ {
const Array<OneD,const NekDouble>& coeffs = const Array<OneD,const NekDouble>& coeffs =
m_bndCondExpansions[i]->GetCoeffs(); m_bndCondExpansions[i]->GetCoeffs();
for(j = 0; j < (m_bndCondExpansions[i])->GetNcoeffs(); ++j) for(j = 0; j < (m_bndCondExpansions[i])->GetNcoeffs(); ++j)
{ {
...@@ -683,7 +683,7 @@ namespace Nektar ...@@ -683,7 +683,7 @@ namespace Nektar
bndcnt += m_bndCondExpansions[i]->GetNcoeffs(); bndcnt += m_bndCondExpansions[i]->GetNcoeffs();
} }
} }
Vmath::Vcopy(nDir, tmp, 1, outarray, 1); Vmath::Vcopy(nDir, tmp, 1, outarray, 1);
} }
......
...@@ -279,7 +279,7 @@ namespace Nektar ...@@ -279,7 +279,7 @@ namespace Nektar
// Evaluate convection terms // Evaluate convection terms
m_advObject->DoAdvection(m_fields, m_nConvectiveFields, m_velocity, m_advObject->DoAdvection(m_fields, m_nConvectiveFields, m_velocity,
inarray, outarray, m_time); inarray, outarray, m_time);
// Smooth advection // Smooth advection
if(m_SmoothAdvection) if(m_SmoothAdvection)
{ {
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment