I encountered a problem when running a 3DH2D incompressible solver. Since both boundaries are Neumann type for the pressure field, the solver fails with:
Fatal: Level 0 assertion violation
Where: LibUtilities/LinearAlgebra/NekLinSys.hpp[661]
Message: ERROR: The leading minor of order 5 is not positive definite from dpbtrf
This appears to be caused by the singular ((k_y,k_z)=(0,0)) Fourier mode of the
Poisson solver.
Maybe try to default to a generic banded case for that specific check to see if it solve the problem?
The singularity happens in the (ky=0, kz=0) mode. I asked codex to remove the singularity in this mode, and there is a solution. We can test this after the server recovers.
diff --git a/library/MultiRegions/ContField3DHomogeneous2D.cpp b/library/MultiRegions/ContField3DHomogeneous2D.cpp
index 3374f8dc69..a7f69deab2 100644
--- a/library/MultiRegions/ContField3DHomogeneous2D.cpp
+++ b/library/MultiRegions/ContField3DHomogeneous2D.cpp
@@ -49,12 +49,11 @@ ContField3DHomogeneous2D::ContField3DHomogeneous2D(
: DisContField3DHomogeneous2D(In, false)
{
- ContFieldSharedPtr zero_line =
- std::dynamic_pointer_cast<ContField>(In.m_lines[0]);
-
for (int n = 0; n < m_lines.size(); ++n)
{
- m_lines[n] = MemoryManager<ContField>::AllocateSharedPtr(*zero_line);
+ ContFieldSharedPtr line =
+ std::dynamic_pointer_cast<ContField>(In.m_lines[n]);
+ m_lines[n] = MemoryManager<ContField>::AllocateSharedPtr(*line);
}
SetCoeffPhys();
@@ -70,7 +69,8 @@ ContField3DHomogeneous2D::ContField3DHomogeneous2D(
const LibUtilities::BasisKey &HomoBasis_z, const NekDouble lhom_y,
const NekDouble lhom_z, const bool useFFT, const bool dealiasing,
const SpatialDomains::MeshGraphSharedPtr &graph1D,
- const std::string &variable, const Collections::ImplementationType ImpType)
+ const std::string &variable, const bool CheckIfSingularSystem,
+ const Collections::ImplementationType ImpType)
: DisContField3DHomogeneous2D(pSession, HomoBasis_y, HomoBasis_z, lhom_y,
lhom_z, useFFT, dealiasing, ImpType)
{
@@ -79,7 +79,7 @@ ContField3DHomogeneous2D::ContField3DHomogeneous2D(
SpatialDomains::BoundaryConditions bcs(pSession, graph1D);
m_lines[0] = line_zero = MemoryManager<ContField>::AllocateSharedPtr(
- pSession, graph1D, variable, false, false, ImpType);
+ pSession, graph1D, variable, false, CheckIfSingularSystem, ImpType);
m_exp = MemoryManager<LocalRegions::ExpansionVector>::AllocateSharedPtr();
nel = m_lines[0]->GetExpSize();
@@ -196,6 +196,15 @@ GlobalLinSysKey ContField3DHomogeneous2D::v_HelmSolve(
{
for (m = 0; m < nhom_modes_y; ++m, l++)
{
+ // Index 1 is the null Fourier mode and is set to zero by the
+ // forward and backward transforms.
+ if (m == 1 || n == 1)
+ {
+ cnt += m_lines[l]->GetTotPoints();
+ cnt1 += m_lines[l]->GetNcoeffs();
+ continue;
+ }
+
beta_z = 2 * M_PI * (n / 2) / m_lhom_z;
beta_y = 2 * M_PI * (m / 2) / m_lhom_y;
beta = beta_y * beta_y + beta_z * beta_z;
diff --git a/library/MultiRegions/ContField3DHomogeneous2D.h b/library/MultiRegions/ContField3DHomogeneous2D.h
index 304dd14d6a..b23465f6a3 100644
--- a/library/MultiRegions/ContField3DHomogeneous2D.h
+++ b/library/MultiRegions/ContField3DHomogeneous2D.h
@@ -52,7 +52,7 @@ public:
const LibUtilities::BasisKey &HomoBasis_z, const NekDouble lhom_y,
const NekDouble lhom_z, const bool useFFT, const bool dealiasing,
const SpatialDomains::MeshGraphSharedPtr &graph1D,
- const std::string &variable,
+ const std::string &variable, const bool CheckIfSingularSystem = false,
const Collections::ImplementationType ImpType =
Collections::eNoImpType);
diff --git a/library/SolverUtils/EquationSystem.cpp b/library/SolverUtils/EquationSystem.cpp
index de6bfea655..aef0f15e42 100644
--- a/library/SolverUtils/EquationSystem.cpp
+++ b/library/SolverUtils/EquationSystem.cpp
@@ -330,7 +330,8 @@ void EquationSystem::v_InitObject(bool DeclareFields)
AllocateSharedPtr(m_session, BkeyY, BkeyZ,
m_LhomY, m_LhomZ, m_useFFT,
m_homogen_dealiasing, m_graph,
- m_session->GetVariable(i));
+ m_session->GetVariable(i),
+ m_checkIfSystemSingular[i]);
}
}
else
```