Nektar++
All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
NekNonlinSysNewton.cpp
Go to the documentation of this file.
1///////////////////////////////////////////////////////////////////////////////
2//
3// File: NekNonlinSysNewton.cpp
4//
5// For more information, please see: http://www.nektar.info
6//
7// The MIT License
8//
9// Copyright (c) 2006 Division of Applied Mathematics, Brown University (USA),
10// Department of Aeronautics, Imperial College London (UK), and Scientific
11// Computing and Imaging Institute, University of Utah (USA).
12//
13// License for the specific language governing rights and limitations under
14// Permission is hereby granted, free of charge, to any person obtaining a
15// copy of this software and associated documentation files (the "Software"),
16// to deal in the Software without restriction, including without limitation
17// the rights to use, copy, modify, merge, publish, distribute, sublicense,
18// and/or sell copies of the Software, and to permit persons to whom the
19// Software is furnished to do so, subject to the following conditions:
20//
21// The above copyright notice and this permission notice shall be included
22// in all copies or substantial portions of the Software.
23//
24// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
25// OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
26// FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
27// THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
28// LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
29// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
30// DEALINGS IN THE SOFTWARE.
31//
32// Description: NekNonlinSysNewton definition
33//
34///////////////////////////////////////////////////////////////////////////////
35
38
39using namespace std;
40
41namespace Nektar
42{
43namespace LibUtilities
44{
45/**
46 * @class NekNonlinSysNewton
47 *
48 * Solves a nonlinear system using iterative methods.
49 */
52 "Newton", NekNonlinSysNewton::create, "NekNonlinSysNewton solver.");
53
56 const LibUtilities::CommSharedPtr &vRowComm, const int nscale,
57 const NekSysKey &pKey)
58 : NekNonlinSys(pSession, vRowComm, nscale, pKey)
59{
60}
61
63{
67}
68
70{
71}
72
73/**
74 *
75 **/
77 const int nGlobal, const Array<OneD, const NekDouble> &pInput,
78 Array<OneD, NekDouble> &pOutput, const int nDir, const NekDouble tol,
79 const NekDouble factor)
80{
81 boost::ignore_unused(factor);
82
83 v_SetupNekNonlinSystem(nGlobal, pInput, pInput, nDir);
84
85 m_Solution = pOutput;
86 Vmath::Vcopy(nGlobal - nDir, pInput, 1, m_Solution, 1);
87
88 int ntotal = nGlobal - nDir;
90 m_converged = false;
91
92 NekDouble resnormOld = 0.0;
93 int NttlNonlinIte = 0;
94 for (; NttlNonlinIte < m_maxiter; ++NttlNonlinIte)
95 {
97
98 m_converged = v_ConvergenceCheck(NttlNonlinIte, m_Residual, tol);
99 if (m_converged)
100 {
101 break;
102 }
103
104 NekDouble LinSysRelativeIteTol =
105 CalcInexactNewtonForcing(NttlNonlinIte, resnormOld, m_SysResNorm);
106 NekDouble LinSysTol = LinSysRelativeIteTol * sqrt(m_SysResNorm);
107 resnormOld = m_SysResNorm;
108
109 int ntmpLinSysIts =
110 m_linsol->SolveSystem(ntotal, m_Residual, m_DeltSltn, 0, LinSysTol);
111 m_NtotLinSysIts += ntmpLinSysIts;
112
113 Vmath::Vsub(ntotal, m_Solution, 1, m_DeltSltn, 1, m_Solution, 1);
114 }
115
116 if (((!m_converged) || m_verbose) && m_root && m_FlagWarnings)
117 {
118 int nwidthcolm = 11;
119
121 " # Nonlinear solver not converge in DoImplicitSolve");
122 cout << right << scientific << setw(nwidthcolm)
123 << setprecision(nwidthcolm - 6)
124 << " * Newton-Its converged (RES=" << sqrt(m_SysResNorm)
125 << " Res/(DtRHS): " << sqrt(m_SysResNorm / m_SysResNorm0)
126 << " with " << setw(3) << NttlNonlinIte << " Non-Its)" << endl;
127 }
128
129 return NttlNonlinIte;
130}
131
133 const int nIteration, const Array<OneD, const NekDouble> &Residual,
134 const NekDouble tol)
135{
136 bool converged = false;
137 NekDouble resratio = 1.0;
138 int ntotal = Residual.size();
139
140 m_SysResNorm = Vmath::Dot(ntotal, Residual, Residual);
142
143 if (nIteration == 0)
144 {
146 resratio = 1.0;
147 }
148 else
149 {
150 resratio = m_SysResNorm / m_SysResNorm0;
151 }
152
154 m_SysResNorm < tol * tol)
155 {
156 converged = true;
157 }
158
159 return converged;
160}
161
163 const int &k, const NekDouble &resnormOld, const NekDouble &resnorm)
164{
165 if (k == 0 || !m_InexactNewtonForcing)
166 {
168 }
169 else
170 {
171 NekDouble tmpForc =
172 m_forcingGamma * pow((resnorm / resnormOld), m_forcingAlpha);
173 return max(min(m_LinSysRelativeTolInNonlin, tmpForc), 1.0E-6);
174 }
175}
176
178 const int nGlobal, const Array<OneD, const NekDouble> &pInput,
179 const Array<OneD, const NekDouble> &pSource, const int nDir)
180{
181 boost::ignore_unused(pInput);
182
183 ASSERTL0(0 == nDir, "0 != nDir not tested");
184 ASSERTL0(m_SysDimen == nGlobal, "m_SysDimen!=nGlobal");
185
186 m_SourceVec = pSource;
187
188 m_linsol->SetSysOperators(m_operator);
189}
190
191} // namespace LibUtilities
192} // namespace Nektar
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:215
#define WARNINGL0(condition, msg)
Definition: ErrorUtil.hpp:222
tKey RegisterCreatorFunction(tKey idKey, CreatorFunction classCreator, std::string pDesc="")
Register a class with the factory.
Definition: NekFactory.hpp:198
Array< OneD, NekDouble > m_Residual
Definition: NekNonlinSys.h:133
Array< OneD, NekDouble > m_DeltSltn
Definition: NekNonlinSys.h:134
Array< OneD, NekDouble > m_SourceVec
Definition: NekNonlinSys.h:135
NekLinSysIterSharedPtr m_linsol
Definition: NekNonlinSys.h:122
Array< OneD, NekDouble > m_Solution
Definition: NekNonlinSys.h:132
static NekNonlinSysSharedPtr create(const LibUtilities::SessionReaderSharedPtr &pSession, const LibUtilities::CommSharedPtr &vRowComm, const int nDimen, const NekSysKey &pKey)
virtual void v_SetupNekNonlinSystem(const int nGlobal, const Array< OneD, const NekDouble > &pInput, const Array< OneD, const NekDouble > &pSource, const int nDir) override
NekDouble CalcInexactNewtonForcing(const int &k, const NekDouble &resnormOld, const NekDouble &resnorm)
NekNonlinSysNewton(const LibUtilities::SessionReaderSharedPtr &pSession, const LibUtilities::CommSharedPtr &vRowComm, const int nscale, const NekSysKey &pKey)
virtual int v_SolveSystem(const int nGlobal, const Array< OneD, const NekDouble > &pInput, Array< OneD, NekDouble > &pOutput, const int nDir, const NekDouble tol, const NekDouble factor) override
virtual bool v_ConvergenceCheck(const int nIteration, const Array< OneD, const NekDouble > &Residual, const NekDouble tol) override
bool m_root
Root if parallel.
Definition: NekSys.h:306
LibUtilities::CommSharedPtr m_rowComm
Communicate.
Definition: NekSys.h:302
virtual void v_InitObject()
Definition: NekSys.h:315
bool m_verbose
Verbose.
Definition: NekSys.h:308
NekSysOperators m_operator
Operators.
Definition: NekSys.h:311
int m_SysDimen
The dimension of the system.
Definition: NekSys.h:313
bool m_converged
Whether the iteration has been converged.
Definition: NekSys.h:304
int m_maxiter
Maximum iterations.
Definition: NekSys.h:298
void DoNekSysResEval(InArrayType &inarray, OutArrayType &outarray, const bool &flag=false) const
Definition: NekSys.h:136
NekNonlinSysFactory & GetNekNonlinSysFactory()
std::shared_ptr< SessionReader > SessionReaderSharedPtr
std::shared_ptr< Comm > CommSharedPtr
Pointer to a Communicator object.
Definition: Comm.h:57
The above copyright notice and this permission notice shall be included.
Definition: CoupledSolver.h:2
double NekDouble
T Dot(int n, const T *w, const T *x)
dot (vector times vector): z = w*x
Definition: Vmath.cpp:1095
void Vcopy(int n, const T *x, const int incx, T *y, const int incy)
Definition: Vmath.cpp:1191
void Vsub(int n, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
Subtract vector z = x-y.
Definition: Vmath.cpp:414
scalarT< T > sqrt(scalarT< T > in)
Definition: scalar.hpp:294