Nektar++
FilterModalEnergy.cpp
Go to the documentation of this file.
1///////////////////////////////////////////////////////////////////////////////
2//
3// File: FilterModalEnergy.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: Output values of the modal energy
32//
33///////////////////////////////////////////////////////////////////////////////
34
35#include <iomanip>
36
40
41using namespace std;
42
43namespace Nektar::SolverUtils
44{
48
49/**
50 * Constructor.
51 */
54 const std::shared_ptr<EquationSystem> &pEquation, const ParamMap &pParams)
55 : Filter(pSession, pEquation)
56{
57 // OutputFile
58 std::string ext = ".mdl";
59 m_outputFile = Filter::SetupOutput(ext, pParams);
60
61 // OutputFrequency
62 auto it = pParams.find("OutputFrequency");
63 if (it == pParams.end())
64 {
66 }
67 else
68 {
69 LibUtilities::Equation equ(m_session->GetInterpreter(), it->second);
70 m_outputFrequency = round(equ.Evaluate());
71 }
72
73 m_session->MatchSolverInfo("Homogeneous", "1D", m_isHomogeneous1D, false);
74 m_session->MatchSolverInfo("Homogeneous", "2D", m_isHomogeneous2D, false);
75 m_session->MatchSolverInfo("CalculatePerturbationEnergy", "True",
76 m_PertEnergy, false);
77 m_session->LoadParameter("NumQuadPointsError", m_NumQuadPointsError, 0);
78 m_EqTypeStr = m_session->GetSolverInfo("EQTYPE");
79
80 // OutputPlane
82 {
83 m_session->LoadParameter("LZ", m_LhomZ);
84
85 it = pParams.find("OutputPlane");
86 if (it == pParams.end())
87 {
88 m_outputPlane = 0;
89 }
90 else
91 {
92 LibUtilities::Equation equ(m_session->GetInterpreter(), it->second);
93 m_outputPlane = round(equ.Evaluate());
94 }
95 }
96
98}
99
100/**
101 * Destructor.
102 */
104{
105}
106
107/**
108 * Initialize the parallel communication and the output stream.
109 */
112 [[maybe_unused]] const NekDouble &time)
113{
114 LibUtilities::CommSharedPtr vComm = pFields[0]->GetComm();
115
116 if (vComm->GetRank() == 0)
117 {
118 // Open output stream
119 bool adaptive;
120 m_session->MatchSolverInfo("Driver", "Adaptive", adaptive, false);
121 if (adaptive)
122 {
123 m_outputStream.open(m_outputFile.c_str(), ofstream::app);
124 }
125 else
126 {
127 m_outputStream.open(m_outputFile.c_str());
128 }
130 {
131 m_outputStream << "# Time, Fourier Mode, Energy ";
132 m_outputStream << endl;
133 }
134 else
135 {
136 m_outputStream << "# Time, Energy ";
137 m_outputStream << endl;
138 }
139 }
140
141 m_index = 0;
143 {
144 v_Update(pFields, time);
145 }
146}
147
148/**
149 * Update the modal energy every m_outputFrequency.
150 */
153 const NekDouble &time)
154{
155 // Only output every m_outputFrequency
156 if ((m_index++) % m_outputFrequency)
157 {
158 return;
159 }
160
161 LibUtilities::CommSharedPtr vComm = pFields[0]->GetComm();
162
163 // Homogeneous 1D implementation
165 {
166 int colrank = vComm->GetColumnComm()->GetRank();
167 int nproc = vComm->GetColumnComm()->GetSize();
168 m_npointsZ = (m_session->GetParameter("HomModesZ"));
169 int locsize = m_npointsZ / nproc / 2;
170
171 Array<OneD, NekDouble> energy(locsize, 0.0);
172 Array<OneD, NekDouble> energy_tmp(locsize, 0.0);
174
175 // Calculate the energy of the perturbation for stability
176 // analysis
177 if (m_PertEnergy)
178 {
179 // Compressible Flow Solver
180 if (m_EqTypeStr == "EulerCFE" || m_EqTypeStr == "EulerADCFE" ||
181 m_EqTypeStr == "NavierStokesCFE")
182 {
183 ASSERTL0(false, "Stability analysis module not "
184 "implemented for the Compressible Flow "
185 "Solver. Please remove the function BaseFlow "
186 "from your .xml file");
187 }
188 // Incompressible Navier-Stokes Solver
189 else
190 {
193 SetUpBaseFields(graphShrPtr);
194 string file = m_session->GetFunctionFilename("BaseFlow", 0);
195 ImportFldBase(file);
196
197 for (int i = 0; i < pFields.size() - 1; ++i)
198 {
199 Vmath::Vsub(pFields[i]->GetNcoeffs(),
200 pFields[i]->GetCoeffs(), 1,
201 m_base[i]->GetCoeffs(), 1,
202 pFields[i]->UpdateCoeffs(), 1);
203
204 energy_tmp = pFields[i]->HomogeneousEnergy();
205 Vmath::Vadd(locsize, energy_tmp, 1, energy, 1, energy, 1);
206
207 Vmath::Vadd(pFields[i]->GetNcoeffs(),
208 pFields[i]->GetCoeffs(), 1,
209 m_base[i]->GetCoeffs(), 1,
210 pFields[i]->UpdateCoeffs(), 1);
211 }
212 }
213 }
214 // Calculate the modal energy for general simulation
215 else
216 {
217 // Compressible Flow Solver
218 if (m_EqTypeStr == "EulerCFE" || m_EqTypeStr == "EulerADCFE" ||
219 m_EqTypeStr == "NavierStokesCFE")
220 {
221 // Extracting kinetic energy
222 for (int i = 1; i < pFields.size() - 1; ++i)
223 {
224 energy_tmp = pFields[i]->HomogeneousEnergy();
225 Vmath::Vadd(locsize, energy_tmp, 1, energy, 1, energy, 1);
226 }
227 }
228 // Incompressible Navier-Stokes Solver
229 else
230 {
231 // Extracting kinetic energy
232 for (int i = 0; i < pFields.size() - 1; ++i)
233 {
234 energy_tmp = pFields[i]->HomogeneousEnergy();
235 Vmath::Vadd(locsize, energy_tmp, 1, energy, 1, energy, 1);
236 }
237 }
238 }
239
240 // Send to root process
241 if (colrank == 0)
242 {
243 int j, m = 0;
244
245 for (j = 0; j < energy.size(); ++j, ++m)
246 {
247 m_outputStream << setw(10) << time << setw(5) << m << setw(18)
248 << energy[j] << endl;
249 }
250
251 for (int i = 1; i < nproc; ++i)
252 {
253 vComm->GetColumnComm()->Recv(i, energy);
254
255 for (j = 0; j < energy.size(); ++j, ++m)
256 {
257 m_outputStream << setw(10) << time << setw(5) << m
258 << setw(18) << energy[j] << endl;
259 }
260 }
261 }
262 else
263 {
264 vComm->GetColumnComm()->Send(0, energy);
265 }
266 }
267 // Homogeneous 2D implementation
268 else if (m_isHomogeneous2D)
269 {
270 ASSERTL0(false, "3D Homogeneous 2D energy "
271 "dumping not implemented yet");
272 }
273 // General implementation
274 else
275 {
276 // Compressible Flow Solver
277 if (m_EqTypeStr == "EulerCFE" || m_EqTypeStr == "EulerADCFE" ||
278 m_EqTypeStr == "NavierStokesCFE")
279 {
280 // Total energy
281 NekDouble energy = 0.0;
282 for (int i = 1; i < pFields.size() - 1; ++i)
283 {
284 pFields[i]->SetPhysState(true);
285 NekDouble norm = L2Error(pFields, i, time);
286 energy += norm * norm;
287 }
288
289 m_outputStream << setprecision(6) << time;
290 m_outputStream.width(25);
291 m_outputStream << setprecision(8) << 0.5 * energy;
292 m_outputStream << endl;
293 }
294 // Incompressible Navier-Stokes Solver
295 else
296 {
297 // Kinetic energy
298 NekDouble energy = 0.0;
299 for (int i = 0; i < pFields.size() - 1; ++i)
300 {
301 pFields[i]->SetPhysState(true);
302 NekDouble norm = L2Error(pFields, i, time);
303 energy += norm * norm;
304 }
305 m_outputStream << setprecision(6) << time;
306 m_outputStream.width(25);
307 m_outputStream << setprecision(8) << 0.5 * energy;
308 m_outputStream << endl;
309 }
310 }
311}
312
313/**
314 * Close the output stream.
315 */
318 [[maybe_unused]] const NekDouble &time)
319{
320 if (pFields[0]->GetComm()->GetRank() == 0)
321 {
322 m_outputStream.close();
323 }
324}
325
326/**
327 * Calculate the L2 norm of a given field for calculating the
328 * modal energy.
329 */
332 unsigned int field, [[maybe_unused]] const NekDouble &time)
333{
334 NekDouble L2error = -1.0;
335 LibUtilities::CommSharedPtr vComm = pFields[0]->GetComm();
336
337 if (m_NumQuadPointsError == 0)
338 {
339 if (pFields[field]->GetPhysState() == false)
340 {
341 pFields[field]->BwdTrans(pFields[field]->GetCoeffs(),
342 pFields[field]->UpdatePhys());
343 }
344 }
345
346 L2error = pFields[field]->L2(pFields[field]->GetPhys());
347 return L2error;
348}
349
350/**
351 * Setup the base fields in case of stability analyses.
352 */
355{
356 int i;
357 int m_expdim = graphShrPtr->GetMeshDimension();
358
359 // definition of the projection tipe:
360 if (m_session->DefinesSolverInfo("PROJECTION"))
361 {
362 std::string ProjectStr = m_session->GetSolverInfo("PROJECTION");
363
364 if ((ProjectStr == "Continuous") || (ProjectStr == "Galerkin") ||
365 (ProjectStr == "CONTINUOUS") || (ProjectStr == "GALERKIN"))
366 {
368 }
369 else if ((ProjectStr == "MixedCGDG") ||
370 (ProjectStr == "Mixed_CG_Discontinuous"))
371 {
373 }
374 else if (ProjectStr == "DisContinuous")
375 {
377 }
378 else
379 {
380 ASSERTL0(false, "PROJECTION value not recognised");
381 }
382 }
383 else
384 {
385 cerr << "Projection type not specified in SOLVERINFO,"
386 "defaulting to continuous Galerkin"
387 << endl;
389 }
390
391 if (m_session->DefinesSolverInfo("ModeType"))
392 {
393 m_session->MatchSolverInfo("ModeType", "SingleMode", m_SingleMode,
394 false);
395 m_session->MatchSolverInfo("ModeType", "HalfMode", m_HalfMode, false);
396 m_session->MatchSolverInfo("ModeType", "MultipleModes", m_MultipleModes,
397 false);
398 }
399
400 m_session->MatchSolverInfo("USEFFT", "FFTW", m_useFFT, false);
401 m_session->MatchSolverInfo("DEALIASING", "True", m_homogen_dealiasing,
402 false);
403
404 // Stability Analysis flags
405 if (m_session->DefinesSolverInfo("ModeType"))
406 {
407 if (m_SingleMode)
408 {
409 m_npointsZ = 2;
410 }
411 else if (m_HalfMode)
412 {
413 m_npointsZ = 1;
414 }
415 else if (m_MultipleModes)
416 {
417 m_npointsZ = m_session->GetParameter("HomModesZ");
418 }
419 else
420 {
421 ASSERTL0(false, "SolverInfo ModeType not valid");
422 }
423 }
424 else
425 {
426 m_npointsZ = m_session->GetParameter("HomModesZ");
427 }
428
431 {
432 switch (m_expdim)
433 {
434 case 1:
435 {
436 for (i = 0; i < m_base.size(); i++)
437 {
439 AllocateSharedPtr(m_session, graphShrPtr,
440 m_session->GetVariable(0));
441 }
442 }
443 break;
444 case 2:
445 {
447 {
448 if (m_SingleMode)
449 {
450 const LibUtilities::PointsKey PkeyZ(
452 const LibUtilities::BasisKey BkeyZ(
454
455 for (i = 0; i < m_base.size(); i++)
456 {
459 AllocateSharedPtr(
460 m_session, BkeyZ, m_LhomZ, m_useFFT,
461 m_homogen_dealiasing, graphShrPtr,
462 m_session->GetVariable(i));
463
464 m_base[i]->SetWaveSpace(true);
465 }
466 }
467 else if (m_HalfMode)
468 {
469 // 1 plane field (half mode expansion)
470 const LibUtilities::PointsKey PkeyZ(
472 const LibUtilities::BasisKey BkeyZ(
474 PkeyZ);
475
476 for (i = 0; i < m_base.size(); i++)
477 {
480 AllocateSharedPtr(
481 m_session, BkeyZ, m_LhomZ, m_useFFT,
482 m_homogen_dealiasing, graphShrPtr,
483 m_session->GetVariable(i));
484
485 m_base[i]->SetWaveSpace(true);
486 }
487 }
488 else
489 {
490 const LibUtilities::PointsKey PkeyZ(
492 const LibUtilities::BasisKey BkeyZ(
494
495 for (i = 0; i < m_base.size(); i++)
496 {
499 AllocateSharedPtr(
500 m_session, BkeyZ, m_LhomZ, m_useFFT,
501 m_homogen_dealiasing, graphShrPtr,
502 m_session->GetVariable(i));
503
504 m_base[i]->SetWaveSpace(false);
505 }
506 }
507 }
508 else
509 {
510 i = 0;
513 AllocateSharedPtr(m_session, graphShrPtr,
514 m_session->GetVariable(i));
515
516 m_base[0] = firstbase;
517
518 for (i = 1; i < m_base.size(); i++)
519 {
521 AllocateSharedPtr(*firstbase, graphShrPtr,
522 m_session->GetVariable(i));
523 }
524 }
525 }
526 break;
527 case 3:
528 {
531 m_session, graphShrPtr, m_session->GetVariable(0));
532 m_base[0] = firstbase;
533 for (i = 1; i < m_base.size(); i++)
534 {
536 AllocateSharedPtr(*firstbase, graphShrPtr,
537 m_session->GetVariable(0));
538 }
539 }
540 break;
541 default:
543 "Expansion dimension not recognised");
544 break;
545 }
546 }
547 else
548 {
549 switch (m_expdim)
550 {
551 case 1:
552 {
553 // need to use zero for variable as may be more base
554 // flows than variables
555 for (i = 0; i < m_base.size(); i++)
556 {
558 AllocateSharedPtr(m_session, graphShrPtr,
559 m_session->GetVariable(0));
560 }
561 break;
562 }
563 case 2:
564 {
565 for (i = 0; i < m_base.size(); i++)
566 {
568 AllocateSharedPtr(m_session, graphShrPtr,
569 m_session->GetVariable(0));
570 }
571 break;
572 }
573 case 3:
574 NEKERROR(ErrorUtil::efatal, "3D not set up");
575 break;
576 default:
578 "Expansion dimension not recognised");
579 break;
580 }
581 }
582}
583
584/**
585 * Import the base flow fld file.
586 */
587void FilterModalEnergy::ImportFldBase(std::string pInfile)
588{
589 std::vector<LibUtilities::FieldDefinitionsSharedPtr> FieldDef;
590 std::vector<std::vector<NekDouble>> FieldData;
591
592 // Get Homogeneous
593 m_fld->Import(pInfile, FieldDef, FieldData);
594
595 int nvar = m_session->GetVariables().size();
596 if (m_session->DefinesSolverInfo("HOMOGENEOUS"))
597 {
598 std::string HomoStr = m_session->GetSolverInfo("HOMOGENEOUS");
599 }
600 // Copy FieldData into m_fields
601 for (int j = 0; j < nvar; ++j)
602 {
603 for (int i = 0; i < FieldDef.size(); ++i)
604 {
605 bool flag = FieldDef[i]->m_fields[j] == m_session->GetVariable(j);
606
607 ASSERTL0(flag, (std::string("Order of ") + pInfile +
608 std::string(" data and that defined in "
609 "m_boundaryconditions differs"))
610 .c_str());
611
612 m_base[j]->ExtractDataToCoeffs(FieldDef[i], FieldData[i],
613 FieldDef[i]->m_fields[j],
614 m_base[j]->UpdateCoeffs());
615 }
616 }
617}
618
619/**
620 * Flag for time-dependent flows.
621 */
623{
624 return true;
625}
626} // 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
Describes the specification for a Basis.
Definition: Basis.h:45
static std::shared_ptr< FieldIO > CreateDefault(const LibUtilities::SessionReaderSharedPtr session)
Returns an object for the default FieldIO method.
Definition: FieldIO.cpp:194
tKey RegisterCreatorFunction(tKey idKey, CreatorFunction classCreator, std::string pDesc="")
Register a class with the factory.
Defines a specification for a set of points.
Definition: Points.h:50
General purpose memory allocation routines with the ability to allocate from thread specific memory p...
static std::shared_ptr< DataType > AllocateSharedPtr(const Args &...args)
Allocate a shared pointer from the memory pool.
SOLVER_UTILS_EXPORT std::string SetupOutput(const std::string ext, const ParamMap &pParams)
Definition: Filter.h:139
LibUtilities::SessionReaderSharedPtr m_session
Definition: Filter.h:93
std::map< std::string, std::string > ParamMap
Definition: Filter.h:66
static FilterSharedPtr create(const LibUtilities::SessionReaderSharedPtr &pSession, const std::shared_ptr< EquationSystem > &pEquation, const ParamMap &pParams)
void v_Initialise(const Array< OneD, const MultiRegions::ExpListSharedPtr > &pFields, const NekDouble &time) override
void v_Finalise(const Array< OneD, const MultiRegions::ExpListSharedPtr > &pFields, const NekDouble &time) override
void v_Update(const Array< OneD, const MultiRegions::ExpListSharedPtr > &pFields, const NekDouble &time) override
enum MultiRegions::ProjectionType m_projectionType
SOLVER_UTILS_EXPORT FilterModalEnergy(const LibUtilities::SessionReaderSharedPtr &pSession, const std::shared_ptr< EquationSystem > &pEquation, const ParamMap &pParams)
SOLVER_UTILS_EXPORT ~FilterModalEnergy() override
LibUtilities::FieldIOSharedPtr m_fld
NekDouble L2Error(const Array< OneD, const MultiRegions::ExpListSharedPtr > &pFields, unsigned int field, const NekDouble &time)
Array< OneD, MultiRegions::ExpListSharedPtr > m_base
void SetUpBaseFields(SpatialDomains::MeshGraphSharedPtr &mesh)
static MeshGraphSharedPtr Read(const LibUtilities::SessionReaderSharedPtr pSession, LibUtilities::DomainRangeShPtr rng=LibUtilities::NullDomainRangeShPtr, bool fillGraph=true, SpatialDomains::MeshGraphSharedPtr partitionedGraph=nullptr)
Definition: MeshGraphIO.cpp:51
std::shared_ptr< SessionReader > SessionReaderSharedPtr
@ eFourierEvenlySpaced
1D Evenly-spaced points using Fourier Fit
Definition: PointsType.h:74
@ eFourierSingleModeSpaced
1D Non Evenly-spaced points for Single Mode analysis
Definition: PointsType.h:75
std::shared_ptr< Comm > CommSharedPtr
Pointer to a Communicator object.
Definition: Comm.h:55
@ eFourierHalfModeRe
Fourier Modified expansions with just the real part of the first mode .
Definition: BasisType.h:66
@ eFourier
Fourier Expansion .
Definition: BasisType.h:55
std::shared_ptr< ContField > ContFieldSharedPtr
Definition: ContField.h:268
FilterFactory & GetFilterFactory()
std::shared_ptr< MeshGraph > MeshGraphSharedPtr
Definition: MeshGraph.h:174
double NekDouble
void Vadd(int n, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
Add vector z = x+y.
Definition: Vmath.hpp:180
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.hpp:220
STL namespace.