35#include <boost/algorithm/string.hpp>
58 return int(x > 0.0 ? x + 0.5 : x - 0.5);
65 int size = inarray.size();
66 ASSERTL1(outarray.size() >= size,
"Array sizes not compatible");
69 for (
int i = 0; i < size; i++)
72 outarray[i] = int(x > 0.0 ? x + 0.5 : x - 0.5);
80 : m_session(), m_comm(), m_hash(0), m_numLocalBndCoeffs(0),
81 m_numGlobalBndCoeffs(0), m_numLocalDirBndCoeffs(0),
82 m_numGlobalDirBndCoeffs(0), m_solnType(
eNoSolnType),
83 m_bndSystemBandWidth(0), m_successiveRHS(0),
84 m_linSysIterSolver(
"ConjugateGradient"), m_gsh(nullptr), m_bndGsh(nullptr)
90 const std::string variable)
91 : m_session(pSession), m_comm(comm), m_variable(variable), m_hash(0),
92 m_numLocalBndCoeffs(0), m_numGlobalBndCoeffs(0),
93 m_numLocalDirBndCoeffs(0), m_numGlobalDirBndCoeffs(0),
94 m_bndSystemBandWidth(0), m_successiveRHS(0),
95 m_linSysIterSolver(
"ConjugateGradient"), m_gsh(nullptr), m_bndGsh(nullptr)
102 if (pSession->DefinesGlobalSysSolnInfo(variable,
"GlobalSysSoln"))
104 std::string sysSoln =
105 pSession->GetGlobalSysSolnInfo(variable,
"GlobalSysSoln");
107 "GlobalSysSoln", sysSoln);
112 if (pSession->DefinesGlobalSysSolnInfo(variable,
"Preconditioner"))
115 pSession->GetGlobalSysSolnInfo(variable,
"Preconditioner");
117 else if (pSession->DefinesSolverInfo(
"Preconditioner"))
119 m_preconType = pSession->GetSolverInfo(
"Preconditioner");
127 if (pSession->DefinesGlobalSysSolnInfo(variable,
128 "IterativeSolverTolerance"))
131 pSession->GetGlobalSysSolnInfo(variable,
"IterativeSolverTolerance")
136 pSession->LoadParameter(
"IterativeSolverTolerance",
141 if (pSession->DefinesGlobalSysSolnInfo(variable,
"AbsoluteTolerance"))
144 pSession->GetGlobalSysSolnInfo(variable,
"AbsoluteTolerance");
146 boost::iequals(boost::to_upper_copy(abstol),
"TRUE");
153 if (pSession->DefinesGlobalSysSolnInfo(variable,
"SuccessiveRHS"))
156 pSession->GetGlobalSysSolnInfo(variable,
"SuccessiveRHS").c_str());
163 if (pSession->DefinesGlobalSysSolnInfo(variable,
"LinSysIterSolver"))
166 pSession->GetGlobalSysSolnInfo(variable,
"LinSysIterSolver");
168 else if (pSession->DefinesSolverInfo(
"LinSysIterSolver"))
185 : m_session(oldLevelMap->m_session), m_comm(oldLevelMap->GetComm()),
186 m_hash(0), m_solnType(oldLevelMap->m_solnType),
187 m_preconType(oldLevelMap->m_preconType),
188 m_iterativeTolerance(oldLevelMap->m_iterativeTolerance),
189 m_successiveRHS(oldLevelMap->m_successiveRHS),
190 m_linSysIterSolver(oldLevelMap->m_linSysIterSolver),
191 m_gsh(oldLevelMap->m_gsh), m_bndGsh(oldLevelMap->m_bndGsh),
192 m_lowestStaticCondLevel(oldLevelMap->m_lowestStaticCondLevel)
214 int newLevel = staticCondLevelOld + 1;
226 numGlobalDirBndCoeffsOld);
236 multiLevelGraph->MaskPatches(newLevel, globHomPatchMask);
244 locPatchMask_NekDouble);
256 int numPatchesWithIntNew =
257 multiLevelGraph->GetNpatchesWithInterior(newLevel);
258 int numPatchesNew = numPatchesWithIntNew;
262 std::map<int, int> numLocalBndCoeffsPerPatchNew;
263 for (
int i = 0; i < numPatchesNew; i++)
265 numLocalBndCoeffsPerPatchNew[i] = 0;
271 for (i = cnt = 0; i < numPatchesOld; i++)
285 *min_element(&locPatchMask[cnt],
286 &locPatchMask[cnt] + numLocalBndCoeffsPerPatchOld[i]);
288 *max_element(&locPatchMask[cnt],
289 &locPatchMask[cnt] + numLocalBndCoeffsPerPatchOld[i]);
290 ASSERTL0((minval == maxval) || (minval == -1),
291 "These values should never be the same");
295 curPatch = numPatchesNew;
296 numLocalBndCoeffsPerPatchNew[curPatch] = 0;
304 for (j = 0; j < numLocalBndCoeffsPerPatchOld[i]; j++)
306 ASSERTL0((locPatchMask[cnt] == maxval) ||
307 (locPatchMask[cnt] == minval),
308 "These values should never be the same");
309 if (locPatchMask[cnt] == -1)
311 ++numLocalBndCoeffsPerPatchNew[curPatch];
329 (
unsigned int)numLocalBndCoeffsPerPatchNew[i];
341 "This method should only be called for in "
342 "case of multi-level static condensation.");
387 numLocalBndCoeffsPerPatchOffset[i] +=
388 numLocalBndCoeffsPerPatchOffset[i - 1] +
389 numLocalBndCoeffsPerPatchNew[i - 1];
392 int additionalPatchCnt = numPatchesWithIntNew;
399 for (i = cnt = 0; i < numPatchesOld; i++)
402 *min_element(&locPatchMask[cnt],
403 &locPatchMask[cnt] + numLocalBndCoeffsPerPatchOld[i]);
405 *max_element(&locPatchMask[cnt],
406 &locPatchMask[cnt] + numLocalBndCoeffsPerPatchOld[i]);
407 ASSERTL0((minval == maxval) || (minval == -1),
408 "These values should never be the same");
412 curPatch = additionalPatchCnt;
413 additionalPatchCnt++;
420 for (j = 0; j < numLocalBndCoeffsPerPatchOld[i]; j++)
422 ASSERTL0((locPatchMask[cnt] == maxval) ||
423 (locPatchMask[cnt] == minval),
424 "These values should never be the same");
428 if (locPatchMask[cnt] == -1)
430 newid = numLocalBndCoeffsPerPatchOffset[curPatch];
440 blockid = bndDofPerPatchCnt[curPatch];
443 numLocalBndCoeffsPerPatchOffset[curPatch]++;
444 bndDofPerPatchCnt[curPatch]++;
454 multiLevelGraph->GetInteriorOffset(newLevel, curPatch);
516 for (j = 0; j < locSize; j++)
531 bwidth = (bwidth > (maxId - minId)) ? bwidth : (maxId - minId);
552 [[maybe_unused]]
const int i)
const
580 [[maybe_unused]]
const int i)
const
596 [[maybe_unused]]
bool useComm)
const
604 [[maybe_unused]]
bool useComm)
const
653 [[maybe_unused]]
int offset)
const
715 [[maybe_unused]]
const ExpList &locexp,
719 static std::shared_ptr<AssemblyMap> result;
840 if (global.data() ==
loc.data())
860 if (global.data() ==
loc.data())
879 if (global.data() ==
loc.data())
1067 "Local vector is not of correct dimension");
1069 "Global vector is not of correct dimension");
1074 tmp.get() + offset, 1);
1092 "Local vector is not of correct dimension");
1094 "Global vector is not of correct dimension");
1097 if (global.data() ==
loc.data())
1123 "Local vector is not of correct dimension");
1125 "Global vector is not of correct dimension");
1155 "Local vector is not of correct dimension");
1157 "Global vector is not of correct dimension");
1179 "LocBnd vector is not of correct dimension");
1181 "Local vector is not of correct dimension");
1191 "Locint vector is not of correct dimension");
1193 "Local vector is not of correct dimension");
1203 "LocBnd vector is not of correct dimension");
1205 "Local vector is not of correct dimension");
1215 "LocBnd vector is not of correct dimension");
1217 "Local vector is not of correct dimension");
1239 "Local array is not of correct dimension");
1241 "Global array is not of correct dimension");
1263 "Local vector is not of correct dimension");
1265 "Global vector is not of correct dimension");
1389 "Local vector is not of correct dimension");
1391 "Global vector is not of correct dimension");
1398 bool printHeader)
const
1401 bool isRoot =
m_session->GetComm()->IsParallelInTime()
1403 : vRowComm->GetRank() == 0;
1404 int n = vRowComm->GetSize();
1408 int globBndCnt = 0, globDirCnt = 0;
1433 int totGlobDof = globCnt;
1434 int totGlobBndDof = globBndCnt;
1435 int totGlobDirDof = globDirCnt;
1440 int meanValence = 0;
1442 int minValence = 10000000;
1450 if (tmpGlob[i] > maxValence)
1452 maxValence = tmpGlob[i];
1454 if (tmpGlob[i] < minValence)
1456 minValence = tmpGlob[i];
1458 meanValence += tmpGlob[i];
1471 meanValence /= totGlobBndDof;
1477 out <<
"Assembly map statistics for field " << variable <<
":"
1481 out <<
" - Number of local/global dof : " << totLocalDof
1482 <<
" " << totGlobDof << endl;
1483 out <<
" - Number of local/global boundary dof : " << totLocalBndDof
1484 <<
" " << totGlobBndDof << endl;
1485 out <<
" - Number of local/global Dirichlet dof : " << totLocalDirDof
1486 <<
" " << totGlobDirDof << endl;
1487 out <<
" - dof valency (min/max/mean) : " << minValence
1488 <<
" " << maxValence <<
" " << meanValence << endl;
1496 for (i = 1; i < n; ++i)
1498 vRowComm->Recv(i, tmp);
1500 mean2 += tmp[0] * tmp[0];
1502 if (tmp[0] > maxval)
1506 if (tmp[0] < minval)
1514 out <<
" - Local dof dist. (min/max/mean/dev) : " << minval
1515 <<
" " << maxval <<
" " << (mean / n) <<
" "
1516 <<
sqrt(mean2 / n - mean * mean / n / n) << endl;
1522 mean2 = mean * mean;
1524 for (i = 1; i < n; ++i)
1526 vRowComm->Recv(i, tmp);
1528 mean2 += tmp[0] * tmp[0];
1530 if (tmp[0] > maxval)
1534 if (tmp[0] < minval)
1540 out <<
" - Local bnd dof dist. (min/max/mean/dev) : " << minval
1541 <<
" " << maxval <<
" " << (mean / n) <<
" "
1542 <<
sqrt(mean2 / n - mean * mean / n / n) << endl;
1549 vRowComm->Send(0, tmp);
1552 vRowComm->Send(0, tmp);
1564 while (tmp->m_nextLevelLocalToGlobalMap)
1566 tmp = tmp->m_nextLevelLocalToGlobalMap;
1579 for (i = 1; i < n; ++i)
1581 vRowComm->Recv(i, tmpRecv);
1583 mean2 += tmpRecv[0] * tmpRecv[0];
1585 if (tmpRecv[0] > maxval)
1587 maxval = (int)(tmpRecv[0] + 0.5);
1589 if (tmpRecv[0] < minval)
1591 minval = (int)(tmpRecv[0] + 0.5);
1595 out <<
" - M-level sc. dist. (min/max/mean/dev) : " << minval
1596 <<
" " << maxval <<
" " << (mean / n) <<
" "
1597 <<
sqrt(mean2 / n - mean * mean / n / n) << endl;
1603 vRowComm->Send(0, tmpSend);
1610 out <<
" - Number of static cond. levels : " <<
level
1617 out <<
"Stats at lowest static cond. level:" << endl;
1619 tmp->PrintStats(out, variable,
false);
#define ASSERTL0(condition, msg)
#define NEKERROR(type, msg)
Assert Level 0 – Fundamental assert which is used whether in FULLDEBUG, DEBUG or OPT compilation mode...
#define ASSERTL1(condition, msg)
Assert Level 1 – Debugging which is used whether in FULLDEBUG or DEBUG compilation mode....
#define sign(a, b)
return the sign(b)*a
static std::shared_ptr< DataType > AllocateSharedPtr(const Args &...args)
Allocate a shared pointer from the memory pool.
Base class for constructing local to global mapping of degrees of freedom.
void PatchLocalToGlobal(const Array< OneD, const NekDouble > &loc, Array< OneD, NekDouble > &global) const
const Array< OneD, const int > & GetExtraDirEdges()
GlobalSysSolnType m_solnType
The solution type of the global system.
int GetNumPatches() const
Returns the number of patches in this static condensation level.
const Array< OneD, const int > & GetGlobalToUniversalMapUnique()
int GetNumGlobalCoeffs() const
Returns the total number of global coefficients.
int m_numLocalCoeffs
Total number of local coefficients.
Array< OneD, int > m_bndCondCoeffsToLocalTraceMap
Integer map of bnd cond coeff to local trace coeff.
const Array< OneD, const unsigned int > & GetNumLocalBndCoeffsPerPatch()
Returns the number of local boundary coefficients in each patch.
Array< OneD, int > m_bndCondCoeffsToLocalCoeffsMap
Integer map of bnd cond coeffs to local coefficients.
virtual const Array< OneD, const int > & v_GetExtraDirEdges()
int GetNumDirFaces() const
int GetNumLocalCoeffs() const
Returns the total number of local coefficients.
std::string GetVariable()
Retrieves the variable string.
const Array< OneD, const unsigned int > & GetNumLocalIntCoeffsPerPatch()
Returns the number of local interior coefficients in each patch.
int GetNumNonDirFaceModes() const
bool IsAbsoluteTolerance() const
virtual void v_Assemble(const Array< OneD, const NekDouble > &loc, Array< OneD, NekDouble > &global) const
bool m_signChange
Flag indicating if modes require sign reversal.
void Assemble(const Array< OneD, const NekDouble > &loc, Array< OneD, NekDouble > &global) const
NekDouble GetLocalToGlobalBndSign(const int i) const
Retrieve the sign change of a given local boundary mode.
virtual std::shared_ptr< AssemblyMap > v_LinearSpaceMap(const ExpList &locexp, GlobalSysSolnType solnType)
Generate a linear space mapping from existing mapping.
void AssembleBnd(const NekVector< NekDouble > &loc, NekVector< NekDouble > &global, int offset) const
virtual void v_LocalToGlobal(const Array< OneD, const NekDouble > &loc, Array< OneD, NekDouble > &global, bool useComm) const
PatchMapSharedPtr m_patchMapFromPrevLevel
Mapping information for previous level in MultiLevel Solver.
Array< OneD, int > m_localToLocalIntMap
Integer map of local boundary coeffs to local interior system numbering.
virtual const Array< OneD, const int > & v_GetGlobalToUniversalMapUnique()
NekDouble m_iterativeTolerance
Tolerance for iterative solver.
NekDouble GetIterativeTolerance() const
int m_numGlobalCoeffs
Total number of global coefficients.
std::string m_linSysIterSolver
Iterative solver: Conjugate Gradient, GMRES.
bool AtLastLevel() const
Returns true if this is the last level in the multi-level static condensation.
void UniversalAssembleBnd(Array< OneD, NekDouble > &pGlobal) const
const Array< OneD, const int > & GetLocalToGlobalMap()
const Array< OneD, const int > & GetBndCondCoeffsToLocalTraceMap()
Retrieves the local indices corresponding to the boundary expansion modes to global trace.
int GetNumNonDirEdgeModes() const
virtual const Array< OneD, const int > & v_GetGlobalToUniversalMap()
Array< OneD, int > m_globalToUniversalBndMap
Integer map of process coeffs to universal space.
void PatchAssemble(const Array< OneD, const NekDouble > &loc, Array< OneD, NekDouble > &global) const
int m_successiveRHS
sucessive RHS for iterative solver
std::shared_ptr< AssemblyMap > LinearSpaceMap(const ExpList &locexp, GlobalSysSolnType solnType)
virtual int v_GetNumNonDirVertexModes() const
void CalculateBndSystemBandWidth()
Calculates the bandwidth of the boundary system.
void LocalToGlobal(const Array< OneD, const NekDouble > &loc, Array< OneD, NekDouble > &global, bool useComm=true) const
Array< OneD, NekDouble > m_localToGlobalBndSign
Integer sign of local boundary coeffs to global space.
int GetNumNonDirEdges() const
void UniversalAbsMaxBnd(Array< OneD, NekDouble > &bndvals)
Array< OneD, unsigned int > m_numLocalBndCoeffsPerPatch
The number of bnd dofs per patch.
virtual const Array< OneD, const int > & v_GetLocalToGlobalMap()
const Array< OneD, const int > & GetGlobalToUniversalMap()
int GetNumNonDirFaces() const
std::string GetLinSysIterSolver() const
virtual int v_GetNumNonDirEdgeModes() const
virtual int v_GetNumDirFaces() const
LibUtilities::SessionReaderSharedPtr m_session
Session object.
virtual const Array< OneD, NekDouble > & v_GetLocalToGlobalSign() const
int GetNumGlobalDirBndCoeffs() const
Returns the number of global Dirichlet boundary coefficients.
const Array< OneD, NekDouble > & GetBndCondCoeffsToLocalCoeffsSign()
Returns the modal sign associated with a given boundary expansion mode.
int m_numLocalBndCoeffs
Number of local boundary coefficients.
void LocalToLocalBnd(const Array< OneD, const NekDouble > &local, Array< OneD, NekDouble > &locbnd) const
const Array< OneD, const int > & GetBndCondCoeffsToLocalCoeffsMap()
Retrieves the local indices corresponding to the boundary expansion modes.
void PrintStats(std::ostream &out, std::string variable, bool printHeader=true) const
AssemblyMapSharedPtr m_nextLevelLocalToGlobalMap
Map from the patches of the previous level to the patches of the current level.
const Array< OneD, const int > & GetLocalToGlobalBndMap()
Retrieve the global indices of the local boundary modes.
int m_staticCondLevel
The level of recursion in the case of multi-level static condensation.
bool m_isAbsoluteTolerance
int GetSuccessiveRHS() const
int m_bndSystemBandWidth
The bandwith of the global bnd system.
const Array< OneD, const int > & GetGlobalToUniversalBndMapUnique()
void GlobalToLocal(const Array< OneD, const NekDouble > &global, Array< OneD, NekDouble > &loc) const
GlobalSysSolnType GetGlobalSysSolnType() const
Returns the method of solving global systems.
int m_numLocalDirBndCoeffs
Number of Local Dirichlet Boundary Coefficients.
bool GetSingularSystem() const
Retrieves if the system is singular (true) or not (false)
std::string m_variable
Variable string identifier.
virtual int v_GetNumNonDirFaceModes() const
void PatchGlobalToLocal(const Array< OneD, const NekDouble > &global, Array< OneD, NekDouble > &loc) const
virtual int v_GetFullSystemBandWidth() const
int m_numGlobalDirBndCoeffs
Number of Global Dirichlet Boundary Coefficients.
LibUtilities::CommSharedPtr GetComm()
Retrieves the communicator.
Array< OneD, unsigned int > m_numLocalIntCoeffsPerPatch
The number of int dofs per patch.
const AssemblyMapSharedPtr GetNextLevelLocalToGlobalMap() const
Returns the local to global mapping for the next level in the multi-level static condensation.
int GetBndSystemBandWidth() const
Returns the bandwidth of the boundary system.
virtual void v_UniversalAssemble(Array< OneD, NekDouble > &pGlobal) const
int GetStaticCondLevel() const
Returns the level of static condensation for this map.
void GlobalToLocalBndWithoutSign(const Array< OneD, const NekDouble > &global, Array< OneD, NekDouble > &loc)
int GetNumNonDirVertexModes() const
bool m_systemSingular
Flag indicating if the system is singular or not.
Array< OneD, int > m_localToGlobalBndMap
Integer map of local coeffs to global Boundary Dofs.
size_t GetHash() const
Retrieves the hash of this map.
int GetLocalToGlobalBndMap(const int i) const
Retrieve the global index of a given local boundary mode.
virtual ~AssemblyMap()
Destructor.
Gs::gs_data * m_dirBndGsh
gs gather communication to impose Dirhichlet BCs.
virtual void v_GlobalToLocal(const Array< OneD, const NekDouble > &global, Array< OneD, NekDouble > &loc) const
const Array< OneD, NekDouble > & GetLocalToGlobalSign() const
int GetNumLocalDirBndCoeffs() const
Returns the number of local Dirichlet boundary coefficients.
int GetNumDirEdges() const
size_t m_hash
Hash for map.
Array< OneD, const NekDouble > GetLocalToGlobalBndSign() const
Retrieve the sign change for all local boundary modes.
void LocalBndToLocal(const Array< OneD, const NekDouble > &locbnd, Array< OneD, NekDouble > &local) const
AssemblyMap()
Default constructor.
Array< OneD, int > m_globalToUniversalBndMapUnique
Integer map of unique process coeffs to universal space (signed)
virtual int v_GetNumNonDirEdges() const
virtual int v_GetNumDirEdges() const
const PatchMapSharedPtr & GetPatchMapFromPrevLevel(void) const
Returns the patch map from the previous level of the multi-level static condensation.
LibUtilities::CommSharedPtr m_comm
Communicator.
int GetNumGlobalBndCoeffs() const
Returns the total number of global boundary coefficients.
int GetFullSystemBandWidth() const
const Array< OneD, const int > & GetGlobalToUniversalBndMap()
Array< OneD, int > m_localToLocalBndMap
Integer map of local boundary coeffs to local boundary system numbering.
std::string GetPreconType() const
Array< OneD, NekDouble > m_bndCondCoeffsToLocalCoeffsSign
Integer map of sign of bnd cond coeffs to local coefficients.
const Array< OneD, const int > & GetBndCondIDToGlobalTraceID()
void LocalBndToGlobal(const Array< OneD, const NekDouble > &loc, Array< OneD, NekDouble > &global, int offset, bool UseComm=true) const
void UniversalAssemble(Array< OneD, NekDouble > &pGlobal) const
int m_numPatches
The number of patches (~elements) in the current level.
void LocalToLocalInt(const Array< OneD, const NekDouble > &local, Array< OneD, NekDouble > &locint) const
virtual int v_GetNumNonDirFaces() const
bool GetSignChange()
Returns true if using a modal expansion requiring a change of sign of some modes.
void LocalIntToLocal(const Array< OneD, const NekDouble > &locbnd, Array< OneD, NekDouble > &local) const
void GlobalToLocalBnd(const NekVector< NekDouble > &global, NekVector< NekDouble > &loc, int offset) const
std::string m_preconType
Type type of preconditioner to use in iterative solver.
int GetNumLocalBndCoeffs() const
Returns the total number of local boundary coefficients.
int m_numGlobalBndCoeffs
Total number of global boundary coefficients.
Array< OneD, int > m_bndCondIDToGlobalTraceID
Integer map of bnd cond trace number to global trace number.
Base class for all multi-elemental spectral/hp expansions.
Array< OneD, DataType > & GetPtr()
static void Gather(Nektar::Array< OneD, NekDouble > pU, gs_op pOp, gs_data *pGsh, Nektar::Array< OneD, NekDouble > pBuffer=NullNekDouble1DArray)
Performs a gather-scatter operation of the provided values.
std::shared_ptr< SessionReader > SessionReaderSharedPtr
std::shared_ptr< Comm > CommSharedPtr
Pointer to a Communicator object.
@ eIterativeMultiLevelStaticCond
@ eNoSolnType
No Solution type specified.
@ eDirectMultiLevelStaticCond
@ eXxtMultiLevelStaticCond
@ ePETScMultiLevelStaticCond
std::shared_ptr< BottomUpSubStructuredGraph > BottomUpSubStructuredGraphSharedPtr
std::shared_ptr< AssemblyMap > AssemblyMapSharedPtr
std::shared_ptr< PatchMap > PatchMapSharedPtr
int RoundNekDoubleToInt(NekDouble x)
Rounds a double precision number to an integer.
static const NekDouble kNekIterativeTol
void Gathr(I n, const T *x, const I *y, T *z)
Gather vector z[i] = x[y[i]].
void Scatr(int n, const T *x, const int *y, T *z)
Scatter vector z[y[i]] = x[i].
void Assmb(int n, const T *x, const int *y, T *z)
Assemble z[y[i]] += x[i]; z should be zero'd first.
void Zero(int n, T *x, const int incx)
Zero vector.
void Vcopy(int n, const T *x, const int incx, T *y, const int incy)
scalarT< T > sqrt(scalarT< T > in)