Commit 2657c623 authored by Chris Cantwell's avatar Chris Cantwell
Browse files

Reimplementation of CG solver to resolve numerical instability issues.

Reverted regression test updates.
Reverted default solver tolerance back to 1e-09.


git-svn-id: https://gforge.sci.utah.edu/svn/nektar/trunk@4004 305cdda6-5ce1-45b3-a98d-dfc68c8b3305
parent acc00705
......@@ -48,7 +48,7 @@ namespace Nektar
static const NekDouble kNekZeroTol = 1.0e-12;
static const NekDouble kGeomRightAngleTol = 1e-14;
static const NekDouble kNekSqrtTol = 1.0e-16;
static const NekDouble kNekIterativeTol = 1e-06;
static const NekDouble kNekIterativeTol = 1e-09;
}
} //end of namespace
......
......@@ -74,7 +74,7 @@ namespace Nektar
* parallel execution of the solver.
*
* The implemented algorithm uses a reduced-communication reordering of
* the standard PCG method (Amin, Sadayappan, Gudavalli, IEEE 1994)
* the standard PCG method (Demmel, Heath and Vorst, 1993)
*
* @param pInput Input vector of all DOFs.
* @param pOutput Solution vector of all DOFs.
......@@ -102,33 +102,47 @@ namespace Nektar
int nNonDir = nGlobal - nDir;
// Allocate array storage
Array<OneD, NekDouble> p_A (nGlobal, 0.0);
Array<OneD, NekDouble> q_A (nGlobal, 0.0);
Array<OneD, NekDouble> y_A (nNonDir, 0.0);
Array<OneD, NekDouble> w_A (nGlobal, 0.0);
Array<OneD, NekDouble> s_A (nGlobal, 0.0);
Array<OneD, NekDouble> p_A (nNonDir, 0.0);
Array<OneD, NekDouble> r_A (nNonDir, 0.0);
Array<OneD, NekDouble> z_A (nNonDir, 0.0);
Array<OneD, NekDouble> q_A (nNonDir, 0.0);
Array<OneD, NekDouble> tmp;
// Create NekVector wrappers for linear algebra operations
NekVector<NekDouble> in (nNonDir,pInput + nDir, eWrapper);
NekVector<NekDouble> in (nNonDir,pInput + nDir,eWrapper);
NekVector<NekDouble> out(nNonDir,pOutput + nDir,eWrapper);
NekVector<NekDouble> p (nNonDir,p_A + nDir, eWrapper);
NekVector<NekDouble> q (nNonDir,q_A + nDir, eWrapper);
NekVector<NekDouble> y (nNonDir,y_A, eWrapper);
NekVector<NekDouble> w (nNonDir,w_A + nDir, eWrapper);
NekVector<NekDouble> s (nNonDir,s_A + nDir, eWrapper);
NekVector<NekDouble> p (nNonDir,p_A, eWrapper);
NekVector<NekDouble> r (nNonDir,r_A, eWrapper);
NekVector<NekDouble> z (nNonDir,z_A, eWrapper);
NekVector<NekDouble> q (nNonDir,q_A, eWrapper);
int k;
NekDouble alpha, beta, r_dot_d, b_dot_b, min_resid;
NekDouble alpha, beta, rho, rho_new, mu, eps, bb_inv, min_resid;
Array<OneD, NekDouble> vExchange(3);
// Initialise with zero as the initial guess.
r = in;
m_precon->DoPreconditioner(r_A,z_A);
p = z;
m_precon->DoPreconditioner(r_A, tmp = w_A + nDir);
v_DoMatrixMultiply(w_A, s_A);
k = 0;
vExchange[0] = Vmath::Dot2(nNonDir, r_A, r_A, m_map + nDir);
vExchange[1] = Vmath::Dot2(nNonDir, r_A, z_A, m_map + nDir);
vExchange[0] = Vmath::Dot2(nNonDir,
r_A,
w_A + nDir,
m_map + nDir);
vExchange[1] = Vmath::Dot2(nNonDir,
s_A + nDir,
w_A + nDir,
m_map + nDir);
vExchange[2] = Vmath::Dot2(nNonDir,
r_A,
r_A,
m_map + nDir);
vComm->AllReduce(vExchange, Nektar::LibUtilities::ReduceSum);
// If input vector is zero, set zero output and skip solve.
......@@ -137,9 +151,14 @@ namespace Nektar
Vmath::Zero(nGlobal, pOutput, 1);
return;
}
b_dot_b = vExchange[0];
r_dot_d = vExchange[1];
min_resid = b_dot_b;
rho = vExchange[0];
mu = vExchange[1];
beta = 0.0;
alpha = rho/mu;
eps = 0.0;
bb_inv = 1.0/vExchange[2];
min_resid = bb_inv;
// Continue until convergence
while (true)
......@@ -147,63 +166,207 @@ namespace Nektar
ASSERTL0(k < 20000,
"Exceeded maximum number of iterations (20000)");
ASSERTL0(vExchange[0] <= b_dot_b || k < 10,
ASSERTL0(eps*bb_inv <= 1.0 || k < 10,
"Conjugate gradient diverged. Tolerance too small?"
"Minimum residual achieved: "
+ boost::lexical_cast<std::string>(sqrt(min_resid)));
// Perform the method-specific matrix-vector multiply operation.
v_DoMatrixMultiply(p_A, q_A);
// Compute new search direction p_k, q_k
p = w + beta * p;
q = s + beta * q;
// Update solution x_{k+1}
out = out + alpha * p;
// Update residual vector r_{k+1}
r = r - alpha * q;
// Apply preconditioner
m_precon->DoPreconditioner(q_A + nDir, y_A);
m_precon->DoPreconditioner(r_A, tmp = w_A + nDir);
// Perform the method-specific matrix-vector multiply operation.
v_DoMatrixMultiply(w_A, s_A);
// <r_k, r_k>
// <r_{k+1}, w_{k+1}>
vExchange[0] = Vmath::Dot2(nNonDir,
r_A,
r_A,
w_A + nDir,
m_map + nDir);
// <p_k, q_k>
// <s_{k+1}, w_{k+1}>
vExchange[1] = Vmath::Dot2(nNonDir,
p_A + nDir,
q_A + nDir,
s_A + nDir,
w_A + nDir,
m_map + nDir);
// <q_k, y_k>
// <r_{k+1}, r_{k+1}>
vExchange[2] = Vmath::Dot2(nNonDir,
q_A + nDir,
y_A,
r_A,
r_A,
m_map + nDir);
// Perform exchange
// Perform inner-product exchanges
vComm->AllReduce(vExchange, Nektar::LibUtilities::ReduceSum);
min_resid = min(min_resid, vExchange[0]);
rho_new = vExchange[0];
mu = vExchange[1];
eps = vExchange[2];
// test if norm is within tolerance
if (vExchange[0] < m_tolerance * m_tolerance)
if (eps*bb_inv < m_tolerance * m_tolerance)
{
break;
}
min_resid = min(min_resid, eps);
// Compute search direction and solution coefficients
beta = rho_new/rho;
alpha = rho_new/(mu - rho_new*beta/alpha);
rho = rho_new;
k++;
}
}
#if 0
/**
* Solve a global linear system using the conjugate gradient method.
* We solve only for the non-Dirichlet modes. The operator is evaluated
* using the local-matrix representation. Distributed math routines are
* used to support parallel execution of the solver.
*
* Standard CG algorithm.
*
* @param pInput Input vector of non-Dirichlet DOFs.
* @param pOutput Solution vector of non-Dirichlet DOFs.
*/
void GlobalLinSysIterative::v_SolveLinearSystem(
const int nGlobal,
const Array<OneD,const NekDouble> &pInput,
Array<OneD, NekDouble> &pOutput,
const AssemblyMapSharedPtr &plocToGloMap,
const int nDir)
{
// Check if preconditioner has been computed and compute if needed.
if (!m_precon)
{
v_UniqueMap();
m_precon = MemoryManager<Preconditioner>::AllocateSharedPtr(GetSharedThisPtr(),plocToGloMap);
}
// Get the communicator for performing data exchanges
LibUtilities::CommSharedPtr vComm = m_expList.lock()->GetComm()->GetRowComm();
// Get vector sizes
int nNonDir = nGlobal - nDir;
// Allocate array storage
Array<OneD, NekDouble> d_A (nGlobal, 0.0);
Array<OneD, NekDouble> p_A (nGlobal, 0.0);
Array<OneD, NekDouble> z_A (nNonDir, 0.0);
Array<OneD, NekDouble> z_new_A(nNonDir, 0.0);
Array<OneD, NekDouble> r_A (nNonDir, 0.0);
Array<OneD, NekDouble> r_new_A(nNonDir, 0.0);
// Create NekVector wrappers for linear algebra operations
NekVector<NekDouble> in(nNonDir,pInput + nDir,eWrapper);
NekVector<NekDouble> out(nNonDir,pOutput + nDir,eWrapper);
NekVector<NekDouble> r(nNonDir,r_A,eWrapper);
NekVector<NekDouble> r_new(nNonDir,r_new_A,eWrapper);
NekVector<NekDouble> z(nNonDir,z_A,eWrapper);
NekVector<NekDouble> z_new(nNonDir,z_new_A,eWrapper);
NekVector<NekDouble> d(nNonDir,d_A + nDir, eWrapper);
NekVector<NekDouble> p(nNonDir,p_A + nDir,eWrapper);
int k;
NekDouble alpha, beta, normsq, r_dot_z_old;
Array<OneD, NekDouble> vExchange(2);
// INVERSE of preconditioner matrix.
//const DNekMat &M = (*m_preconditioner);
// Initialise with zero as the initial guess.
r = in;
m_precon->DoPreconditioner(r_A,z_A);
d = z;
k = 0;
vExchange[0] = Vmath::Dot2(nNonDir, r_A, r_A, m_map + nDir);
vExchange[1] = Vmath::Dot2(nNonDir, r_A, z_A, m_map + nDir);
vComm->AllReduce(vExchange, Nektar::LibUtilities::ReduceSum);
// If input vector is zero, set zero output and skip solve.
if (vExchange[0] < NekConstants::kNekZeroTol)
{
Vmath::Zero(nGlobal, pOutput, 1);
return;
}
r_dot_z_old = vExchange[1];
// Continue until convergence
while (true)
{
// Perform the method-specific matrix-vector multiply operation.
v_DoMatrixMultiply(d_A, p_A);
// compute step length alpha
// alpha denominator
vExchange[0] = Vmath::Dot2(nNonDir,
d_A + nDir,
p_A + nDir,
m_map + nDir);
// perform exchange
vComm->AllReduce(vExchange, Nektar::LibUtilities::ReduceSum);
// compute alpha
alpha = r_dot_d / vExchange[1];
beta = alpha * vExchange[2] / vExchange[1] - 1.0;
r_dot_d *= beta;
alpha = r_dot_z_old/vExchange[0];
// approximate solution
out = out + alpha*d;
// Update residual vector
r = r - alpha * q;
z = z - alpha * y;
// compute residual
r_new = r - alpha*p;
// Update solution
out = out + alpha * p;
// Apply preconditioner to new residual
m_precon->DoPreconditioner(r_new_A,z_new_A);
// beta
vExchange[0] = Vmath::Dot2(nNonDir,
r_new_A,
z_new_A,
m_map + nDir);
// residual
vExchange[1] = Vmath::Dot2(nNonDir,
r_new_A,
r_new_A,
m_map + nDir);
// perform exchange
vComm->AllReduce(vExchange, Nektar::LibUtilities::ReduceSum);
// extract values for beta and norm
beta = vExchange[0]/r_dot_z_old;
r_dot_z_old = vExchange[0];
normsq = vExchange[1];
// test if norm is within tolerance
if (normsq < m_tolerance * m_tolerance)
{
break;
}
// Compute new search direction
p = z + beta * p;
d = z_new + beta*d;
// Next step
r = r_new;
z = z_new;
k++;
ASSERTL0(k < 20000,
"Exceeded maximum number of iterations (20000)");
}
}
#endif
}
}
......@@ -160,7 +160,6 @@
<CONDITIONS>
<PARAMETERS>
<P> Lambda = 1 </P>
<P> IterativeSolverTolerance = 5e-08 </P>
</PARAMETERS>
<VARIABLES>
......
......@@ -24,7 +24,6 @@
<P> Kinvis = 1 </P>
<P> HomModesZ = 4 </P>
<P> LZ = 1.0 </P>
<P> IterativeSolverTolerance = 5e-08 </P>
</PARAMETERS>
<VARIABLES>
......
......@@ -25,7 +25,6 @@
<P> HomModesZ = 4 </P>
<P> LZ = 1.0 </P>
<P> PROC_Z = 2 </P>
<P> IterativeSolverTolerance = 5e-08 </P>
</PARAMETERS>
<VARIABLES>
......
......@@ -22,7 +22,6 @@
<P> IO_CheckSteps = 1000 </P>
<P> IO_InfoSteps = 1000 </P>
<P> Kinvis = 1 </P>
<P> IterativeSolverTolerance = 5e-08 </P>
</PARAMETERS>
<VARIABLES>
......
......@@ -163,7 +163,7 @@
<P> IO_CheckSteps = 10 </P>
<P> IO_InfoSteps = 10 </P>
<P> Kinvis = 1 </P>
<P> IterativeSolverTolerance = 2e-07 </P>
<P> IterativeSolverTolerance = 1e-10 </P>
</PARAMETERS>
<VARIABLES>
......
......@@ -331,7 +331,7 @@
<P> IO_CheckSteps = 10 </P>
<P> IO_InfoSteps = 10 </P>
<P> Kinvis = 1 </P>
<P> IterativeSolverTolerance = 5e-07 </P>
<P> IterativeSolverTolerance = 1e-10 </P>
</PARAMETERS>
<VARIABLES>
......
......@@ -655,9 +655,8 @@
<P> TimeStep = 0.001 </P>
<P> NumSteps = 10 </P>
<P> IO_CheckSteps = 100 </P>
<P> IO_InfoSteps = 1 </P>
<P> IO_InfoSteps = 10 </P>
<P> Kinvis = 1 </P>
<P> IterativeSolverTolerance = 1e-07 </P>
</PARAMETERS>
<VARIABLES>
......
IncNavierStokesSolver
==========================================
Test_ChanFlow_m3.xml
L 2 error (variable u) : 1.06003e-16
L inf error (variable u) : 8.32667e-16
L 2 error (variable u) : 4.80863e-16
L inf error (variable u) : 4.32987e-15
L 2 error (variable v) : 0
L inf error (variable v) : 1.69891e-16
L 2 error (variable p) : 1.3259e-14
L inf error (variable p) : 2.73115e-14
L inf error (variable v) : 2.54809e-16
L 2 error (variable p) : 8.40089e-15
L inf error (variable p) : 4.9738e-14
----------------------------------------
Test_ChanFlow_m8_BodyForce.xml
L 2 error (variable u) : 9.90432e-06
L inf error (variable u) : 1.40068e-05
L 2 error (variable v) : 0
L inf error (variable v) : 1.97885e-16
L 2 error (variable p) : 2.98061e-14
L inf error (variable p) : 5.26916e-14
L inf error (variable v) : 5.41638e-16
L 2 error (variable p) : 7.49767e-15
L inf error (variable p) : 4.90237e-14
----------------------------------------
Test_ChanFlow_m8.xml
L 2 error (variable u) : 2.65566e-16
L inf error (variable u) : 8.07687e-15
L 2 error (variable v) : 1.62009e-16
L inf error (variable v) : 6.28378e-16
L 2 error (variable p) : 1.11902e-14
L inf error (variable p) : 3.30846e-14
L 2 error (variable u) : 5.32742e-16
L inf error (variable u) : 1.21847e-14
L 2 error (variable v) : 1.411e-16
L inf error (variable v) : 7.88839e-16
L 2 error (variable p) : 2.10864e-14
L inf error (variable p) : 5.55112e-14
----------------------------------------
Test_ChanFlow2D_bcsfromfiles.xml
L 2 error (variable u) : 1.57986e-13
L inf error (variable u) : 3.30402e-13
L 2 error (variable v) : 1.10243e-13
L inf error (variable v) : 1.05568e-13
L 2 error (variable p) : 1.32919e-11
L inf error (variable p) : 1.87586e-11
L 2 error (variable u) : 1.77702e-13
L inf error (variable u) : 7.49623e-13
L 2 error (variable v) : 2.65129e-13
L inf error (variable v) : 2.49849e-13
L 2 error (variable p) : 4.09367e-11
L inf error (variable p) : 1.01518e-10
----------------------------------------
Test_KovaFlow_m3.xml
L 2 error (variable u) : 0.343606
......@@ -71,18 +71,18 @@ L 2 error (variable v) : 0.000172531
L inf error (variable v) : 0.000353995
----------------------------------------
Test_ChanFlow_m8_singular.xml
L 2 error (variable u) : 2.7648e-15
L inf error (variable u) : 1.43358e-14
L 2 error (variable v) : 3.16374e-15
L inf error (variable v) : 1.24069e-14
L 2 error (variable p) : 1.48815e-12
L inf error (variable p) : 1.75965e-12
L 2 error (variable u) : 3.08173e-15
L inf error (variable u) : 1.60011e-14
L 2 error (variable v) : 3.79662e-15
L inf error (variable v) : 1.58614e-14
L 2 error (variable p) : 3.58907e-12
L inf error (variable p) : 4.26148e-12
----------------------------------------
Test_ChanFlow_LinNS_m8.xml
L 2 error (variable u) : 3.6131e-14
L inf error (variable u) : 4.10394e-13
L 2 error (variable v) : 2.12739e-14
L inf error (variable v) : 1.26191e-13
L 2 error (variable u) : 2.51232e-14
L inf error (variable u) : 3.04645e-13
L 2 error (variable v) : 1.2183e-14
L inf error (variable v) : 6.61758e-14
----------------------------------------
ChanStability.xml
L 2 error (variable u) : 0.0372818
......@@ -147,54 +147,54 @@ L 2 error (variable w) : 228.772
L inf error (variable w) : 0
----------------------------------------
Test_SinCos_LinNS_3DHom1D.xml
L 2 error (variable u) : 5.96531e-16
L inf error (variable u) : 5.99392e-15
L 2 error (variable v) : 7.1524e-16
L inf error (variable v) : 3.9968e-15
L 2 error (variable u) : 1.12396e-15
L inf error (variable u) : 5.01525e-15
L 2 error (variable v) : 1.35063e-15
L inf error (variable v) : 5.51422e-15
L 2 error (variable w) : 1.15463e-05
L inf error (variable w) : 1.11177e-05
----------------------------------------
Test_2DFlow_lineforcing_bcfromfile.xml
L 2 error (variable u) : 0.612387
L inf error (variable u) : 1.5
L 2 error (variable v) : 7.64016e-15
L inf error (variable v) : 3.2962e-14
L 2 error (variable v) : 4.19869e-14
L inf error (variable v) : 2.18145e-13
----------------------------------------
Test_ChanFlow_3DH1D_MVM.xml
L 2 error (variable u) : 7.28096e-16
L inf error (variable u) : 4.88498e-15
L 2 error (variable v) : 1.9485e-16
L inf error (variable v) : 3.6538e-16
L 2 error (variable u) : 3.61387e-16
L inf error (variable u) : 3.77476e-15
L 2 error (variable v) : 9.73834e-17
L inf error (variable v) : 3.84658e-16
L 2 error (variable w) : 0
L inf error (variable w) : 4.36612e-17
L 2 error (variable p) : 3.81525e-14
L inf error (variable p) : 8.79297e-14
L inf error (variable w) : 1.67056e-16
L 2 error (variable p) : 5.37526e-14
L inf error (variable p) : 9.76996e-14
----------------------------------------
Test_ChanFlow_3DH1D_FFT.xml
L 2 error (variable u) : 6.28245e-16
L inf error (variable u) : 4.27436e-15
L 2 error (variable v) : 1.13998e-16
L inf error (variable v) : 3.2573e-16
L 2 error (variable u) : 3.61998e-16
L inf error (variable u) : 1.88738e-15
L 2 error (variable v) : 0
L inf error (variable v) : 2.78215e-16
L 2 error (variable w) : 0
L inf error (variable w) : 0
L 2 error (variable p) : 4.27802e-14
L inf error (variable p) : 8.23785e-14
L 2 error (variable p) : 2.815e-14
L inf error (variable p) : 1.4988e-13
----------------------------------------
Test_ChanFlow_3DH2D_MVM.xml
L 2 error (variable u) : 0
L inf error (variable u) : 1.94533e-17
L 2 error (variable v) : 0
L inf error (variable v) : 5.55112e-17
L inf error (variable u) : 8.04776e-17
L 2 error (variable v) : 4.60587e-16
L inf error (variable v) : 1.08247e-15
L 2 error (variable w) : 0
L inf error (variable w) : 1.96598e-33
L 2 error (variable p) : 5.75425e-16
L inf error (variable p) : 8.98783e-16
L inf error (variable w) : 2.42788e-17
L 2 error (variable p) : 1.75622e-15
L inf error (variable p) : 3.58664e-15
----------------------------------------
Test_ChanFlow_3DH2D_FFT.xml
L 2 error (variable u) : 0
L inf error (variable u) : 0
L 2 error (variable v) : 0
L inf error (variable v) : 5.55112e-17
L inf error (variable v) : 1.11022e-16
L 2 error (variable w) : 0
L inf error (variable w) : 0
L 2 error (variable p) : 0
......@@ -211,12 +211,12 @@ L 2 error (variable p) : 8.55847e-05
L inf error (variable p) : 0.000513632
----------------------------------------
Test_ChanFlow_m3_SKS.xml
L 2 error (variable u) : 1.84051e-16
L inf error (variable u) : 1.44329e-15
L 2 error (variable u) : 4.62106e-16
L inf error (variable u) : 4.27436e-15
L 2 error (variable v) : 0
L inf error (variable v) : 1.58395e-16
L 2 error (variable p) : 1.1153e-14
L inf error (variable p) : 2.68674e-14
L inf error (variable v) : 2.79367e-16
L 2 error (variable p) : 1.10533e-14
L inf error (variable p) : 6.68354e-14
----------------------------------------
Test_KovaFlow_3DH1D_P5_20modes_SKS_MVM.xml
L 2 error (variable u) : 4.2366e-05
......@@ -230,160 +230,160 @@ L inf error (variable p) : 0.000565207
----------------------------------------
Test_Hex_channel_m3.xml
L 2 error (variable u) : 0
L inf error (variable u) : 1.83553e-16
L inf error (variable u) : 6.32744e-16
L 2 error (variable v) : 0
L inf error (variable v) : 2.01581e-16
L 2 error (variable w) : 3.40931e-16
L inf error (variable w) : 8.96505e-15
L 2 error (variable p) : 5.88883e-15
L inf error (variable p) : 5.26246e-14
L inf error (variable v) : 4.10249e-16
L 2 error (variable w) : 5.19692e-16
L inf error (variable w) : 5.32907e-15