Commit 6e95fa08 authored by Chris Cantwell's avatar Chris Cantwell
Browse files

Merge branch 'fix/IterICs' of localhost:nektar

parents 022709f1 6a7f4185
......@@ -890,6 +890,32 @@ namespace Nektar
}
Vmath::Vcopy(m_numGlobalBndCoeffs-offset, tmp.get()+offset, 1, global.get(), 1);
}
void AssemblyMap::LocalBndToGlobal(
const NekVector<NekDouble>& loc,
NekVector<NekDouble>& global) const
{
LocalBndToGlobal(loc.GetPtr(), global.GetPtr());
}
void AssemblyMap::LocalBndToGlobal(
const Array<OneD, const NekDouble>& loc,
Array<OneD,NekDouble>& global) const
{
ASSERTL1(loc.num_elements() >= m_numLocalBndCoeffs,"Local vector is not of correct dimension");
ASSERTL1(global.num_elements() >= m_numGlobalBndCoeffs,"Global vector is not of correct dimension");
if(m_signChange)
{
Vmath::Scatr(m_numLocalBndCoeffs, m_localToGlobalBndSign.get(), loc.get(), m_localToGlobalBndMap.get(), global.get());
}
else
{
Vmath::Scatr(m_numLocalBndCoeffs, loc.get(), m_localToGlobalBndMap.get(), global.get());
}
}
void AssemblyMap::AssembleBnd(
const NekVector<NekDouble>& loc,
......
......@@ -197,11 +197,19 @@ namespace Nektar
NekVector<NekDouble>& global,
int offset) const;
MULTI_REGIONS_EXPORT void LocalBndToGlobal(
const NekVector<NekDouble>& loc,
NekVector<NekDouble>& global) const;
MULTI_REGIONS_EXPORT void LocalBndToGlobal(
const Array<OneD, const NekDouble>& loc,
Array<OneD,NekDouble>& global,
int offset) const;
MULTI_REGIONS_EXPORT void LocalBndToGlobal(
const Array<OneD, const NekDouble>& loc,
Array<OneD,NekDouble>& global) const;
MULTI_REGIONS_EXPORT void AssembleBnd(const NekVector<NekDouble>& loc,
NekVector<NekDouble>& global, int offset) const;
......
......@@ -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(eps/m_rhs_magnitude) << ")" << 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)
......@@ -510,15 +528,19 @@ namespace Nektar
mu = vExchange[1];
eps = vExchange[2];
m_totalIterations++;
// 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);
......@@ -529,150 +551,19 @@ namespace Nektar
rho = rho_new;
k++;
m_totalIterations++;
}
}
#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
{
......
......@@ -175,6 +175,7 @@ namespace Nektar
{
bool dirForcCalculated = (bool) dirForcing.num_elements();
bool atLastLevel = pLocToGloMap->AtLastLevel();
int scLevel = pLocToGloMap->GetStaticCondLevel();
int nGlobDofs = pLocToGloMap->GetNumGlobalCoeffs();
int nGlobBndDofs = pLocToGloMap->GetNumGlobalBndCoeffs();
......@@ -210,6 +211,13 @@ namespace Nektar
NekVector<NekDouble> V_GlobHomBndTmp(nGlobHomBndDofs,0.0);
// set up normalisation factor for right hand side on first SC level
if(scLevel == 0)
{
Set_Rhs_Magnitude(F_GlobBnd);
}
if(nGlobHomBndDofs)
{
if(pLocToGloMap->GetPreconType() != MultiRegions::eLowEnergy)
......@@ -237,7 +245,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;
......@@ -248,7 +256,6 @@ namespace Nektar
// those lower levels, whilst not sending anything to the
// other partitions, and includes them in the modified right
// hand side vector.
int scLevel = pLocToGloMap->GetStaticCondLevel();
int lcLevel = pLocToGloMap->GetLowestStaticCondLevel();
if(atLastLevel && scLevel < lcLevel)
{
......@@ -314,14 +321,15 @@ namespace Nektar
// solve boundary system
if(atLastLevel)
{
Array<OneD, NekDouble> offsetarray;
Array<OneD, NekDouble> pert(nGlobBndDofs,0.0);
NekVector<NekDouble> Pert(nGlobBndDofs,pert,eWrapper);
Timer t;
t.Start();
// Solve for difference from initial solution given inout;
SolveLinearSystem(nGlobBndDofs, F, out, pLocToGloMap, nDirBndDofs);
SolveLinearSystem(nGlobBndDofs, F, pert, pLocToGloMap, nDirBndDofs);
t.Stop();
//transform back to original basis
......@@ -329,12 +337,15 @@ namespace Nektar
{
DNekScalBlkMat &RT = *m_RTBlk;
pLocToGloMap->GlobalToLocalBnd(V_GlobHomBnd,V_LocBnd, nDirBndDofs);
pLocToGloMap->GlobalToLocalBnd(Pert,V_LocBnd);
V_LocBnd=RT*V_LocBnd;
pLocToGloMap->LocalBndToGlobal(V_LocBnd,V_GlobHomBnd, nDirBndDofs);
pLocToGloMap->LocalBndToGlobal(V_LocBnd,Pert);
}
// Add back initial conditions onto difference
Vmath::Vadd(nGlobHomBndDofs,&out[nDirBndDofs],1,
&pert[nDirBndDofs],1,&out[nDirBndDofs],1);
}
else
{
......
......@@ -89,13 +89,12 @@ namespace Nektar
m_session->MatchSolverInfo("SpectralVanishingViscosity","True",m_useSpecVanVisc,false);
m_session->LoadParameter("SVVCutoffRatio",m_sVVCutoffRatio,0.75);
m_session->LoadParameter("SVVDiffCoeff",m_sVVDiffCoeff,0.1);