Commit c4fddbaa authored by Spencer Sherwin's avatar Spencer Sherwin
Browse files

Updated Helmsolve calls to take initial guess in output vector. Also modifiied...

Updated Helmsolve calls to take initial guess in output vector. Also modifiied Iterative solver to check for zero solution before doing any preconditioning or matrix multiplies
parent 0bfd2a2f
......@@ -829,12 +829,12 @@ namespace Nektar
if(flags.isSet(eUseGlobal))
{
Vmath::Zero(contNcoeffs,outarray,1);
GlobalSolve(key,wsp,outarray,dirForcing);
}
else
{
Array<OneD,NekDouble> tmp(contNcoeffs,0.0);
Array<OneD,NekDouble> tmp(contNcoeffs);
LocalToGlobal(outarray,tmp);
GlobalSolve(key,wsp,tmp,dirForcing);
GlobalToLocal(tmp,outarray);
}
......
......@@ -96,6 +96,10 @@ namespace Nektar
const Array<OneD, const NekDouble> &inarray,
Array<OneD, NekDouble> &outarray) const;
inline void LocalToGlobal(
const Array<OneD, const NekDouble> &inarray,
Array<OneD, NekDouble> &outarray) const;
/// Assembles the global coefficients \f$\boldsymbol{\hat{u}}_g\f$
/// from the local coefficients \f$\boldsymbol{\hat{u}}_l\f$.
inline void Assemble();
......@@ -350,6 +354,13 @@ namespace Nektar
m_locToGloMap->GlobalToLocal(inarray,outarray);
}
inline void ContField2D::LocalToGlobal(
const Array<OneD, const NekDouble> &inarray,
Array<OneD,NekDouble> &outarray) const
{
m_locToGloMap->LocalToGlobal(inarray, outarray);
}
/**
* This operation is evaluated as:
* \f{tabbing}
......
......@@ -522,7 +522,8 @@ namespace Nektar
}
else
{
Array<OneD,NekDouble> tmp(contNcoeffs, 0.0);
Array<OneD,NekDouble> tmp(contNcoeffs);
LocalToGlobal(outarray,tmp);
GlobalSolve(key,wsp,tmp,dirForcing);
GlobalToLocal(tmp,outarray);
}
......
......@@ -199,6 +199,12 @@ namespace Nektar
NekVector<NekDouble> V_GlobHomBndTmp(nGlobHomBndDofs,0.0);
//zero GlobHomBnd so that we ensure we are solving for
//full problem rather than perturbation from initial
//condition in this case
Vmath::Zero(nGlobHomBndDofs,tmp = out+nDirBndDofs,1);
if(nGlobHomBndDofs)
{
if(nIntDofs || ((nDirBndDofs) && (!dirForcCalculated)
......
......@@ -50,10 +50,11 @@ namespace Nektar
const GlobalLinSysKey &pKey,
const boost::weak_ptr<ExpList> &pExpList,
const boost::shared_ptr<AssemblyMap>
&pLocToGloMap)
&pLocToGloMap)
: GlobalLinSys(pKey, pExpList, pLocToGloMap),
m_totalIterations(0),
m_useProjection(false),
m_rhs_magnitude(NekConstants::kNekUnsetDouble),
m_numPrevSols(0)
{
LibUtilities::SessionReaderSharedPtr vSession
......@@ -238,19 +239,19 @@ namespace Nektar
}
}
/**
* 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)
const int nGlobal,
const Array<OneD,const NekDouble> &in,
const int nDir)
{
// Get the communicator for performing data exchanges
LibUtilities::CommSharedPtr vComm
= m_expList.lock()->GetComm()->GetRowComm();
= m_expList.lock()->GetComm()->GetRowComm();
// Get vector sizes
int nNonDir = nGlobal - nDir;
......@@ -273,16 +274,16 @@ namespace Nektar
* Performs normalisation of input vector wrt A-norm.
*/
void GlobalLinSysIterative::UpdateKnownSolutions(
const int nGlobal,
const Array<OneD,const NekDouble> &newX,
const int nDir)
const int nGlobal,
const Array<OneD,const NekDouble> &newX,
const int nDir)
{
// Get vector sizes
int nNonDir = nGlobal - nDir;
// Get the communicator for performing data exchanges
LibUtilities::CommSharedPtr vComm
= m_expList.lock()->GetComm()->GetRowComm();
= m_expList.lock()->GetComm()->GetRowComm();
// Check the solution is non-zero
NekDouble solNorm = Vmath::Dot2(nNonDir,
......@@ -370,15 +371,15 @@ namespace Nektar
* The implemented algorithm uses a reduced-communication reordering of
* the standard PCG method (Demmel, Heath and Vorst, 1993)
*
* @param pInput Input vector of all DOFs.
* @param pInput Input residual of all DOFs.
* @param pOutput Solution vector of all DOFs.
*/
void GlobalLinSysIterative::DoConjugateGradient(
const int nGlobal,
const Array<OneD,const NekDouble> &pInput,
Array<OneD, NekDouble> &pOutput,
const AssemblyMapSharedPtr &plocToGloMap,
const int nDir)
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)
......@@ -393,7 +394,7 @@ namespace Nektar
// Get the communicator for performing data exchanges
LibUtilities::CommSharedPtr vComm
= m_expList.lock()->GetComm()->GetRowComm();
= m_expList.lock()->GetComm()->GetRowComm();
// Get vector sizes
int nNonDir = nGlobal - nDir;
......@@ -416,16 +417,48 @@ namespace Nektar
NekVector<NekDouble> q (nNonDir,q_A, eWrapper);
int k;
NekDouble alpha, beta, rho, rho_new, mu, eps, bb_inv, min_resid;
Array<OneD, NekDouble> vExchange(3);
NekDouble alpha, beta, rho, rho_new, mu, eps, min_resid;
Array<OneD, NekDouble> vExchange(3,0.0);
// Initialise with input initial guess.
// Copy initial residual from input
r = in;
// zero homogeneous out array ready for solution updates
// Should not be earlier in case input vector is same as
// output and above copy has been peformed
Vmath::Zero(nNonDir,tmp = pOutput + nDir,1);
// evaluate initial residual error for exit check
vExchange[2] = Vmath::Dot2(nNonDir,
r_A,
r_A,
m_map + nDir);
vComm->AllReduce(vExchange, Nektar::LibUtilities::ReduceSum);
eps = vExchange[2];
if(m_rhs_magnitude == NekConstants::kNekUnsetDouble)
{
m_rhs_magnitude = 1.0/vExchange[2];
}
// If input residual is less than tolerance skip solve.
if (eps < m_tolerance * m_tolerance * m_rhs_magnitude)
{
if (m_verbose && m_root)
{
cout << "CG iterations made = " << m_totalIterations
<< " using tolerance of " << m_tolerance
<< " (error = " << sqrt(vExchange[2]) << ")" << endl;
}
m_rhs_magnitude = NekConstants::kNekUnsetDouble;
return;
}
m_totalIterations = 1;
// If iteration is progressing calculate first iteration details
m_precon->DoPreconditioner(r_A, tmp = w_A + nDir);
v_DoMatrixMultiply(w_A, s_A);
k = 0;
......@@ -440,29 +473,14 @@ namespace Nektar
w_A + nDir,
m_map + nDir);
vExchange[2] = Vmath::Dot2(nNonDir,
r_A,
r_A,
m_map + nDir);
vComm->AllReduce(vExchange, Nektar::LibUtilities::ReduceSum);
m_totalIterations = 0;
// If input vector is zero, set zero output and skip solve.
if (vExchange[0] < NekConstants::kNekZeroTol)
{
Vmath::Zero(nNonDir, tmp = pOutput+nDir, 1);
return;
}
rho = vExchange[0];
mu = vExchange[1];
min_resid = m_rhs_magnitude;
beta = 0.0;
alpha = rho/mu;
eps = 0.0;
bb_inv = 1.0/vExchange[2];
min_resid = bb_inv;
// Continue until convergence
while (true)
......@@ -511,14 +529,17 @@ namespace Nektar
eps = vExchange[2];
// test if norm is within tolerance
if (eps*bb_inv < m_tolerance * m_tolerance)
if (eps < m_tolerance * m_tolerance * m_rhs_magnitude)
{
if (m_verbose && m_root)
{
cout << "CG iterations made = " << m_totalIterations
<< " using tolerance of " << m_tolerance
<< " (eps = " << sqrt(eps) << ")" << endl;
<< " (error = " << sqrt(eps/m_rhs_magnitude) << ")"
<< " (bb_inv = " << sqrt(1.0/m_rhs_magnitude) << ")"
<< endl;
}
m_rhs_magnitude = NekConstants::kNekUnsetDouble;
break;
}
min_resid = min(min_resid, eps);
......@@ -533,146 +554,16 @@ namespace Nektar
}
}
#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;
void GlobalLinSysIterative::Set_Rhs_Magnitude(const NekVector<NekDouble> &pIn)
{
// 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_z_old/vExchange[0];
// approximate solution
out = out + alpha*d;
// compute residual
r_new = r - 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
d = z_new + beta*d;
// Next step
r = r_new;
z = z_new;
k++;
ASSERTL0(k < 20000,
"Exceeded maximum number of iterations (20000)");
}
Array<OneD, NekDouble> vExchange(1);
vExchange[0] = Vmath::Dot(pIn.GetDimension(),&pIn[0],&pIn[0]);
m_expList.lock()->GetComm()->GetRowComm()->AllReduce(vExchange, Nektar::LibUtilities::ReduceSum);
m_rhs_magnitude = (vExchange[0] > 1e-6)? vExchange[0]:1.0;
}
#endif
}
}
......@@ -68,11 +68,14 @@ namespace Nektar
/// Tolerance of iterative solver.
NekDouble m_tolerance;
/// dot product of rhs to normalise stopping criterion
NekDouble m_rhs_magnitude;
PreconditionerSharedPtr m_precon;
MultiRegions::PreconditionerType m_precontype;
int m_totalIterations;
/// Whether to apply projection technique
......@@ -106,6 +109,8 @@ namespace Nektar
const int pNumDir);
void Set_Rhs_Magnitude(const NekVector<NekDouble> &pIn);
private:
void printArray(
......
......@@ -100,8 +100,8 @@ namespace Nektar
*
* @param pInput RHS of linear system, \f$b\f$.
* @param pOutput On input, values of dirichlet degrees
* of freedom. On output, the solution
* \f$x\f$.
* of freedom with initial guess on other values.
* On output, the solution \f$x\f$.
* @param pLocToGloMap Local to global mapping.
* @param pDirForcing Precalculated Dirichlet forcing.
*/
......@@ -155,8 +155,10 @@ namespace Nektar
}
if (vCG)
{
SolveLinearSystem(
nGlobDofs, tmp, pOutput, pLocToGloMap, nDirDofs);
Array<OneD, NekDouble> out(nGlobDofs,0.0);
// solve for perturbation from intiial guess in pOutput
SolveLinearSystem(nGlobDofs, tmp, out, pLocToGloMap, nDirDofs);
Vmath::Vadd(nGlobDofs,out,1,pOutput,1,out,1);
}
else
{
......
......@@ -210,6 +210,9 @@ namespace Nektar
NekVector<NekDouble> V_GlobHomBndTmp(nGlobHomBndDofs,0.0);
// set up normalisation factor for right hand side
Set_Rhs_Magnitude(F_GlobBnd);
if(nGlobHomBndDofs)
{
if(pLocToGloMap->GetPreconType() != MultiRegions::eLowEnergy)
......@@ -237,7 +240,7 @@ namespace Nektar
DNekScalBlkMat &BinvD = *m_BinvD;
V_LocBnd = BinvD*F_Int;
}
pLocToGloMap->AssembleBnd(V_LocBnd,V_GlobHomBndTmp,
nDirBndDofs);
F_HomBnd = F_HomBnd - V_GlobHomBndTmp;
......@@ -314,14 +317,17 @@ namespace Nektar
// solve boundary system
if(atLastLevel)
{
Array<OneD, NekDouble> offsetarray;
Array<OneD, NekDouble> pert(nGlobBndDofs,0.0);
Timer t;
t.Start();
// Solve for difference from initial solution given inout;
SolveLinearSystem(nGlobBndDofs, F, out, pLocToGloMap, nDirBndDofs);
SolveLinearSystem(nGlobBndDofs, F, pert, pLocToGloMap, nDirBndDofs);
// Add back initial conditions onto difference
Vmath::Vadd(nGlobHomBndDofs,&out[nDirBndDofs],1,
&pert[nDirBndDofs],1,&out[nDirBndDofs],1);
t.Stop();
//transform back to original basis
......
......@@ -87,18 +87,15 @@ namespace Nektar
ASSERTL0(i != (int) LibUtilities::SIZE_TimeIntegrationMethod, "Invalid time integration type.");
m_session->MatchSolverInfo("SpectralVanishingViscosity","True",m_useSpecVanVisc,false);
if(m_useSpecVanVisc)
{
m_session->LoadParameter("SVVCutoffRatio",m_sVVCutoffRatio,0.75);
m_session->LoadParameter("SVVDiffCoeff",m_sVVDiffCoeff,0.1);
}
m_session->LoadParameter("SVVCutoffRatio",m_sVVCutoffRatio,0.75);
m_session->LoadParameter("SVVDiffCoeff",m_sVVDiffCoeff,0.1);
m_session->MatchSolverInfo("SpectralVanishingViscosityHomo1D","True",m_useHomo1DSpecVanVisc,false);
if(m_HomogeneousType == eHomogeneous1D)
{
ASSERTL0(m_nConvectiveFields > 2,"Expect to have three velcoity fields with homogenous expansion");
m_session->MatchSolverInfo("SpectralVanishingViscosityHomo1D","True",m_useHomo1DSpecVanVisc,false);
if(m_useHomo1DSpecVanVisc == false)
{
......
Markdown is supported
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