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 std::map<std::string, std::string> &pParams) :
54  Filter(pSession)
55  {
56  if (pParams.find("OutputFile") == pParams.end())
57  {
58  m_outputFile = m_session->GetSessionName();
59  }
60  else
61  {
62  ASSERTL0(!(pParams.find("OutputFile")->second.empty()),
63  "Missing parameter 'OutputFile'.");
64  m_outputFile = pParams.find("OutputFile")->second;
65  }
66  if (!(m_outputFile.length() >= 4
67  && m_outputFile.substr(m_outputFile.length() - 4) == ".mdl"))
68  {
69  m_outputFile += ".mdl";
70  }
71 
72  if (pParams.find("OutputFrequency") == pParams.end())
73  {
75  }
76  else
77  {
79  atoi(pParams.find("OutputFrequency")->second.c_str());
80  }
81 
82 
83  m_session->MatchSolverInfo("Homogeneous", "1D",
84  m_isHomogeneous1D, false);
85  m_session->MatchSolverInfo("Homogeneous", "2D",
86  m_isHomogeneous2D, false);
87  m_session->MatchSolverInfo("CalculatePerturbationEnergy", "True",
88  m_PertEnergy, false);
89  m_session->LoadParameter ("NumQuadPointsError",
91 
92  m_EqTypeStr = m_session->GetSolverInfo("EQTYPE");
93 
95  {
96  m_session->LoadParameter("LZ", m_LhomZ);
97 
98  if (pParams.find("OutputPlane") == pParams.end())
99  {
100  m_outputPlane = 0;
101  }
102  else
103  {
104  m_outputPlane =
105  atoi(pParams.find("OutputPlane")->second.c_str());
106  }
107  }
108 
110  AllocateSharedPtr(pSession->GetComm());
111 
112  }
113 
114  /**
115  * Destructor.
116  */
118  {
119 
120  }
121 
122  /**
123  * Initialize the parallel communication and the output stream.
124  */
127  const NekDouble &time)
128  {
129  LibUtilities::CommSharedPtr vComm = pFields[0]->GetComm();
130 
131  if (vComm->GetRank() == 0)
132  {
133  // Open output stream
135  {
136  m_outputStream.open(m_outputFile.c_str());
137  m_outputStream << "# Time, Fourier Mode, Energy ";
138  m_outputStream << endl;
139  }
140  else
141  {
142  m_outputStream.open(m_outputFile.c_str());
143  m_outputStream << "# Time, Energy ";
144  m_outputStream << endl;
145  }
146 
147  }
148 
149  m_index = 0;
150  v_Update(pFields, time);
151  }
152 
153 
154  /**
155  * Update the modal energy every m_outputFrequency.
156  */
159  const NekDouble &time)
160  {
161  // Only output every m_outputFrequency
162  if ((m_index++) % m_outputFrequency)
163  {
164  return;
165  }
166 
167  LibUtilities::CommSharedPtr vComm = pFields[0]->GetComm();
168 
169  // Homogeneous 1D implementation
170  if (m_isHomogeneous1D)
171  {
172  int colrank = vComm->GetColumnComm()->GetRank();
173  int nproc = vComm->GetColumnComm()->GetSize();
174  m_npointsZ = (m_session->GetParameter("HomModesZ"));
175  int locsize = m_npointsZ/nproc/2;
176 
177  Array<OneD, NekDouble> energy (locsize, 0.0);
178  Array<OneD, NekDouble> energy_tmp(locsize, 0.0);
180 
181  // Calculate the energy of the perturbation for stability
182  // analysis
183  if (m_PertEnergy)
184  {
185  // Compressible Flow Solver
186  if (m_EqTypeStr=="EulerCFE" ||
187  m_EqTypeStr=="EulerADCFE" ||
188  m_EqTypeStr=="NavierStokesCFE")
189  {
190  ASSERTL0(false, "Stability analysis module not "
191  "implemented for the Compressible Flow "
192  "Solver. Please remove the function BaseFlow "
193  "from your .xml file");
194  }
195  // Incompressible Navier-Stokes Solver
196  else
197  {
200  SetUpBaseFields(graphShrPtr);
201  string file = m_session->
202  GetFunctionFilename("BaseFlow", 0);
203  ImportFldBase(file, graphShrPtr);
204 
205  for (int i = 0; i < pFields.num_elements()-1; ++i)
206  {
207  Vmath::Vsub(pFields[i]->GetNcoeffs(),
208  pFields[i]->GetCoeffs(), 1,
209  m_base [i]->GetCoeffs(), 1,
210  pFields[i]->UpdateCoeffs(), 1);
211 
212  energy_tmp = pFields[i]->HomogeneousEnergy();
213  Vmath::Vadd(locsize, energy_tmp, 1,
214  energy, 1, energy, 1);
215 
216  Vmath::Vadd(pFields[i]->GetNcoeffs(),
217  pFields[i]->GetCoeffs(), 1,
218  m_base[i]->GetCoeffs(), 1,
219  pFields[i]->UpdateCoeffs(), 1);
220  }
221  }
222  }
223  // Calculate the modal energy for general simulation
224  else
225  {
226  // Compressible Flow Solver
227  if (m_EqTypeStr=="EulerCFE" ||
228  m_EqTypeStr=="EulerADCFE" ||
229  m_EqTypeStr=="NavierStokesCFE")
230  {
231  // Extracting kinetic energy
232  for (int i = 1; i < pFields.num_elements()-1; ++i)
233  {
234  energy_tmp = pFields[i]->HomogeneousEnergy();
235  Vmath::Vadd(locsize, energy_tmp, 1,
236  energy, 1, energy, 1);
237  }
238  }
239  // Incompressible Navier-Stokes Solver
240  else
241  {
242  // Extracting kinetic energy
243  for (int i = 0; i < pFields.num_elements()-1; ++i)
244  {
245  energy_tmp = pFields[i]->HomogeneousEnergy();
246  Vmath::Vadd(locsize, energy_tmp, 1,
247  energy, 1, energy, 1);
248  }
249  }
250  }
251 
252  // Send to root process
253  if (colrank == 0)
254  {
255  int j, m = 0;
256 
257  for (j = 0; j < energy.num_elements(); ++j, ++m)
258  {
259  m_outputStream << setw(10) << time
260  << setw(5) << m
261  << setw(18) << energy[j] << endl;
262  }
263 
264  for (int i = 1; i < nproc; ++i)
265  {
266  vComm->GetColumnComm()->Recv(i, energy);
267 
268  for (j = 0; j < energy.num_elements(); ++j, ++m)
269  {
270  m_outputStream << setw(10) << time
271  << setw(5) << m
272  << setw(18) << energy[j] << endl;
273  }
274  }
275  }
276  else
277  {
278  vComm->GetColumnComm()->Send(0, energy);
279  }
280  }
281  // Homogeneous 2D implementation
282  else if (m_isHomogeneous2D)
283  {
284  ASSERTL0(false, "3D Homogeneous 2D energy "
285  "dumping not implemented yet");
286  }
287  // General implementation
288  else
289  {
290  // Compressible Flow Solver
291  if (m_EqTypeStr=="EulerCFE" ||
292  m_EqTypeStr=="EulerADCFE" ||
293  m_EqTypeStr=="NavierStokesCFE")
294  {
295  // Total energy
296  NekDouble energy = 0.0;
297  for (int i = 1; i < pFields.num_elements()-1; ++i)
298  {
299  pFields[i]->SetPhysState(true);
300  NekDouble norm = L2Error(pFields, i, time);
301  energy += norm * norm;
302  }
303 
304  m_outputStream << setprecision(6) << time;
305  m_outputStream.width(25);
306  m_outputStream << setprecision(8) << 0.5*energy;
307  m_outputStream << endl;
308  }
309  // Incompressible Navier-Stokes Solver
310  else
311  {
312  // Kinetic energy
313  NekDouble energy = 0.0;
314  for (int i = 0; i < pFields.num_elements()-1; ++i)
315  {
316  pFields[i]->SetPhysState(true);
317  NekDouble norm = L2Error(pFields, i, time);
318  energy += norm * norm;
319  }
320  m_outputStream << setprecision(6) << time;
321  m_outputStream.width(25);
322  m_outputStream << setprecision(8) << 0.5*energy;
323  m_outputStream << endl;
324  }
325  }
326  }
327 
328  /**
329  * Close the output stream.
330  */
333  const NekDouble &time)
334  {
335  if (pFields[0]->GetComm()->GetRank() == 0)
336  {
337  m_outputStream.close();
338  }
339  }
340 
341  /**
342  * Calculate the L2 norm of a given field for calculating the
343  * modal energy.
344  */
347  unsigned int field,
348  const NekDouble &time)
349  {
350  NekDouble L2error = -1.0;
351  LibUtilities::CommSharedPtr vComm = pFields[0]->GetComm();
352 
353  if (m_NumQuadPointsError == 0)
354  {
355  if (pFields[field]->GetPhysState() == false)
356  {
357  pFields[field]->BwdTrans(pFields[field]->GetCoeffs(),
358  pFields[field]->UpdatePhys());
359  }
360  }
361 
362  L2error = pFields[field]->L2(pFields[field]->GetPhys());
363  return L2error;
364  }
365 
366  /**
367  * Setup the base fields in case of stability analyses.
368  */
371  {
372  int i;
373  int m_expdim = graphShrPtr->GetMeshDimension();
374 
375  //definition of the projection tipe:
376  if(m_session->DefinesSolverInfo("PROJECTION"))
377  {
378  std::string ProjectStr = m_session->GetSolverInfo("PROJECTION");
379 
380  if ((ProjectStr == "Continuous") ||
381  (ProjectStr == "Galerkin") ||
382  (ProjectStr == "CONTINUOUS") ||
383  (ProjectStr == "GALERKIN"))
384  {
386  }
387  else if ((ProjectStr == "MixedCGDG") ||
388  (ProjectStr == "Mixed_CG_Discontinuous"))
389  {
391  }
392  else if(ProjectStr == "DisContinuous")
393  {
395  }
396  else
397  {
398  ASSERTL0(false, "PROJECTION value not recognised");
399  }
400  }
401  else
402  {
403  cerr << "Projection type not specified in SOLVERINFO,"
404  "defaulting to continuous Galerkin" << endl;
406  }
407 
408  if (m_session->DefinesSolverInfo("ModeType"))
409  {
410  m_session->MatchSolverInfo("ModeType", "SingleMode",
411  m_SingleMode, false);
412  m_session->MatchSolverInfo("ModeType", "HalfMode",
413  m_HalfMode, false);
414  m_session->MatchSolverInfo("ModeType", "MultipleModes",
415  m_MultipleModes, false);
416  }
417 
418  m_session->MatchSolverInfo("USEFFT","FFTW", m_useFFT, false);
419  m_session->MatchSolverInfo("DEALIASING", "True",
420  m_homogen_dealiasing, false);
421 
422  if (m_homogen_dealiasing == false)
423  {
424  m_session->MatchSolverInfo("DEALIASING", "On",
425  m_homogen_dealiasing, false);
426  }
427 
428  // Stability Analysis flags
429  if (m_session->DefinesSolverInfo("ModeType"))
430  {
431  if (m_SingleMode)
432  {
433  m_npointsZ = 2;
434  }
435  else if (m_HalfMode)
436  {
437  m_npointsZ = 1;
438  }
439  else if (m_MultipleModes)
440  {
441  m_npointsZ = m_session->GetParameter("HomModesZ");
442  }
443  else
444  {
445  ASSERTL0(false, "SolverInfo ModeType not valid");
446  }
447  }
448  else
449  {
450  m_npointsZ = m_session->GetParameter("HomModesZ");
451  }
452 
455  {
456  switch (m_expdim)
457  {
458  case 1:
459  {
460  for(i = 0; i < m_base.num_elements(); i++)
461  {
463  ::AllocateSharedPtr(m_session, graphShrPtr,
464  m_session->GetVariable(0));
465  }
466  }
467  break;
468  case 2:
469  {
470  if (m_isHomogeneous1D)
471  {
472  if (m_SingleMode)
473  {
474  const LibUtilities::PointsKey PkeyZ(
475  m_npointsZ,
477  const LibUtilities::BasisKey BkeyZ(
479  m_npointsZ, PkeyZ);
480 
481  for (i = 0 ; i < m_base.num_elements(); i++)
482  {
485  AllocateSharedPtr(
486  m_session, BkeyZ, m_LhomZ,
488  graphShrPtr,
489  m_session->GetVariable(i));
490 
491  m_base[i]->SetWaveSpace(true);
492  }
493  }
494  else if (m_HalfMode)
495  {
496  //1 plane field (half mode expansion)
497  const LibUtilities::PointsKey PkeyZ(
498  m_npointsZ,
500  const LibUtilities::BasisKey BkeyZ(
502  m_npointsZ,PkeyZ);
503 
504  for (i = 0 ; i < m_base.num_elements(); i++)
505  {
508  AllocateSharedPtr(
509  m_session, BkeyZ, m_LhomZ,
511  graphShrPtr,
512  m_session->GetVariable(i));
513 
514  m_base[i]->SetWaveSpace(true);
515  }
516  }
517  else
518  {
519  const LibUtilities::PointsKey PkeyZ(
520  m_npointsZ,
522  const LibUtilities::BasisKey BkeyZ(
524 
525  for (i = 0 ; i < m_base.num_elements(); i++)
526  {
529  AllocateSharedPtr(
530  m_session, BkeyZ, m_LhomZ,
532  graphShrPtr,
533  m_session->GetVariable(i));
534 
535  m_base[i]->SetWaveSpace(false);
536  }
537  }
538  }
539  else
540  {
541  i = 0;
545  m_session,graphShrPtr,
546  m_session->GetVariable(i));
547 
548  m_base[0] = firstbase;
549 
550  for (i = 1 ; i < m_base.num_elements(); i++)
551  {
552  m_base[i] = MemoryManager<MultiRegions::
553  ContField2D>::AllocateSharedPtr(
554  *firstbase, graphShrPtr,
555  m_session->GetVariable(i));
556  }
557  }
558  }
559  break;
560  case 3:
561  {
564  AllocateSharedPtr(m_session, graphShrPtr,
565  m_session->GetVariable(0));
566  m_base[0] = firstbase;
567  for (i = 1 ; i < m_base.num_elements(); i++)
568  {
569  m_base[i] = MemoryManager<MultiRegions::
570  ContField3D>::AllocateSharedPtr(
571  *firstbase, graphShrPtr,
572  m_session->GetVariable(0));
573  }
574  }
575  break;
576  default:
577  ASSERTL0(false, "Expansion dimension not recognised");
578  break;
579  }
580  }
581  else
582  {
583  switch (m_expdim)
584  {
585  case 1:
586  {
587  // need to use zero for variable as may be more base
588  // flows than variables
589  for (i = 0 ; i < m_base.num_elements(); i++)
590  {
592  DisContField1D>::AllocateSharedPtr(
593  m_session, graphShrPtr,
594  m_session->GetVariable(0));
595  }
596  break;
597  }
598  case 2:
599  {
600  for (i = 0 ; i < m_base.num_elements(); i++)
601  {
603  DisContField2D>::AllocateSharedPtr(
604  m_session, graphShrPtr,
605  m_session->GetVariable(0));
606  }
607  break;
608  }
609  case 3:
610  ASSERTL0(false, "3D not set up");
611  default:
612  ASSERTL0(false, "Expansion dimension not recognised");
613  break;
614  }
615  }
616  }
617 
618  /**
619  * Import the base flow fld file.
620  */
622  std::string pInfile,
624  {
625  std::vector<LibUtilities::FieldDefinitionsSharedPtr> FieldDef;
626  std::vector<std::vector<NekDouble> > FieldData;
627 
628  // Get Homogeneous
629  m_fld->Import(pInfile,FieldDef,FieldData);
630 
631  int nvar = m_session->GetVariables().size();
632  if(m_session->DefinesSolverInfo("HOMOGENEOUS"))
633  {
634  std::string HomoStr = m_session->GetSolverInfo("HOMOGENEOUS");
635  }
636  // Copy FieldData into m_fields
637  for (int j = 0; j < nvar; ++j)
638  {
639  for (int i = 0; i < FieldDef.size(); ++i)
640  {
641  bool flag =
642  FieldDef[i]->m_fields[j] == m_session->GetVariable(j);
643 
644  ASSERTL0(flag, (std::string("Order of ") + pInfile
645  + std::string(" data and that defined in "
646  "m_boundaryconditions differs")).c_str());
647 
648  m_base[j]->ExtractDataToCoeffs(FieldDef[i], FieldData[i],
649  FieldDef[i]->m_fields[j],
650  m_base[j]->UpdateCoeffs());
651  }
652  }
653  }
654 
655  /**
656  * Flag for time-dependent flows.
657  */
659  {
660  return true;
661  }
662  }
663 }
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:135
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
static FilterSharedPtr create(const LibUtilities::SessionReaderSharedPtr &pSession, const std::map< std::string, std::string > &pParams)
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
SOLVER_UTILS_EXPORT FilterModalEnergy(const LibUtilities::SessionReaderSharedPtr &pSession, const std::map< std::string, std::string > &pParams)
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)
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)
LibUtilities::SessionReaderSharedPtr m_session
Definition: Filter.h:76
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
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