Commit 508f3779 authored by Pavel Burovskiy's avatar Pavel Burovskiy
Browse files

Adding projection method for speeding up successive linear solves

git-svn-id: https://gforge.sci.utah.edu/svn/nektar/trunk@4039 305cdda6-5ce1-45b3-a98d-dfc68c8b3305
parent c4810613
......@@ -43,7 +43,7 @@ namespace Nektar
/**
* @class GlobalLinSysIterative
*
* Solves a linear system using direct methods.
* Solves a linear system using iterative methods.
*/
/// Constructor for full direct matrix solve.
......@@ -52,7 +52,10 @@ namespace Nektar
const boost::weak_ptr<ExpList> &pExpList,
const boost::shared_ptr<AssemblyMap>
&pLocToGloMap)
: GlobalLinSys(pKey, pExpList, pLocToGloMap)
: GlobalLinSys(pKey, pExpList, pLocToGloMap),
m_totalIterations(0),
m_useProjection(false),
m_numPrevSols(0)
{
LibUtilities::SessionReaderSharedPtr vSession
= pExpList.lock()->GetSession();
......@@ -60,12 +63,267 @@ namespace Nektar
m_tolerance,
NekConstants::kNekIterativeTol);
std::string successiveRhs;
vSession->LoadSolverInfo("SuccessiveRHS", successiveRhs );
try
{
int solutionsToStore = boost::lexical_cast<int>(successiveRhs);
m_prevLinSol.set_capacity(solutionsToStore);
m_useProjection = true;
std::cout << "Using successive rhs projection with " << solutionsToStore << " solutions to be stored" << std::endl;
}
catch(...)
{
// lexical_cast will throw bad_lexical_cast if successiveRhs is not integer-valued
}
}
GlobalLinSysIterative::~GlobalLinSysIterative()
{
}
/**
*
*/
void GlobalLinSysIterative::v_SolveLinearSystem(
const int nGlobal,
const Array<OneD,const NekDouble> &pInput,
Array<OneD, NekDouble> &pOutput,
const AssemblyMapSharedPtr &plocToGloMap,
const int nDir)
{
// Get vector sizes
int nNonDir = nGlobal - nDir;
if (m_useProjection)
{
DoAconjugateProjection(nGlobal, pInput, pOutput, plocToGloMap, nDir);
if (0)
{
// check correctness: solve the same system with plain CG and compare
Array<OneD, NekDouble> cg_s (nGlobal, 0.0);
NekVector<NekDouble> cg (nNonDir, cg_s + nDir, eWrapper);
NekVector<NekDouble> x (nNonDir, pOutput + nDir, eWrapper);
DoConjugateGradient(nGlobal, pInput, cg_s, plocToGloMap, nDir);
cg -= x;
NekDouble norm = CalculateAnorm(nGlobal, cg_s, nDir);
std::cout << "norm of solutions difference is = " << norm << std::endl;
}
}
else
{
// applying plain Conjugate Gradient
DoConjugateGradient(nGlobal, pInput, pOutput, plocToGloMap, nDir);
}
//std::cout << "CG iterations made = " << m_totalIterations << std::endl << std::endl;
}
/**
* This method implements A-conjugate projection technique
* in order to speed up successive linear solves with
* right-hand sides arising from time-dependent discretisations.
* (P.F.Fischer, Comput. Methods Appl. Mech. Engrg. 163, 1998)
*
*/
void GlobalLinSysIterative::DoAconjugateProjection(
const int nGlobal,
const Array<OneD,const NekDouble> &pInput,
Array<OneD, NekDouble> &pOutput,
const AssemblyMapSharedPtr &plocToGloMap,
const int nDir)
{
// Get vector sizes
int nNonDir = nGlobal - nDir;
if (0 == m_numPrevSols)
{
// no previous solutions found, call CG
DoConjugateGradient(nGlobal, pInput, pOutput, plocToGloMap, nDir);
UpdateKnownSolutions(nGlobal, pOutput, nDir);
}
else
{
// Create NekVector wrappers for linear algebra operations
NekVector<NekDouble> b (nNonDir, pInput + nDir, eWrapper);
NekVector<NekDouble> x (nNonDir, pOutput + nDir, eWrapper);
// check the input vector (rhs) is not zero
NekDouble rhsNorm = Vmath::Dot(nNonDir,
pInput + nDir,
pInput + nDir);
if (rhsNorm < NekConstants::kNekZeroTol)
{
Array<OneD, NekDouble> tmp = pOutput+nDir;
Vmath::Zero(nNonDir, tmp, 1);
return;
}
// Allocate array storage
Array<OneD, NekDouble> px_s (nGlobal, 0.0);
Array<OneD, NekDouble> pb_s (nGlobal, 0.0);
Array<OneD, NekDouble> tmpAx_s (nGlobal, 0.0);
Array<OneD, NekDouble> tmpx_s (nGlobal, 0.0);
NekVector<NekDouble> pb (nNonDir, pb_s + nDir, eWrapper);
NekVector<NekDouble> px (nNonDir, px_s + nDir, eWrapper);
NekVector<NekDouble> tmpAx (nNonDir, tmpAx_s + nDir, eWrapper);
NekVector<NekDouble> tmpx (nNonDir, tmpx_s + nDir, eWrapper);
// notation follows the paper cited:
// \alpha_i = \tilda{x_i}^T b^n
// projected x, px = \sum \alpha_i \tilda{x_i}
for (int i = 0; i < m_prevLinSol.size(); i++)
{
NekDouble alphai = Vmath::Dot(nNonDir,
m_prevLinSol[i],
pInput + nDir);
if (alphai < NekConstants::kNekZeroTol)
{
continue;
}
NekVector<NekDouble> xi (nNonDir, m_prevLinSol[i], eWrapper);
px += alphai * xi;
}
// pb = b^n - A px
Vmath::Vcopy(nNonDir,
pInput.get() + nDir, 1,
pb_s.get() + nDir, 1);
v_DoMatrixMultiply(px_s, tmpAx_s);
pb -= tmpAx;
// solve the system with projected rhs
DoConjugateGradient(nGlobal, pb_s, tmpx_s, plocToGloMap, nDir);
// remainder solution + projection of previous solutions
x = tmpx + px;
// save the auxiliary solution to prev. known solutions
UpdateKnownSolutions(nGlobal, tmpx_s, nDir);
}
}
/**
* Calculating A-norm of an input vector,
* A-norm(x) := sqrt( < x, Ax > )
*/
NekDouble GlobalLinSysIterative::CalculateAnorm(
const int nGlobal,
const Array<OneD,const NekDouble> &in,
const int nDir)
{
// Get vector sizes
int nNonDir = nGlobal - nDir;
// Allocate array storage
Array<OneD, NekDouble> tmpAx_s (nGlobal, 0.0);
v_DoMatrixMultiply(in, tmpAx_s);
const NekDouble anorm_sq = Vmath::Dot(nNonDir,
in + nDir,
tmpAx_s + nDir);
return std::sqrt(anorm_sq);
}
/**
* Updates the storage of previously known solutions.
* Performs normalisation of input vector wrt A-norm.
*/
void GlobalLinSysIterative::UpdateKnownSolutions(
const int nGlobal,
const Array<OneD,const NekDouble> &newX,
const int nDir)
{
// Get vector sizes
int nNonDir = nGlobal - nDir;
// Check the solution is non-zero
NekDouble solNorm = Vmath::Dot(nNonDir,
newX + nDir,
newX + nDir);
if (solNorm < NekConstants::kNekZeroTol)
{
return;
}
// Allocate array storage
Array<OneD, NekDouble> tmpAx_s (nGlobal, 0.0);
Array<OneD, NekDouble> px_s (nGlobal, 0.0);
Array<OneD, NekDouble> tmp1, tmp2;
// Create NekVector wrappers for linear algebra operations
NekVector<NekDouble> px (nNonDir, px_s + nDir, eWrapper);
NekVector<NekDouble> tmpAx (nNonDir, tmpAx_s + nDir, eWrapper);
// calculating \tilda{x} - sum \alpha_i\tilda{x}_i
Vmath::Vcopy(nNonDir,
tmp1 = newX + nDir, 1,
tmp2 = px_s + nDir, 1);
if (m_prevLinSol.size() > 0)
{
v_DoMatrixMultiply(newX, tmpAx_s);
}
for (int i = 0; i < m_prevLinSol.size(); i++)
{
NekDouble alphai = Vmath::Dot(nNonDir,
m_prevLinSol[i],
tmpAx_s + nDir);
if (alphai < NekConstants::kNekZeroTol)
{
continue;
}
NekVector<NekDouble> xi (nNonDir, m_prevLinSol[i], eWrapper);
px -= alphai * xi;
}
// Some solutions generated by CG are identical zeros, see
// solutions generated for Test_Tet_equitri.xml (IncNavierStokesSolver).
// Not going to store identically zero solutions.
NekDouble anorm = CalculateAnorm(nGlobal, px_s, nDir);
if (anorm < NekConstants::kNekZeroTol)
{
return;
}
// normalisation of new solution
Vmath::Smul(nNonDir, 1.0/anorm, px_s.get() + nDir, 1, px_s.get() + nDir, 1);
// updating storage with non-Dirichlet-dof part of new solution vector
m_prevLinSol.push_back(px_s + nDir);
m_numPrevSols++;
}
/**
* Solve a global linear system using the conjugate gradient method.
* We solve only for the non-Dirichlet modes. The operator is evaluated
......@@ -79,7 +337,7 @@ namespace Nektar
* @param pInput Input vector of all DOFs.
* @param pOutput Solution vector of all DOFs.
*/
void GlobalLinSysIterative::v_SolveLinearSystem(
void GlobalLinSysIterative::DoConjugateGradient(
const int nGlobal,
const Array<OneD,const NekDouble> &pInput,
Array<OneD, NekDouble> &pOutput,
......@@ -224,6 +482,7 @@ namespace Nektar
rho = rho_new;
k++;
m_totalIterations++;
}
}
......
......@@ -40,6 +40,8 @@
#include <MultiRegions/Preconditioner.h>
#include <MultiRegions/AssemblyMap/AssemblyMapCG.h>
#include <boost/circular_buffer.hpp>
namespace Nektar
{
namespace MultiRegions
......@@ -72,9 +74,58 @@ namespace Nektar
PreconditionerSharedPtr m_precon;
MultiRegions::PreconditionerType m_precontype;
MultiRegions::PreconditionerType m_precontype;
int m_totalIterations;
/// Whether to apply projection technique
bool m_useProjection;
/// Storage for solutions to previous linear problems
boost::circular_buffer<Array<OneD, NekDouble> > m_prevLinSol;
/// Total counter of previous solutions
int m_numPrevSols;
/// A-conjugate projection technique
void DoAconjugateProjection(
const int pNumRows,
const Array<OneD,const NekDouble> &pInput,
Array<OneD, NekDouble> &pOutput,
const AssemblyMapSharedPtr &locToGloMap,
const int pNumDir);
/// Actual iterative solve
void DoConjugateGradient(
const int pNumRows,
const Array<OneD,const NekDouble> &pInput,
Array<OneD, NekDouble> &pOutput,
const AssemblyMapSharedPtr &locToGloMap,
const int pNumDir);
private:
void printArray(
const std::string& msg,
const Array<OneD, const NekDouble> &in,
const int len,
const int offset);
void UpdateKnownSolutions(
const int pGlobalBndDofs,
const Array<OneD,const NekDouble> &pSolution,
const int pNumDirBndDofs);
NekDouble CalculateAnorm(
const int nGlobal,
const Array<OneD,const NekDouble> &in,
const int nDir);
/// Solve the matrix system
virtual void v_SolveLinearSystem(
const int pNumRows,
......
......@@ -649,6 +649,7 @@
<I PROPERTY="Projection" VALUE="Galerkin" />
<I PROPERTY="TimeIntegrationMethod" VALUE="IMEXOrder1" />
<I PROPERTY="GlobalSysSoln" VALUE="IterativeStaticCond" />
<\I PROPERTY="SuccessiveRHS" VALUE="20"/>
</SOLVERINFO>
<PARAMETERS>
......
Supports Markdown
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