Commit bf725ec4 authored by Chris Cantwell's avatar Chris Cantwell
Browse files

Rewrote conjugate gradient algorithm to do only one communication per iteration

Added divergence detection to conjugate gradient solver.


git-svn-id: https://gforge.sci.utah.edu/svn/nektar/trunk@4001 305cdda6-5ce1-45b3-a98d-dfc68c8b3305
parent 19664e9f
......@@ -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-09;
static const NekDouble kNekIterativeTol = 1e-06;
}
} //end of namespace
......
......@@ -69,10 +69,15 @@ namespace Nektar
/**
* 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.
* @param pInput Input vector of non-Dirichlet DOFs.
* @param pOutput Solution vector of non-Dirichlet DOFs.
* using an auxiliary function v_DoMatrixMultiply defined by the
* specific solver. Distributed math routines are used to support
* parallel execution of the solver.
*
* The implemented algorithm uses a reduced-communication reordering of
* the standard PCG method (Amin, Sadayappan, Gudavalli, IEEE 1994)
*
* @param pInput Input vector of all DOFs.
* @param pOutput Solution vector of all DOFs.
*/
void GlobalLinSysIterative::v_SolveLinearSystem(
const int nGlobal,
......@@ -84,47 +89,42 @@ namespace Nektar
// Check if preconditioner has been computed and compute if needed.
if (!m_precon)
{
//v_ComputePreconditioner();
v_UniqueMap();
m_precon = MemoryManager<Preconditioner>::AllocateSharedPtr(GetSharedThisPtr(),plocToGloMap);
v_UniqueMap();
m_precon = MemoryManager<Preconditioner>::AllocateSharedPtr(
GetSharedThisPtr(),plocToGloMap);
}
// Get the communicator for performing data exchanges
LibUtilities::CommSharedPtr vComm = m_expList.lock()->GetComm()->GetRowComm();
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> q_A (nGlobal, 0.0);
Array<OneD, NekDouble> y_A (nNonDir, 0.0);
Array<OneD, NekDouble> r_A (nNonDir, 0.0);
Array<OneD, NekDouble> r_new_A(nNonDir, 0.0);
Array<OneD, NekDouble> z_A (nNonDir, 0.0);
// 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> 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);
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> r (nNonDir,r_A, eWrapper);
NekVector<NekDouble> z (nNonDir,z_A, 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);
NekDouble alpha, beta, r_dot_d, b_dot_b, min_resid;
Array<OneD, NekDouble> vExchange(3);
// Initialise with zero as the initial guess.
r = in;
m_precon->DoPreconditioner(r_A,z_A);
//z = M * r;
d = z;
m_precon->DoPreconditioner(r_A,z_A);
p = z;
k = 0;
vExchange[0] = Vmath::Dot2(nNonDir, r_A, r_A, m_map + nDir);
......@@ -137,74 +137,73 @@ namespace Nektar
Vmath::Zero(nGlobal, pOutput, 1);
return;
}
r_dot_z_old = vExchange[1];
b_dot_b = vExchange[0];
r_dot_d = vExchange[1];
min_resid = b_dot_b;
// 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_z_old/vExchange[0];
// approximate solution
out = out + alpha*d;
ASSERTL0(k < 20000,
"Exceeded maximum number of iterations (20000)");
// compute residual
r_new = r - alpha*p;
ASSERTL0(vExchange[0] <= b_dot_b || k < 10,
"Conjugate gradient diverged. Tolerance too small?"
"Minimum residual achieved: "
+ boost::lexical_cast<std::string>(sqrt(min_resid)));
m_precon->DoPreconditioner(r_new_A,z_new_A);
// Perform the method-specific matrix-vector multiply operation.
v_DoMatrixMultiply(p_A, q_A);
// Apply preconditioner to new residual
//z_new = M * r_new;
// Apply preconditioner
m_precon->DoPreconditioner(q_A + nDir, y_A);
// beta
// <r_k, r_k>
vExchange[0] = Vmath::Dot2(nNonDir,
r_new_A,
z_new_A,
r_A,
r_A,
m_map + nDir);
// residual
// <p_k, q_k>
vExchange[1] = Vmath::Dot2(nNonDir,
r_new_A,
r_new_A,
p_A + nDir,
q_A + nDir,
m_map + nDir);
// perform exchange
// <q_k, y_k>
vExchange[2] = Vmath::Dot2(nNonDir,
q_A + nDir,
y_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];
min_resid = min(min_resid, vExchange[0]);
// test if norm is within tolerance
if (normsq < m_tolerance * m_tolerance)
if (vExchange[0] < m_tolerance * m_tolerance)
{
break;
}
// compute alpha
alpha = r_dot_d / vExchange[1];
beta = alpha * vExchange[2] / vExchange[1] - 1.0;
r_dot_d *= beta;
// Update residual vector
r = r - alpha * q;
z = z - alpha * y;
// Update solution
out = out + alpha * p;
// Compute new search direction
d = z_new + beta*d;
p = z + beta * p;
// Next step
r = r_new;
z = z_new;
k++;
ASSERTL0(k < 20000,
"Exceeded maximum number of iterations (20000)");
}
}
}
}
......@@ -160,6 +160,7 @@
<CONDITIONS>
<PARAMETERS>
<P> Lambda = 1 </P>
<P> IterativeSolverTolerance = 5e-08 </P>
</PARAMETERS>
<VARIABLES>
......
......@@ -24,6 +24,7 @@
<P> Kinvis = 1 </P>
<P> HomModesZ = 4 </P>
<P> LZ = 1.0 </P>
<P> IterativeSolverTolerance = 5e-08 </P>
</PARAMETERS>
<VARIABLES>
......
......@@ -25,6 +25,7 @@
<P> HomModesZ = 4 </P>
<P> LZ = 1.0 </P>
<P> PROC_Z = 2 </P>
<P> IterativeSolverTolerance = 5e-08 </P>
</PARAMETERS>
<VARIABLES>
......
......@@ -22,6 +22,7 @@
<P> IO_CheckSteps = 1000 </P>
<P> IO_InfoSteps = 1000 </P>
<P> Kinvis = 1 </P>
<P> IterativeSolverTolerance = 5e-08 </P>
</PARAMETERS>
<VARIABLES>
......
......@@ -163,6 +163,7 @@
<P> IO_CheckSteps = 10 </P>
<P> IO_InfoSteps = 10 </P>
<P> Kinvis = 1 </P>
<P> IterativeSolverTolerance = 2e-07 </P>
</PARAMETERS>
<VARIABLES>
......
......@@ -331,6 +331,7 @@
<P> IO_CheckSteps = 10 </P>
<P> IO_InfoSteps = 10 </P>
<P> Kinvis = 1 </P>
<P> IterativeSolverTolerance = 5e-07 </P>
</PARAMETERS>
<VARIABLES>
......
......@@ -655,8 +655,9 @@
<P> TimeStep = 0.001 </P>
<P> NumSteps = 10 </P>
<P> IO_CheckSteps = 100 </P>
<P> IO_InfoSteps = 10 </P>
<P> IO_InfoSteps = 1 </P>
<P> Kinvis = 1 </P>
<P> IterativeSolverTolerance = 1e-07 </P>
</PARAMETERS>
<VARIABLES>
......
IncNavierStokesSolver
==========================================
Test_ChanFlow_m3.xml
L 2 error (variable u) : 4.80863e-16
L inf error (variable u) : 4.32987e-15
L 2 error (variable u) : 1.06003e-16
L inf error (variable u) : 8.32667e-16
L 2 error (variable v) : 0
L inf error (variable v) : 2.54809e-16
L 2 error (variable p) : 8.40089e-15
L inf error (variable p) : 4.9738e-14
L inf error (variable v) : 1.69891e-16
L 2 error (variable p) : 1.3259e-14
L inf error (variable p) : 2.73115e-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) : 5.41638e-16
L 2 error (variable p) : 7.49767e-15
L inf error (variable p) : 4.90237e-14
L inf error (variable v) : 1.97885e-16
L 2 error (variable p) : 2.98061e-14
L inf error (variable p) : 5.26916e-14
----------------------------------------
Test_ChanFlow_m8.xml
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
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
----------------------------------------
Test_ChanFlow2D_bcsfromfiles.xml
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
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
----------------------------------------
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) : 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
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
----------------------------------------
Test_ChanFlow_LinNS_m8.xml
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
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
----------------------------------------
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) : 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 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 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) : 4.19869e-14
L inf error (variable v) : 2.18145e-13
L 2 error (variable v) : 7.64016e-15
L inf error (variable v) : 3.2962e-14
----------------------------------------
Test_ChanFlow_3DH1D_MVM.xml
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 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 w) : 0
L inf error (variable w) : 1.67056e-16
L 2 error (variable p) : 5.37526e-14
L inf error (variable p) : 9.76996e-14
L inf error (variable w) : 4.36612e-17
L 2 error (variable p) : 3.81525e-14
L inf error (variable p) : 8.79297e-14
----------------------------------------
Test_ChanFlow_3DH1D_FFT.xml
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 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 w) : 0
L inf error (variable w) : 0
L 2 error (variable p) : 2.815e-14
L inf error (variable p) : 1.4988e-13
L 2 error (variable p) : 4.27802e-14
L inf error (variable p) : 8.23785e-14
----------------------------------------
Test_ChanFlow_3DH2D_MVM.xml
L 2 error (variable u) : 0
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 inf error (variable u) : 1.94533e-17
L 2 error (variable v) : 0
L inf error (variable v) : 5.55112e-17
L 2 error (variable w) : 0
L inf error (variable w) : 2.42788e-17
L 2 error (variable p) : 1.75622e-15
L inf error (variable p) : 3.58664e-15
L inf error (variable w) : 1.96598e-33
L 2 error (variable p) : 5.75425e-16
L inf error (variable p) : 8.98783e-16
----------------------------------------
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) : 1.11022e-16
L inf error (variable v) : 5.55112e-17
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) : 4.62106e-16
L inf error (variable u) : 4.27436e-15
L 2 error (variable u) : 1.84051e-16
L inf error (variable u) : 1.44329e-15
L 2 error (variable v) : 0
L inf error (variable v) : 2.79367e-16
L 2 error (variable p) : 1.10533e-14
L inf error (variable p) : 6.68354e-14
L inf error (variable v) : 1.58395e-16
L 2 error (variable p) : 1.1153e-14
L inf error (variable p) : 2.68674e-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) : 6.32744e-16
L inf error (variable u) : 1.83553e-16
L 2 error (variable v) : 0
L inf error (variable v) : 4.10249e-16
L 2 error (variable w) : 5.19692e-16
L inf error (variable w) : 5.32907e-15
L 2 error (variable p) : 8.7989e-15
L inf error (variable p) : 4.75175e-14
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
----------------------------------------
Test_Hex_channel_m8.xml
L 2 error (variable u) : 5.84733e-16
L inf error (variable u) : 5.16381e-15
L 2 error (variable v) : 1.23584e-15
L inf error (variable v) : 5.65416e-15
L 2 error (variable w) : 3.55911e-15
L inf error (variable w) : 3.26073e-13
L 2 error (variable p) : 9.91311e-14
L inf error (variable p) : 1.12677e-12
L 2 error (variable u) : 4.23389e-16
L inf error (variable u) : 4.01915e-15
L 2 error (variable v) : 5.02545e-16
L inf error (variable v) : 3.20812e-15
L 2 error (variable w) : 2.05176e-15
L inf error (variable w) : 3.28015e-13
L 2 error (variable p) : 4.27991e-14
L inf error (variable p) : 6.70908e-13
----------------------------------------
Test_Tet_channel_m3.xml
L 2 error (variable u) : 0
L inf error (variable u) : 2.5194e-16
L inf error (variable u) : 7.05365e-17
L 2 error (variable v) : 0
L inf error (variable v) : 2.83686e-16
L inf error (variable v) : 1.2892e-16
L 2 error (variable w) : 0
L inf error (variable w) : 1.58207e-15
L 2 error (variable p) : 1.33607e-14
L inf error (variable p) : 3.84137e-14
L inf error (variable w) : 3.05311e-16
L 2 error (variable p) : 2.66018e-15
L inf error (variable p) : 1.70974e-14
----------------------------------------
Test_Tet_equitri.xml
L 2 error (variable u) : 3.63545e-11
L inf error (variable u) : 6.40291e-10
L 2 error (variable v) : 1.04004e-10
L inf error (variable v) : 1.56871e-09
L 2 error (variable w) : 7.75326e-09
L inf error (variable w) : 2.59194e-08
L 2 error (variable p) : 7.59863e-08
L inf error (variable p) : 3.2549e-07
L 2 error (variable u) : 2.96004e-10
L inf error (variable u) : 3.17249e-09
L 2 error (variable v) : 3.44962e-10
L inf error (variable v) : 3.00739e-09
L 2 error (variable w) : 1.2646e-08
L inf error (variable w) : 1.98356e-07
L 2 error (variable p) : 5.02338e-07
L inf error (variable p) : 2.58776e-06
----------------------------------------
Test_Prism_channel_m6.xml
L 2 error (variable u) : 0
L inf error (variable u) : 5.25265e-16
L inf error (variable u) : 7.71125e-16
L 2 error (variable v) : 0
L inf error (variable v) : 4.32855e-16
L 2 error (variable w) : 3.10451e-16
L inf error (variable w) : 7.85483e-15
L 2 error (variable p) : 7.38999e-15
L inf error (variable p) : 6.99441e-14
L inf error (variable v) : 9.02258e-16
L 2 error (variable w) : 4.97271e-16
L inf error (variable w) : 1.42664e-14
L 2 error (variable p) : 1.27331e-14
L inf error (variable p) : 9.75886e-14
----------------------------------------
Test_Hex_channel_m8_par.xml
L 2 error (variable u) : 1.03412e-11
L inf error (variable u) : 9.31369e-11
L 2 error (variable v) : 7.62679e-12
L inf error (variable v) : 7.31505e-11
L 2 error (variable w) : 8.14841e-11
L inf error (variable w) : 5.83095e-09
L 2 error (variable p) : 2.21238e-09
L inf error (variable p) : 9.42423e-08
L 2 error (variable u) : 2.59382e-09
L inf error (variable u) : 2.61441e-08
L 2 error (variable v) : 2.92001e-09
L inf error (variable v) : 5.88714e-08
L 2 error (variable w) : 1.86053e-08
L inf error (variable w) : 1.5756e-06
L 2 error (variable p) : 4.93257e-07
L inf error (variable p) : 5.37924e-05
----------------------------------------
Test_Tet_channel_m8_par.xml
L 2 error (variable u) : 1.04057e-11
L inf error (variable u) : 6.08492e-11
L 2 error (variable v) : 9.66528e-12
L inf error (variable v) : 7.06681e-11
L 2 error (variable w) : 1.09503e-10
L inf error (variable w) : 2.66938e-09
L 2 error (variable p) : 3.55555e-09
L inf error (variable p) : 5.6127e-08
L 2 error (variable u) : 5.00228e-09
L inf error (variable u) : 3.07381e-08
L 2 error (variable v) : 5.97896e-09
L inf error (variable v) : 4.76502e-08
L 2 error (variable w) : 6.55109e-08
L inf error (variable w) : 2.19381e-06
L 2 error (variable p) : 1.78493e-06
L inf error (variable p) : 7.82052e-05
----------------------------------------
Test_ChanFlow_m3_par.xml
L 2 error (variable u) : 8.42518e-14
L inf error (variable u) : 6.34659e-13
L 2 error (variable v) : 4.35393e-13
L inf error (variable v) : 6.83331e-13
L 2 error (variable p) : 6.56835e-11
L inf error (variable p) : 1.68574e-10
L 2 error (variable u) : 1.70764e-12
L inf error (variable u) : 8.91182e-12
L 2 error (variable v) : 6.46485e-12
L inf error (variable v) : 1.01549e-11