Nektar++
DriverArnoldi.cpp
Go to the documentation of this file.
1 ///////////////////////////////////////////////////////////////////////////////
2 //
3 // File DriverArnoldi.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: Base Driver class for the stability solver
33 //
34 ///////////////////////////////////////////////////////////////////////////////
35 #include <iomanip>
36 
38 
39 namespace Nektar
40 {
41 namespace SolverUtils
42 {
43 
44 /**
45  * Constructor
46  */
48  : Driver(pSession)
49 {
50  m_session->LoadParameter("IO_InfoSteps", m_infosteps, 1);
51 };
52 
53 
54 /**
55  * Destructor
56  */
58 {
59 };
60 
61 
62 /**
63  * Arnoldi driver initialisation
64  */
65 void DriverArnoldi::v_InitObject(ostream &out)
66 {
68  m_session->MatchSolverInfo("SolverType",
69  "VelocityCorrectionScheme",
71 
73  {
74  m_period = m_session->GetParameter("TimeStep")
75  * m_session->GetParameter("NumSteps");
76  m_nfields = m_equ[0]->UpdateFields().num_elements() - 1;
77 
78  }
79  else
80  {
81  m_period = 1.0;
82  m_nfields = m_equ[0]->UpdateFields().num_elements();
83  }
84 
85  if(m_session->DefinesSolverInfo("ModeType") &&
86  (boost::iequals(m_session->GetSolverInfo("ModeType"),
87  "SingleMode")||
88  boost::iequals(m_session->GetSolverInfo("ModeType"),
89  "HalfMode")))
90  {
91  for(int i = 0; i < m_nfields; ++i)
92  {
93  m_equ[0]->UpdateFields()[i]->SetWaveSpace(true);
94  }
95  }
96  m_negatedOp = m_equ[0]->v_NegatedOp();
97 
98  m_session->LoadParameter("kdim", m_kdim, 16);
99  m_session->LoadParameter("nvec", m_nvec, 2);
100  m_session->LoadParameter("nits", m_nits, 500);
101  m_session->LoadParameter("evtol", m_evtol, 1e-06);
102 
103  m_session->LoadParameter("realShift", m_realShift, 0.0);
104  m_equ[0]->SetLambda(m_realShift);
105 
106  m_session->LoadParameter("imagShift", m_imagShift, 0.0);
107 
108 }
109 
110 void DriverArnoldi::ArnoldiSummary(std::ostream &out)
111 {
112  if (m_comm->GetRank() == 0)
113  {
114  if(m_session->DefinesSolverInfo("ModeType") &&
115  boost::iequals(m_session->GetSolverInfo("ModeType"),
116  "SingleMode"))
117  {
118  out << "\tSingle Fourier mode : true " << endl;
119  ASSERTL0(m_session->DefinesSolverInfo("Homogeneous"),
120  "Expected a homogeneous expansion to be defined "
121  "with single mode");
122  }
123  else
124  {
125  out << "\tSingle Fourier mode : false " << endl;
126  }
127  if(m_session->DefinesSolverInfo("BetaZero"))
128  {
129  out << "\tBeta set to Zero : true (overrides LHom)"
130  << endl;
131  }
132  else
133  {
134  out << "\tBeta set to Zero : false " << endl;
135  }
136 
138  {
139  out << "\tEvolution operator : "
140  << m_session->GetSolverInfo("EvolutionOperator")
141  << endl;
142  }
143  else
144  {
145  out << "\tShift (Real,Imag) : " << m_realShift
146  << "," << m_imagShift << endl;
147  }
148  out << "\tKrylov-space dimension : " << m_kdim << endl;
149  out << "\tNumber of vectors : " << m_nvec << endl;
150  out << "\tMax iterations : " << m_nits << endl;
151  out << "\tEigenvalue tolerance : " << m_evtol << endl;
152  out << "======================================================"
153  << endl;
154  }
155 }
156 
157 /**
158  * Copy Arnoldi array to field variables which depend from
159  * either the m_fields or m_forces
160  */
162 {
163 
164  Array<OneD, MultiRegions::ExpListSharedPtr>& fields = m_equ[0]->UpdateFields();
165  int nq = fields[0]->GetNcoeffs();
166 
167  for (int k = 0; k < m_nfields; ++k)
168  {
169  Vmath::Vcopy(nq, &array[k*nq], 1, &fields[k]->UpdateCoeffs()[0], 1);
170  fields[k]->SetPhysState(false);
171  }
172 };
173 
174 /**
175  * Copy field variables which depend from either the m_fields
176  * or m_forces array the Arnoldi array
177  */
179 {
180 
182 
184  {
185  // This matters for the Adaptive SFD method because
186  // m_equ[1] is the nonlinear problem with non
187  // homogeneous BCs.
188  fields = m_equ[0]->UpdateFields();
189  }
190  else
191  {
192  fields = m_equ[m_nequ-1]->UpdateFields();
193  }
194 
195  for (int k = 0; k < m_nfields; ++k)
196  {
197  int nq = fields[0]->GetNcoeffs();
198  Vmath::Vcopy(nq, &fields[k]->GetCoeffs()[0], 1, &array[k*nq], 1);
199  fields[k]->SetPhysState(false);
200 
201  }
202 };
203 
204 
205 /**
206  * Initialisation for the transient growth
207  */
209 {
211 
213  {
214  fields = m_equ[0]->UpdateFields();
215  int nq = fields[0]->GetNcoeffs();
216 
217  for (int k=0 ; k < m_nfields; ++k)
218  {
219  Vmath::Vcopy(nq,
220  &fields[k]->GetCoeffs()[0], 1,
221  &m_equ[1]->UpdateFields()[k]->UpdateCoeffs()[0], 1);
222 
223  }
224  }
225  else
226  {
227  ASSERTL0(false,"Transient Growth non available for Coupled Solver");
228 
229  }
230 };
231 
232 void DriverArnoldi::WriteFld(std::string file, std::vector<Array<OneD, NekDouble> > coeffs)
233 {
234 
235  std::vector<std::string> variables(m_nfields);
236 
237  ASSERTL1(coeffs.size() >= m_nfields, "coeffs is not of the correct length");
238  for(int i = 0; i < m_nfields; ++i)
239  {
240  variables[i] = m_equ[0]->GetVariable(i);
241  }
242 
243  m_equ[0]->WriteFld(file,m_equ[0]->UpdateFields()[0], coeffs, variables);
244 }
245 
246 
247 void DriverArnoldi::WriteFld(std::string file, Array<OneD, NekDouble> coeffs)
248 {
249 
250  std::vector<std::string> variables(m_nfields);
251  std::vector<Array<OneD, NekDouble> > fieldcoeffs(m_nfields);
252 
253  int ncoeffs = m_equ[0]->UpdateFields()[0]->GetNcoeffs();
254  ASSERTL1(coeffs.num_elements() >= ncoeffs*m_nfields,"coeffs is not of sufficient size");
255 
256  for(int i = 0; i < m_nfields; ++i)
257  {
258  variables[i] = m_equ[0]->GetVariable(i);
259  fieldcoeffs[i] = coeffs + i*ncoeffs;
260  }
261 
262  m_equ[0]->WriteFld(file,m_equ[0]->UpdateFields()[0], fieldcoeffs, variables);
263 }
264 
266  ostream &evlout,
267  const int i,
268  const NekDouble re_ev,
269  const NekDouble im_ev,
270  NekDouble resid,
271  bool DumpInverse)
272 {
274  {
275  NekDouble abs_ev = hypot (re_ev, im_ev);
276  NekDouble ang_ev = atan2 (im_ev, re_ev);
277 
278  evlout << "EV: " << setw(2) << i
279  << setw(12) << abs_ev
280  << setw(12) << ang_ev
281  << setw(12) << log (abs_ev) / m_period
282  << setw(12) << ang_ev / m_period;
283 
284  if(resid != NekConstants::kNekUnsetDouble)
285  {
286  evlout << setw(12) << resid;
287  }
288  evlout << endl;
289  }
290  else
291  {
292  NekDouble invmag = 1.0/(re_ev*re_ev + im_ev*im_ev);
293  NekDouble sign;
294  if(m_negatedOp)
295  {
296  sign = -1.0;
297  }
298  else
299  {
300  sign = 1.0;
301  }
302 
303  evlout << "EV: " << setw(2) << i
304  << setw(14) << sign*re_ev
305  << setw(14) << sign*im_ev;
306 
307  if(DumpInverse)
308  {
309  evlout << setw(14) << sign*re_ev*invmag + m_realShift
310  << setw(14) << sign*im_ev*invmag + m_imagShift;
311  }
312 
313  if(resid != NekConstants::kNekUnsetDouble)
314  {
315  evlout << setw(12) << resid;
316  }
317  evlout << endl;
318  }
319 }
320 
321 }
322 }
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:161
DriverArnoldi(const LibUtilities::SessionReaderSharedPtr pSession)
Constructor.
void CopyFwdToAdj()
Copy the forward field to the adjoint system in transient growth calculations.
NekDouble m_evtol
Maxmum number of iterations.
Definition: DriverArnoldi.h:58
virtual SOLVER_UTILS_EXPORT void v_InitObject(ostream &out=cout)
Definition: Driver.cpp:90
void CopyArnoldiArrayToField(Array< OneD, NekDouble > &array)
Copy Arnoldi storage to fields.
#define sign(a, b)
return the sign(b)*a
Definition: Polylib.cpp:22
SOLVER_UTILS_EXPORT void ArnoldiSummary(std::ostream &out)
void CopyFieldToArnoldiArray(Array< OneD, NekDouble > &array)
Copy fields to Arnoldi storage.
boost::shared_ptr< SessionReader > SessionReaderSharedPtr
Definition: MeshPartition.h:50
void WriteFld(std::string file, std::vector< Array< OneD, NekDouble > > coeffs)
Write coefficients to file.
virtual ~DriverArnoldi()
Destructor.
int m_nvec
Dimension of Krylov subspace.
Definition: DriverArnoldi.h:56
enum EvolutionOperatorType m_EvolutionOperator
Evolution Operator.
Definition: Driver.h:100
double NekDouble
static const NekDouble kNekUnsetDouble
int m_nfields
interval to dump information if required.
Definition: DriverArnoldi.h:64
bool m_timeSteppingAlgorithm
Period of time stepping algorithm.
Definition: DriverArnoldi.h:60
void WriteEvs(ostream &evlout, const int k, const NekDouble real, const NekDouble imag, NekDouble resid=NekConstants::kNekUnsetDouble, bool DumpInverse=true)
Array< OneD, EquationSystemSharedPtr > m_equ
Equation system to solve.
Definition: Driver.h:94
NekDouble m_period
Tolerance of iteratiosn.
Definition: DriverArnoldi.h:59
int m_infosteps
underlying operator is time stepping
Definition: DriverArnoldi.h:62
LibUtilities::CommSharedPtr m_comm
Communication object.
Definition: Driver.h:85
Base class for the development of solvers.
Definition: Driver.h:65
#define ASSERTL1(condition, msg)
Assert Level 1 – Debugging which is used whether in FULLDEBUG or DEBUG compilation mode...
Definition: ErrorUtil.hpp:191
void Vcopy(int n, const T *x, const int incx, T *y, const int incy)
Definition: Vmath.cpp:1038
virtual void v_InitObject(ostream &out=cout)
LibUtilities::SessionReaderSharedPtr m_session
Session reader object.
Definition: Driver.h:88
int m_nequ
number of equations
Definition: Driver.h:97
int m_nits
Number of vectors to test.
Definition: DriverArnoldi.h:57