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 // 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: Output values of the modal energy
33 //
34 ///////////////////////////////////////////////////////////////////////////////
35 
36 #include <iomanip>
37 
40 
41 namespace Nektar
42 {
43 namespace SolverUtils
44 {
46  RegisterCreatorFunction("ModalEnergy", FilterModalEnergy::create);
47 
48 /**
49  * Constructor.
50  */
53  const ParamMap &pParams) :
54  Filter(pSession)
55 {
56  ParamMap::const_iterator it;
57 
58  // OutputFile
59  it = pParams.find("OutputFile");
60  if (it == pParams.end())
61  {
62  m_outputFile = m_session->GetSessionName();
63  }
64  else
65  {
66  ASSERTL0(it->second.length() > 0, "Missing parameter 'OutputFile'.");
67  m_outputFile = it->second;
68  }
69  if (!(m_outputFile.length() >= 4
70  && m_outputFile.substr(m_outputFile.length() - 4) == ".mdl"))
71  {
72  m_outputFile += ".mdl";
73  }
74 
75  // OutputFrequency
76  it = pParams.find("OutputFrequency");
77  if (it == pParams.end())
78  {
80  }
81  else
82  {
83  LibUtilities::Equation equ(m_session, it->second);
84  m_outputFrequency = floor(equ.Evaluate());
85  }
86 
87 
88  m_session->MatchSolverInfo("Homogeneous", "1D", m_isHomogeneous1D, false);
89  m_session->MatchSolverInfo("Homogeneous", "2D", m_isHomogeneous2D, false);
90  m_session->MatchSolverInfo("CalculatePerturbationEnergy", "True",
91  m_PertEnergy, false);
92  m_session->LoadParameter ("NumQuadPointsError", m_NumQuadPointsError, 0);
93  m_EqTypeStr = m_session->GetSolverInfo("EQTYPE");
94 
95  // OutputPlane
97  {
98  m_session->LoadParameter("LZ", m_LhomZ);
99 
100  it = pParams.find("OutputPlane");
101  if (it == pParams.end())
102  {
103  m_outputPlane = 0;
104  }
105  else
106  {
107  LibUtilities::Equation equ(m_session, it->second);
108  m_outputPlane = floor(equ.Evaluate());
109  }
110  }
111 
113  AllocateSharedPtr(pSession->GetComm());
114 
115 }
116 
117 /**
118  * Destructor.
119  */
121 {
122 
123 }
124 
125 /**
126  * Initialize the parallel communication and the output stream.
127  */
130  const NekDouble &time)
131 {
132  LibUtilities::CommSharedPtr vComm = pFields[0]->GetComm();
133 
134  if (vComm->GetRank() == 0)
135  {
136  // Open output stream
138  {
139  m_outputStream.open(m_outputFile.c_str());
140  m_outputStream << "# Time, Fourier Mode, Energy ";
141  m_outputStream << endl;
142  }
143  else
144  {
145  m_outputStream.open(m_outputFile.c_str());
146  m_outputStream << "# Time, Energy ";
147  m_outputStream << endl;
148  }
149 
150  }
151 
152  m_index = 0;
153  v_Update(pFields, time);
154 }
155 
156 
157 /**
158  * Update the modal energy every m_outputFrequency.
159  */
162  const NekDouble &time)
163 {
164  // Only output every m_outputFrequency
165  if ((m_index++) % m_outputFrequency)
166  {
167  return;
168  }
169 
170  LibUtilities::CommSharedPtr vComm = pFields[0]->GetComm();
171 
172  // Homogeneous 1D implementation
173  if (m_isHomogeneous1D)
174  {
175  int colrank = vComm->GetColumnComm()->GetRank();
176  int nproc = vComm->GetColumnComm()->GetSize();
177  m_npointsZ = (m_session->GetParameter("HomModesZ"));
178  int locsize = m_npointsZ/nproc/2;
179 
180  Array<OneD, NekDouble> energy (locsize, 0.0);
181  Array<OneD, NekDouble> energy_tmp(locsize, 0.0);
183 
184  // Calculate the energy of the perturbation for stability
185  // analysis
186  if (m_PertEnergy)
187  {
188  // Compressible Flow Solver
189  if (m_EqTypeStr=="EulerCFE" ||
190  m_EqTypeStr=="EulerADCFE" ||
191  m_EqTypeStr=="NavierStokesCFE")
192  {
193  ASSERTL0(false, "Stability analysis module not "
194  "implemented for the Compressible Flow "
195  "Solver. Please remove the function BaseFlow "
196  "from your .xml file");
197  }
198  // Incompressible Navier-Stokes Solver
199  else
200  {
203  SetUpBaseFields(graphShrPtr);
204  string file = m_session->
205  GetFunctionFilename("BaseFlow", 0);
206  ImportFldBase(file, graphShrPtr);
207 
208  for (int i = 0; i < pFields.num_elements()-1; ++i)
209  {
210  Vmath::Vsub(pFields[i]->GetNcoeffs(),
211  pFields[i]->GetCoeffs(), 1,
212  m_base [i]->GetCoeffs(), 1,
213  pFields[i]->UpdateCoeffs(), 1);
214 
215  energy_tmp = pFields[i]->HomogeneousEnergy();
216  Vmath::Vadd(locsize, energy_tmp, 1,
217  energy, 1, energy, 1);
218 
219  Vmath::Vadd(pFields[i]->GetNcoeffs(),
220  pFields[i]->GetCoeffs(), 1,
221  m_base[i]->GetCoeffs(), 1,
222  pFields[i]->UpdateCoeffs(), 1);
223  }
224  }
225  }
226  // Calculate the modal energy for general simulation
227  else
228  {
229  // Compressible Flow Solver
230  if (m_EqTypeStr=="EulerCFE" ||
231  m_EqTypeStr=="EulerADCFE" ||
232  m_EqTypeStr=="NavierStokesCFE")
233  {
234  // Extracting kinetic energy
235  for (int i = 1; i < pFields.num_elements()-1; ++i)
236  {
237  energy_tmp = pFields[i]->HomogeneousEnergy();
238  Vmath::Vadd(locsize, energy_tmp, 1,
239  energy, 1, energy, 1);
240  }
241  }
242  // Incompressible Navier-Stokes Solver
243  else
244  {
245  // Extracting kinetic energy
246  for (int i = 0; i < pFields.num_elements()-1; ++i)
247  {
248  energy_tmp = pFields[i]->HomogeneousEnergy();
249  Vmath::Vadd(locsize, energy_tmp, 1,
250  energy, 1, energy, 1);
251  }
252  }
253  }
254 
255  // Send to root process
256  if (colrank == 0)
257  {
258  int j, m = 0;
259 
260  for (j = 0; j < energy.num_elements(); ++j, ++m)
261  {
262  m_outputStream << setw(10) << time
263  << setw(5) << m
264  << setw(18) << energy[j] << endl;
265  }
266 
267  for (int i = 1; i < nproc; ++i)
268  {
269  vComm->GetColumnComm()->Recv(i, energy);
270 
271  for (j = 0; j < energy.num_elements(); ++j, ++m)
272  {
273  m_outputStream << setw(10) << time
274  << setw(5) << m
275  << setw(18) << energy[j] << endl;
276  }
277  }
278  }
279  else
280  {
281  vComm->GetColumnComm()->Send(0, energy);
282  }
283  }
284  // Homogeneous 2D implementation
285  else if (m_isHomogeneous2D)
286  {
287  ASSERTL0(false, "3D Homogeneous 2D energy "
288  "dumping not implemented yet");
289  }
290  // General implementation
291  else
292  {
293  // Compressible Flow Solver
294  if (m_EqTypeStr=="EulerCFE" ||
295  m_EqTypeStr=="EulerADCFE" ||
296  m_EqTypeStr=="NavierStokesCFE")
297  {
298  // Total energy
299  NekDouble energy = 0.0;
300  for (int i = 1; i < pFields.num_elements()-1; ++i)
301  {
302  pFields[i]->SetPhysState(true);
303  NekDouble norm = L2Error(pFields, i, time);
304  energy += norm * norm;
305  }
306 
307  m_outputStream << setprecision(6) << time;
308  m_outputStream.width(25);
309  m_outputStream << setprecision(8) << 0.5*energy;
310  m_outputStream << endl;
311  }
312  // Incompressible Navier-Stokes Solver
313  else
314  {
315  // Kinetic energy
316  NekDouble energy = 0.0;
317  for (int i = 0; i < pFields.num_elements()-1; ++i)
318  {
319  pFields[i]->SetPhysState(true);
320  NekDouble norm = L2Error(pFields, i, time);
321  energy += norm * norm;
322  }
323  m_outputStream << setprecision(6) << time;
324  m_outputStream.width(25);
325  m_outputStream << setprecision(8) << 0.5*energy;
326  m_outputStream << endl;
327  }
328  }
329 }
330 
331 /**
332  * Close the output stream.
333  */
336  const NekDouble &time)
337 {
338  if (pFields[0]->GetComm()->GetRank() == 0)
339  {
340  m_outputStream.close();
341  }
342 }
343 
344 /**
345  * Calculate the L2 norm of a given field for calculating the
346  * modal energy.
347  */
350  unsigned int field,
351  const NekDouble &time)
352 {
353  NekDouble L2error = -1.0;
354  LibUtilities::CommSharedPtr vComm = pFields[0]->GetComm();
355 
356  if (m_NumQuadPointsError == 0)
357  {
358  if (pFields[field]->GetPhysState() == false)
359  {
360  pFields[field]->BwdTrans(pFields[field]->GetCoeffs(),
361  pFields[field]->UpdatePhys());
362  }
363  }
364 
365  L2error = pFields[field]->L2(pFields[field]->GetPhys());
366  return L2error;
367 }
368 
369 /**
370  * Setup the base fields in case of stability analyses.
371  */
374 {
375  int i;
376  int m_expdim = graphShrPtr->GetMeshDimension();
377 
378  //definition of the projection tipe:
379  if(m_session->DefinesSolverInfo("PROJECTION"))
380  {
381  std::string ProjectStr = m_session->GetSolverInfo("PROJECTION");
382 
383  if ((ProjectStr == "Continuous") ||
384  (ProjectStr == "Galerkin") ||
385  (ProjectStr == "CONTINUOUS") ||
386  (ProjectStr == "GALERKIN"))
387  {
389  }
390  else if ((ProjectStr == "MixedCGDG") ||
391  (ProjectStr == "Mixed_CG_Discontinuous"))
392  {
394  }
395  else if(ProjectStr == "DisContinuous")
396  {
398  }
399  else
400  {
401  ASSERTL0(false, "PROJECTION value not recognised");
402  }
403  }
404  else
405  {
406  cerr << "Projection type not specified in SOLVERINFO,"
407  "defaulting to continuous Galerkin" << endl;
409  }
410 
411  if (m_session->DefinesSolverInfo("ModeType"))
412  {
413  m_session->MatchSolverInfo("ModeType", "SingleMode",
414  m_SingleMode, false);
415  m_session->MatchSolverInfo("ModeType", "HalfMode",
416  m_HalfMode, false);
417  m_session->MatchSolverInfo("ModeType", "MultipleModes",
418  m_MultipleModes, false);
419  }
420 
421  m_session->MatchSolverInfo("USEFFT","FFTW", m_useFFT, false);
422  m_session->MatchSolverInfo("DEALIASING", "True",
423  m_homogen_dealiasing, false);
424 
425  if (m_homogen_dealiasing == false)
426  {
427  m_session->MatchSolverInfo("DEALIASING", "On",
428  m_homogen_dealiasing, false);
429  }
430 
431  // Stability Analysis flags
432  if (m_session->DefinesSolverInfo("ModeType"))
433  {
434  if (m_SingleMode)
435  {
436  m_npointsZ = 2;
437  }
438  else if (m_HalfMode)
439  {
440  m_npointsZ = 1;
441  }
442  else if (m_MultipleModes)
443  {
444  m_npointsZ = m_session->GetParameter("HomModesZ");
445  }
446  else
447  {
448  ASSERTL0(false, "SolverInfo ModeType not valid");
449  }
450  }
451  else
452  {
453  m_npointsZ = m_session->GetParameter("HomModesZ");
454  }
455 
458  {
459  switch (m_expdim)
460  {
461  case 1:
462  {
463  for(i = 0; i < m_base.num_elements(); i++)
464  {
466  ::AllocateSharedPtr(m_session, graphShrPtr,
467  m_session->GetVariable(0));
468  }
469  }
470  break;
471  case 2:
472  {
473  if (m_isHomogeneous1D)
474  {
475  if (m_SingleMode)
476  {
477  const LibUtilities::PointsKey PkeyZ(
478  m_npointsZ,
480  const LibUtilities::BasisKey BkeyZ(
482  m_npointsZ, PkeyZ);
483 
484  for (i = 0 ; i < m_base.num_elements(); i++)
485  {
488  AllocateSharedPtr(
489  m_session, BkeyZ, m_LhomZ,
491  graphShrPtr,
492  m_session->GetVariable(i));
493 
494  m_base[i]->SetWaveSpace(true);
495  }
496  }
497  else if (m_HalfMode)
498  {
499  //1 plane field (half mode expansion)
500  const LibUtilities::PointsKey PkeyZ(
501  m_npointsZ,
503  const LibUtilities::BasisKey BkeyZ(
505  m_npointsZ,PkeyZ);
506 
507  for (i = 0 ; i < m_base.num_elements(); i++)
508  {
511  AllocateSharedPtr(
512  m_session, BkeyZ, m_LhomZ,
514  graphShrPtr,
515  m_session->GetVariable(i));
516 
517  m_base[i]->SetWaveSpace(true);
518  }
519  }
520  else
521  {
522  const LibUtilities::PointsKey PkeyZ(
523  m_npointsZ,
525  const LibUtilities::BasisKey BkeyZ(
527 
528  for (i = 0 ; i < m_base.num_elements(); i++)
529  {
532  AllocateSharedPtr(
533  m_session, BkeyZ, m_LhomZ,
535  graphShrPtr,
536  m_session->GetVariable(i));
537 
538  m_base[i]->SetWaveSpace(false);
539  }
540  }
541  }
542  else
543  {
544  i = 0;
548  m_session,graphShrPtr,
549  m_session->GetVariable(i));
550 
551  m_base[0] = firstbase;
552 
553  for (i = 1 ; i < m_base.num_elements(); i++)
554  {
555  m_base[i] = MemoryManager<MultiRegions::
556  ContField2D>::AllocateSharedPtr(
557  *firstbase, graphShrPtr,
558  m_session->GetVariable(i));
559  }
560  }
561  }
562  break;
563  case 3:
564  {
567  AllocateSharedPtr(m_session, graphShrPtr,
568  m_session->GetVariable(0));
569  m_base[0] = firstbase;
570  for (i = 1 ; i < m_base.num_elements(); i++)
571  {
572  m_base[i] = MemoryManager<MultiRegions::
573  ContField3D>::AllocateSharedPtr(
574  *firstbase, graphShrPtr,
575  m_session->GetVariable(0));
576  }
577  }
578  break;
579  default:
580  ASSERTL0(false, "Expansion dimension not recognised");
581  break;
582  }
583  }
584  else
585  {
586  switch (m_expdim)
587  {
588  case 1:
589  {
590  // need to use zero for variable as may be more base
591  // flows than variables
592  for (i = 0 ; i < m_base.num_elements(); i++)
593  {
595  DisContField1D>::AllocateSharedPtr(
596  m_session, graphShrPtr,
597  m_session->GetVariable(0));
598  }
599  break;
600  }
601  case 2:
602  {
603  for (i = 0 ; i < m_base.num_elements(); i++)
604  {
606  DisContField2D>::AllocateSharedPtr(
607  m_session, graphShrPtr,
608  m_session->GetVariable(0));
609  }
610  break;
611  }
612  case 3:
613  ASSERTL0(false, "3D not set up");
614  default:
615  ASSERTL0(false, "Expansion dimension not recognised");
616  break;
617  }
618  }
619 }
620 
621 /**
622  * Import the base flow fld file.
623  */
625  std::string pInfile,
627 {
628  std::vector<LibUtilities::FieldDefinitionsSharedPtr> FieldDef;
629  std::vector<std::vector<NekDouble> > FieldData;
630 
631  // Get Homogeneous
632  m_fld->Import(pInfile,FieldDef,FieldData);
633 
634  int nvar = m_session->GetVariables().size();
635  if(m_session->DefinesSolverInfo("HOMOGENEOUS"))
636  {
637  std::string HomoStr = m_session->GetSolverInfo("HOMOGENEOUS");
638  }
639  // Copy FieldData into m_fields
640  for (int j = 0; j < nvar; ++j)
641  {
642  for (int i = 0; i < FieldDef.size(); ++i)
643  {
644  bool flag =
645  FieldDef[i]->m_fields[j] == m_session->GetVariable(j);
646 
647  ASSERTL0(flag, (std::string("Order of ") + pInfile
648  + std::string(" data and that defined in "
649  "m_boundaryconditions differs")).c_str());
650 
651  m_base[j]->ExtractDataToCoeffs(FieldDef[i], FieldData[i],
652  FieldDef[i]->m_fields[j],
653  m_base[j]->UpdateCoeffs());
654  }
655  }
656 }
657 
658 /**
659  * Flag for time-dependent flows.
660  */
662 {
663  return true;
664 }
665 }
666 }
static FilterSharedPtr create(const LibUtilities::SessionReaderSharedPtr &pSession, const ParamMap &pParams)
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:161
enum MultiRegions::ProjectionType m_projectionType
static boost::shared_ptr< MeshGraph > Read(const LibUtilities::SessionReaderSharedPtr &pSession, DomainRangeShPtr &rng=NullDomainRangeShPtr)
Definition: MeshGraph.cpp:119
static boost::shared_ptr< DataType > AllocateSharedPtr()
Allocate a shared pointer from the memory pool.
General purpose memory allocation routines with the ability to allocate from thread specific memory p...
virtual void v_Initialise(const Array< OneD, const MultiRegions::ExpListSharedPtr > &pFields, const NekDouble &time)
virtual void v_Finalise(const Array< OneD, const MultiRegions::ExpListSharedPtr > &pFields, const NekDouble &time)
boost::shared_ptr< ContField2D > ContField2DSharedPtr
Definition: ContField2D.h:291
Array< OneD, MultiRegions::ExpListSharedPtr > m_base
Fourier Expansion .
Definition: BasisType.h:52
boost::shared_ptr< SessionReader > SessionReaderSharedPtr
Definition: MeshPartition.h:50
LibUtilities::FieldIOSharedPtr m_fld
boost::shared_ptr< Comm > CommSharedPtr
Pointer to a Communicator object.
Definition: Comm.h:53
1D Evenly-spaced points using Fourier Fit
Definition: PointsType.h:64
Fourier Modified expansions with just the real part of the first mode .
Definition: BasisType.h:59
This class is the abstraction of a global continuous two- dimensional spectral/hp element expansion w...
Definition: ContField2D.h:56
void SetUpBaseFields(SpatialDomains::MeshGraphSharedPtr &mesh)
NekDouble Evaluate() const
Definition: Equation.h:102
This class is the abstraction of a global discontinuous two- dimensional spectral/hp element expansio...
Defines a specification for a set of points.
Definition: Points.h:58
double NekDouble
void ImportFldBase(std::string pInfile, SpatialDomains::MeshGraphSharedPtr pGraph)
std::map< std::string, std::string > ParamMap
Definition: Filter.h:67
LibUtilities::SessionReaderSharedPtr m_session
Definition: Filter.h:84
NekDouble L2Error(const Array< OneD, const MultiRegions::ExpListSharedPtr > &pFields, unsigned int field, const NekDouble &time)
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:329
virtual void v_Update(const Array< OneD, const MultiRegions::ExpListSharedPtr > &pFields, const NekDouble &time)
1D Non Evenly-spaced points for Single Mode analysis
Definition: PointsType.h:65
SOLVER_UTILS_EXPORT FilterModalEnergy(const LibUtilities::SessionReaderSharedPtr &pSession, const ParamMap &pParams)
FilterFactory & GetFilterFactory()
Definition: Filter.cpp:42
virtual SOLVER_UTILS_EXPORT ~FilterModalEnergy()
boost::shared_ptr< ContField3D > ContField3DSharedPtr
Definition: ContField3D.h:190
boost::shared_ptr< MeshGraph > MeshGraphSharedPtr
Definition: MeshGraph.h:432
Describes the specification for a Basis.
Definition: Basis.h:50
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.cpp:285