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// Permission is hereby granted, free of charge, to any person obtaining a
14// copy of this software and associated documentation files (the "Software"),
15// to deal in the Software without restriction, including without limitation
16// the rights to use, copy, modify, merge, publish, distribute, sublicense,
17// and/or sell copies of the Software, and to permit persons to whom the
18// Software is furnished to do so, subject to the following conditions:
19//
20// The above copyright notice and this permission notice shall be included
21// in all copies or substantial portions of the Software.
22//
23// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
24// OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
25// FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
26// THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
27// LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
28// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
29// DEALINGS IN THE SOFTWARE.
30//
31// Description: Base Driver class for the stability solver
32//
33///////////////////////////////////////////////////////////////////////////////
34
35#include <iomanip>
36
38
39namespace Nektar::SolverUtils
40{
41
42/**
43 * Constructor
44 */
48 : Driver(pSession, pGraph)
49{
50 m_session->LoadParameter("IO_InfoSteps", m_infosteps, 1);
51}
52
53/**
54 * Arnoldi driver initialisation
55 */
56void DriverArnoldi::v_InitObject(std::ostream &out)
57{
59 m_session->MatchSolverInfo("SolverType", "VelocityCorrectionScheme",
61
63 {
64 m_period = m_session->GetParameter("TimeStep") *
65 m_session->GetParameter("NumSteps");
66 m_nfields = m_equ[0]->UpdateFields().size() - 1;
67 }
68 else
69 {
70 m_period = 1.0;
71 m_nfields = m_equ[0]->UpdateFields().size();
72 }
73
74 if (m_session->DefinesSolverInfo("ModeType") &&
75 (boost::iequals(m_session->GetSolverInfo("ModeType"), "SingleMode") ||
76 boost::iequals(m_session->GetSolverInfo("ModeType"), "HalfMode")))
77 {
78 for (int i = 0; i < m_nfields; ++i)
79 {
80 m_equ[0]->UpdateFields()[i]->SetWaveSpace(true);
81 }
82 }
83 m_negatedOp = m_equ[0]->NegatedOp();
84
85 m_session->LoadParameter("kdim", m_kdim, 16);
86 m_session->LoadParameter("nvec", m_nvec, 2);
87 m_session->LoadParameter("nits", m_nits, 500);
88 m_session->LoadParameter("evtol", m_evtol, 1e-06);
89
90 ASSERTL0(m_kdim >= m_nvec, "nvec cannot be larger than kdim.");
91
92 m_session->LoadParameter("realShift", m_realShift, 0.0);
93 m_equ[0]->SetLambda(m_realShift);
94
95 m_session->LoadParameter("imagShift", m_imagShift, 0.0);
96
97 // The imaginary shift is applied at system assembly
98 // Only if using HOMOGENEOUS expansion and ModeType set to SingleMode
99 if (m_imagShift != 0.0)
100 {
101 if (!m_session->DefinesSolverInfo("HOMOGENEOUS") &&
102 !m_session->DefinesSolverInfo("ModeType"))
103 {
105 "Imaginary shift only supported with HOMOGENEOUS "
106 "expansion and ModeType set to SingleMode");
107 }
108 else if (!boost::iequals(m_session->GetSolverInfo("ModeType"),
109 "SingleMode"))
110 {
112 "Imaginary shift only supported with HOMOGENEOUS "
113 "expansion and ModeType set to SingleMode");
114 }
115 }
116 MaskInit();
117}
118
119/**
120 *
121 */
122void DriverArnoldi::v_Execute([[maybe_unused]] std::ostream &out)
123{
124 ASSERTL0(false, "Specific version of Arnoldi driver not implemented");
125}
126
127/**
128 *
129 */
130void DriverArnoldi::ArnoldiSummary(std::ostream &out)
131{
132 if (m_comm->GetRank() == 0)
133 {
134 if (m_session->DefinesSolverInfo("ModeType") &&
135 boost::iequals(m_session->GetSolverInfo("ModeType"), "SingleMode"))
136 {
137 out << "\tSingle Fourier mode : true " << std::endl;
138 ASSERTL0(m_session->DefinesSolverInfo("Homogeneous"),
139 "Expected a homogeneous expansion to be defined "
140 "with single mode");
141 }
142 else
143 {
144 out << "\tSingle Fourier mode : false " << std::endl;
145 }
146 if (m_session->DefinesSolverInfo("BetaZero"))
147 {
148 out << "\tBeta set to Zero : true (overrides LHom)"
149 << std::endl;
150 }
151 else
152 {
153 out << "\tBeta set to Zero : false " << std::endl;
154 }
155
157 {
158 out << "\tEvolution operator : "
159 << m_session->GetSolverInfo("EvolutionOperator") << std::endl;
160 }
161 else
162 {
163 out << "\tShift (Real,Imag) : " << m_realShift << ","
164 << m_imagShift << std::endl;
165 }
166 out << "\tKrylov-space dimension : " << m_kdim << std::endl;
167 out << "\tNumber of vectors : " << m_nvec << std::endl;
168 out << "\tMax iterations : " << m_nits << std::endl;
169 out << "\tEigenvalue tolerance : " << m_evtol << std::endl;
170 out << "======================================================"
171 << std::endl;
172 }
173}
174
175/**
176 * Copy Arnoldi array to field variables which depend from
177 * either the m_fields or m_forces
178 */
180{
181
183 m_equ[0]->UpdateFields();
184 int nq = fields[0]->GetNcoeffs();
185
186 for (int k = 0; k < m_nfields; ++k)
187 {
188 Vmath::Vcopy(nq, &array[k * nq], 1, &fields[k]->UpdateCoeffs()[0], 1);
189 fields[k]->SetPhysState(false);
190 }
191}
192
193/**
194 * Copy field variables which depend from either the m_fields
195 * or m_forces array the Arnoldi array
196 */
198{
199
201
203 {
204 // This matters for the Adaptive SFD method because
205 // m_equ[1] is the nonlinear problem with non
206 // homogeneous BCs.
207 fields = m_equ[0]->UpdateFields();
208 }
209 else
210 {
211 fields = m_equ[m_nequ - 1]->UpdateFields();
212 }
213
214 for (int k = 0; k < m_nfields; ++k)
215 {
216 int nq = fields[0]->GetNcoeffs();
217 Vmath::Vcopy(nq, &fields[k]->GetCoeffs()[0], 1, &array[k * nq], 1);
218 fields[k]->SetPhysState(false);
219 }
220}
221
222/**
223 * Initialisation for the transient growth
224 */
226{
228
230 "Transient Growth non available for Coupled Solver");
231
232 fields = m_equ[0]->UpdateFields();
233 int nq = fields[0]->GetNcoeffs();
234
235 for (int k = 0; k < m_nfields; ++k)
236 {
237 Vmath::Vcopy(nq, &fields[k]->GetCoeffs()[0], 1,
238 &m_equ[1]->UpdateFields()[k]->UpdateCoeffs()[0], 1);
239 }
240}
241
242/**
243 *
244 */
245void DriverArnoldi::WriteFld(std::string file,
246 std::vector<Array<OneD, NekDouble>> coeffs)
247{
248
249 std::vector<std::string> variables(m_nfields);
250
251 ASSERTL1(coeffs.size() >= m_nfields, "coeffs is not of the correct length");
252 for (int i = 0; i < m_nfields; ++i)
253 {
254 variables[i] = m_equ[0]->GetVariable(i);
255 }
256
257 m_equ[0]->WriteFld(file, m_equ[0]->UpdateFields()[0], coeffs, variables);
258}
259
260/**
261 *
262 */
264{
265
266 std::vector<std::string> variables(m_nfields);
267 std::vector<Array<OneD, NekDouble>> fieldcoeffs(m_nfields);
268
269 int ncoeffs = m_equ[0]->UpdateFields()[0]->GetNcoeffs();
270 ASSERTL1(coeffs.size() >= ncoeffs * m_nfields,
271 "coeffs is not of sufficient size");
272
273 for (int i = 0; i < m_nfields; ++i)
274 {
275 variables[i] = m_equ[0]->GetVariable(i);
276 fieldcoeffs[i] = coeffs + i * ncoeffs;
277 }
278
279 m_equ[0]->WriteFld(file, m_equ[0]->UpdateFields()[0], fieldcoeffs,
280 variables);
281}
282
283/**
284 *
285 */
286void DriverArnoldi::WriteEvs(std::ostream &evlout, const int i,
287 const NekDouble re_ev, const NekDouble im_ev,
288 NekDouble resid, bool DumpInverse)
289{
291 {
292 NekDouble abs_ev = hypot(re_ev, im_ev);
293 NekDouble ang_ev = atan2(im_ev, re_ev);
294
295 evlout << "EV: " << std::setw(2) << i << std::setw(12) << abs_ev
296 << std::setw(12) << ang_ev << std::setw(12)
297 << log(abs_ev) / m_period << std::setw(12) << ang_ev / m_period;
298
300 {
301 evlout << std::setw(12) << resid;
302 }
303 evlout << std::endl;
304 }
305 else
306 {
307 NekDouble invmag = 1.0 / (re_ev * re_ev + im_ev * im_ev);
309 if (m_negatedOp)
310 {
311 sign = -1.0;
312 }
313 else
314 {
315 sign = 1.0;
316 }
317
318 evlout << "EV: " << std::setw(2) << i << std::setw(14) << sign * re_ev
319 << std::setw(14) << sign * im_ev;
320
321 if (DumpInverse)
322 {
323 evlout << std::setw(14) << sign * re_ev * invmag + m_realShift
324 << std::setw(14) << sign * im_ev * invmag + m_imagShift;
325 }
326
328 {
329 evlout << std::setw(12) << resid;
330 }
331 evlout << std::endl;
332 }
333}
334
335/**
336 *
337 */
339 std::vector<std::vector<LibUtilities::EquationSharedPtr>> &selectedDomains,
340 std::set<int> &unselectedVariables)
341{
342 selectedDomains.clear();
343 std::string domain("SelectEVCalcDomain0");
344 std::string condition("C0");
345 for (size_t i = 0; i < 10; ++i)
346 {
347 domain[domain.size() - 1] = '0' + i;
348 if (!m_session->DefinesFunction(domain))
349 {
350 break;
351 }
352 for (size_t j = 0; j < 10; ++j)
353 {
354 condition[condition.size() - 1] = '0' + j;
355 if (!m_session->DefinesFunction(domain, condition))
356 {
357 break;
358 }
359 if (j == 0)
360 {
361 selectedDomains.push_back(
362 std::vector<LibUtilities::EquationSharedPtr>());
363 }
364 selectedDomains[selectedDomains.size() - 1].push_back(
365 m_session->GetFunction(domain, condition));
366 }
367 }
368 unselectedVariables.clear();
369 std::string funName("SelectEVCalcVariables");
370 std::vector<std::string> variables = m_session->GetVariables();
371 for (size_t v = 0; v < m_nfields; ++v)
372 {
373 if (!m_session->DefinesFunction(funName, variables[v]))
374 {
375 unselectedVariables.insert(v);
376 }
377 }
378 if (unselectedVariables.size() == m_nfields)
379 {
380 unselectedVariables.clear();
381 }
382}
383
384/**
385 *
386 */
388{
389 std::vector<std::vector<LibUtilities::EquationSharedPtr>> selectedDomains;
390 std::set<int> unselectedVariables;
391 GetMaskInfo(selectedDomains, unselectedVariables);
392 if (selectedDomains.size() == 0 && unselectedVariables.size() == 0)
393 {
394 m_useMask = false;
395 return;
396 }
397 m_useMask = true;
398 MultiRegions::ExpListSharedPtr field = m_equ[0]->UpdateFields()[0];
399 int ncoef = field->GetNcoeffs();
400 int nphys = field->GetNpoints();
403 for (size_t i = 0; i < field->GetExpSize(); ++i)
404 {
406 SpatialDomains::GeometrySharedPtr geom = exp->GetGeom();
407 int nv = geom->GetNumVerts();
408 NekDouble gc[3] = {0., 0., 0.};
409 NekDouble gct[3] = {0., 0., 0.};
410 for (size_t j = 0; j < nv; ++j)
411 {
412 SpatialDomains::PointGeomSharedPtr vertex = geom->GetVertex(j);
413 vertex->GetCoords(gct[0], gct[1], gct[2]);
414 gc[0] += gct[0] / NekDouble(nv);
415 gc[1] += gct[1] / NekDouble(nv);
416 gc[2] += gct[2] / NekDouble(nv);
417 }
418 int unmask = 1;
419 for (size_t m = 0; m < selectedDomains.size(); ++m)
420 {
421 unmask = 0;
422 for (size_t n = 0; n < selectedDomains[m].size(); ++n)
423 {
424 if (selectedDomains[m][n]->Evaluate(gc[0], gc[1], gc[2]) <= 0.)
425 {
426 unmask = 0;
427 break;
428 }
429 else
430 {
431 unmask = 1;
432 }
433 }
434 if (unmask == 1)
435 {
436 break;
437 }
438 }
439 for (int j = 0; j < m_nfields; ++j)
440 {
441 if (unmask == 0 || unselectedVariables.count(j))
442 {
444 exp->GetNcoeffs(), 0.,
445 &m_maskCoeffs[field->GetCoeff_Offset(i) + j * ncoef], 1);
446 Vmath::Fill(exp->GetTotPoints(), 0.,
447 &m_maskPhys[field->GetPhys_Offset(i) + j * nphys],
448 1);
449 }
450 }
451 }
452}
453
454} // namespace Nektar::SolverUtils
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:208
#define NEKERROR(type, msg)
Assert Level 0 – Fundamental assert which is used whether in FULLDEBUG, DEBUG or OPT compilation mode...
Definition: ErrorUtil.hpp:202
#define ASSERTL1(condition, msg)
Assert Level 1 – Debugging which is used whether in FULLDEBUG or DEBUG compilation mode....
Definition: ErrorUtil.hpp:242
#define sign(a, b)
return the sign(b)*a
Definition: Polylib.cpp:47
void CopyFwdToAdj()
Copy the forward field to the adjoint system in transient growth calculations.
void v_InitObject(std::ostream &out=std::cout) override
Virtual function for initialisation implementation.
void WriteFld(std::string file, std::vector< Array< OneD, NekDouble > > coeffs)
Write coefficients to file.
NekDouble m_period
Tolerance of iterations.
Definition: DriverArnoldi.h:64
int m_infosteps
underlying operator is time stepping
Definition: DriverArnoldi.h:67
void v_Execute(std::ostream &out=std::cout) override
Virtual function for solve implementation.
Array< OneD, NekDouble > m_maskCoeffs
Definition: DriverArnoldi.h:78
DriverArnoldi(const LibUtilities::SessionReaderSharedPtr pSession, const SpatialDomains::MeshGraphSharedPtr pGraph)
Constructor.
void CopyFieldToArnoldiArray(Array< OneD, NekDouble > &array)
Copy fields to Arnoldi storage.
int m_nvec
Dimension of Krylov subspace.
Definition: DriverArnoldi.h:61
bool m_timeSteppingAlgorithm
Period of time stepping algorithm.
Definition: DriverArnoldi.h:65
int m_nits
Number of vectors to test.
Definition: DriverArnoldi.h:62
void CopyArnoldiArrayToField(Array< OneD, NekDouble > &array)
Copy Arnoldi storage to fields.
void GetMaskInfo(std::vector< std::vector< LibUtilities::EquationSharedPtr > > &selectedDomains, std::set< int > &unselectedVariables)
Array< OneD, NekDouble > m_maskPhys
Definition: DriverArnoldi.h:79
NekDouble m_evtol
Maxmum number of iterations.
Definition: DriverArnoldi.h:63
int m_nfields
interval to dump information if required.
Definition: DriverArnoldi.h:69
SOLVER_UTILS_EXPORT void ArnoldiSummary(std::ostream &out)
void WriteEvs(std::ostream &evlout, const int k, const NekDouble real, const NekDouble imag, NekDouble resid=NekConstants::kNekUnsetDouble, bool DumpInverse=true)
Base class for the development of solvers.
Definition: Driver.h:65
LibUtilities::SessionReaderSharedPtr m_session
Session reader object.
Definition: Driver.h:83
virtual SOLVER_UTILS_EXPORT void v_InitObject(std::ostream &out=std::cout)
Virtual function for initialisation implementation.
Definition: Driver.cpp:82
LibUtilities::CommSharedPtr m_comm
Communication object.
Definition: Driver.h:80
enum EvolutionOperatorType m_EvolutionOperator
Evolution Operator.
Definition: Driver.h:98
Array< OneD, EquationSystemSharedPtr > m_equ
Equation system to solve.
Definition: Driver.h:92
int m_nequ
number of equations
Definition: Driver.h:95
std::shared_ptr< SessionReader > SessionReaderSharedPtr
std::shared_ptr< Expansion > ExpansionSharedPtr
Definition: Expansion.h:66
std::shared_ptr< ExpList > ExpListSharedPtr
Shared pointer to an ExpList object.
static const NekDouble kNekUnsetDouble
std::shared_ptr< MeshGraph > MeshGraphSharedPtr
Definition: MeshGraph.h:174
std::shared_ptr< PointGeom > PointGeomSharedPtr
Definition: Geometry.h:57
std::shared_ptr< Geometry > GeometrySharedPtr
Definition: Geometry.h:51
double NekDouble
void Fill(int n, const T alpha, T *x, const int incx)
Fill a vector with a constant value.
Definition: Vmath.hpp:54
void Vcopy(int n, const T *x, const int incx, T *y, const int incy)
Definition: Vmath.hpp:825
scalarT< T > log(scalarT< T > in)
Definition: scalar.hpp:294