Nektar++
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
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  */
47  : Driver(pSession)
48  {
49  m_session->LoadParameter("IO_InfoSteps", m_infosteps, 1);
50  };
51 
52 
53  /**
54  *
55  */
57  {
58  };
59 
60 
61  /**
62  *
63  */
64  void DriverArnoldi::v_InitObject(ostream &out)
65  {
67  m_session->MatchSolverInfo("SolverType","VelocityCorrectionScheme",m_timeSteppingAlgorithm, false);
68 
70  {
71  m_period = m_session->GetParameter("TimeStep")* m_session->GetParameter("NumSteps");
72  m_nfields = m_equ[0]->UpdateFields().num_elements() - 1;
73 
74  if(m_session->DefinesSolverInfo("ModeType") &&
75  (m_session->GetSolverInfo("ModeType")=="SingleMode"|| m_session->GetSolverInfo("ModeType")=="HalfMode") )
76  {
77  for(int i = 0; i < m_nfields; ++i)
78  {
79  m_equ[0]->UpdateFields()[i]->SetWaveSpace(true);
80  }
81  }
82 
83  }
84  else
85  {
86  m_period = 1.0;
87  ASSERTL0(m_session->DefinesFunction("BodyForce"),"A BodyForce section needs to be defined for this solver type");
88  m_nfields = m_equ[0]->UpdateFields().num_elements();
89  }
90 
91  m_session->LoadParameter("kdim", m_kdim, 16);
92  m_session->LoadParameter("nvec", m_nvec, 2);
93  m_session->LoadParameter("nits", m_nits, 500);
94  m_session->LoadParameter("evtol", m_evtol, 1e-06);
95 
96  m_session->LoadParameter("realShift", m_realShift, 0.0);
97  m_equ[0]->SetLambda(m_realShift);
98 
99  m_session->LoadParameter("imagShift", m_imagShift, 0.0);
100  }
101 
102  void DriverArnoldi::ArnoldiSummary(std::ostream &out)
103  {
104  if (m_comm->GetRank() == 0)
105  {
106  if(m_session->DefinesSolverInfo("SingleMode"))
107  {
108  out << "\tSingle Fourier mode : true " << endl;
109  ASSERTL0(m_session->DefinesSolverInfo("Homogeneous"),
110  "Expected a homogeneous expansion to be defined "
111  "with single mode");
112  }
113  else
114  {
115  out << "\tSingle Fourier mode : false " << endl;
116  }
117  if(m_session->DefinesSolverInfo("BetaZero"))
118  {
119  out << "\tBeta set to Zero : true (overrides LHom)"
120  << endl;
121  }
122  else
123  {
124  out << "\tBeta set to Zero : false " << endl;
125  }
126 
128  {
129  out << "\tEvolution operator : "
130  << m_session->GetSolverInfo("EvolutionOperator")
131  << endl;
132  }
133  else
134  {
135  out << "\tShift (Real,Imag) : " << m_realShift
136  << "," << m_imagShift << endl;
137  }
138  out << "\tKrylov-space dimension : " << m_kdim << endl;
139  out << "\tNumber of vectors : " << m_nvec << endl;
140  out << "\tMax iterations : " << m_nits << endl;
141  out << "\tEigenvalue tolerance : " << m_evtol << endl;
142  out << "======================================================"
143  << endl;
144  }
145  }
146 
147  /**
148  * Copy Arnoldi array to field variables which depend from
149  * either the m_fields or m_forces
150  */
151  void DriverArnoldi::CopyArnoldiArrayToField(Array<OneD, NekDouble> &array)
152  {
153 
154  Array<OneD, MultiRegions::ExpListSharedPtr>& fields = m_equ[0]->UpdateFields();
155  int nq = fields[0]->GetNcoeffs();
156 
157  for (int k = 0; k < m_nfields; ++k)
158  {
159  Vmath::Vcopy(nq, &array[k*nq], 1, &fields[k]->UpdateCoeffs()[0], 1);
160  fields[k]->SetPhysState(false);
161  }
162  };
163 
164  /**
165  * Copy field variables which depend from either the m_fields
166  * or m_forces array the Arnoldi array
167  */
168  void DriverArnoldi::CopyFieldToArnoldiArray(Array<OneD, NekDouble> &array)
169  {
170 
171  Array<OneD, MultiRegions::ExpListSharedPtr> fields;
172 
174  {
175  //This matters for the Adaptive SFD method
176  //because m_equ[1] is the nonlinear problem with non homogeneous BCs.
177  fields = m_equ[0]->UpdateFields();
178  }
179  else
180  {
181  fields = m_equ[m_nequ-1]->UpdateFields();
182  }
183 
184  for (int k = 0; k < m_nfields; ++k)
185  {
186  int nq = fields[0]->GetNcoeffs();
187  Vmath::Vcopy(nq, &fields[k]->GetCoeffs()[0], 1, &array[k*nq], 1);
188  fields[k]->SetPhysState(false);
189 
190  }
191  };
192 
193 
194  /**
195  * Initialisation for the transient growth
196  */
198  {
199  Array<OneD, MultiRegions::ExpListSharedPtr> fields;
200 
202  {
203  fields = m_equ[0]->UpdateFields();
204  int nq = fields[0]->GetNcoeffs();
205 
206 
207  for (int k=0 ; k < m_nfields; ++k)
208  {
209  Vmath::Vcopy(nq, &fields[k]->GetCoeffs()[0], 1,&m_equ[1]->UpdateFields()[k]->UpdateCoeffs()[0], 1);
210 
211  }
212  }
213  else
214  {
215  ASSERTL0(false,"Transient Growth non available for Coupled Solver");
216 
217  }
218  };
219 
220  void DriverArnoldi::WriteFld(std::string file, std::vector<Array<OneD, NekDouble> > coeffs)
221  {
222 
223  std::vector<std::string> variables(m_nfields);
224 
225  ASSERTL1(coeffs.size() >= m_nfields, "coeffs is not of the correct length");
226  for(int i = 0; i < m_nfields; ++i)
227  {
228  variables[i] = m_equ[0]->GetVariable(i);
229  }
230 
231  m_equ[0]->WriteFld(file,m_equ[0]->UpdateFields()[0], coeffs, variables);
232  }
233 
234 
235  void DriverArnoldi::WriteFld(std::string file, Array<OneD, NekDouble> coeffs)
236  {
237 
238  std::vector<std::string> variables(m_nfields);
239  std::vector<Array<OneD, NekDouble> > fieldcoeffs(m_nfields);
240 
241  int ncoeffs = m_equ[0]->UpdateFields()[0]->GetNcoeffs();
242  ASSERTL1(coeffs.num_elements() >= ncoeffs*m_nfields,"coeffs is not of sufficient size");
243 
244  for(int i = 0; i < m_nfields; ++i)
245  {
246  variables[i] = m_equ[0]->GetVariable(i);
247  fieldcoeffs[i] = coeffs + i*ncoeffs;
248  }
249 
250  m_equ[0]->WriteFld(file,m_equ[0]->UpdateFields()[0], fieldcoeffs, variables);
251  }
252 
253  void DriverArnoldi::WriteEvs(ostream &evlout, const int i, const NekDouble re_ev, const NekDouble im_ev, NekDouble resid)
254  {
256  {
257  NekDouble abs_ev = hypot (re_ev, im_ev);
258  NekDouble ang_ev = atan2 (im_ev, re_ev);
259 
260  evlout << setw(2) << i
261  << setw(12) << abs_ev
262  << setw(12) << ang_ev
263  << setw(12) << log (abs_ev) / m_period
264  << setw(12) << ang_ev / m_period;
265 
266  if(resid != NekConstants::kNekUnsetDouble)
267  {
268  evlout << setw(12) << resid;
269  }
270  evlout << endl;
271  }
272  else
273  {
274  NekDouble invmag = 1.0/(re_ev*re_ev + im_ev*im_ev);
275 
276  evlout << setw(2) << i
277  << setw(14) << re_ev
278  << setw(14) << im_ev
279  << setw(14) << -re_ev*invmag + m_realShift
280  << setw(14) << im_ev*invmag;
281 
282  if(resid != NekConstants::kNekUnsetDouble)
283  {
284  evlout << setw(12) << resid;
285  }
286  evlout << endl;
287  }
288  }
289  }
290 }