Nektar++
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
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 
38 #include <LibUtilities/Memory/NekMemoryManager.hpp>
39 #include <SolverUtils/Filters/FilterModalEnergy.h>
40 
41 namespace Nektar
42 {
43  namespace SolverUtils
44  {
45  std::string FilterModalEnergy::className = GetFilterFactory().
46  RegisterCreatorFunction("ModalEnergy", FilterModalEnergy::create);
47 
48  /**
49  * Constructor.
50  */
51  FilterModalEnergy::FilterModalEnergy(
52  const LibUtilities::SessionReaderSharedPtr &pSession,
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  {
74  m_outputFrequency = 1;
75  }
76  else
77  {
78  m_outputFrequency =
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",
90  m_NumQuadPointsError, 0);
91 
92  m_EqTypeStr = m_session->GetSolverInfo("EQTYPE");
93 
94  if(m_isHomogeneous1D || m_isHomogeneous2D)
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 
109  m_fld = MemoryManager<LibUtilities::FieldIO>::
110  AllocateSharedPtr(pSession->GetComm());
111 
112  }
113 
114  /**
115  * Destructor.
116  */
117  FilterModalEnergy::~FilterModalEnergy()
118  {
119 
120  }
121 
122  /**
123  * Initialize the parallel communication and the output stream.
124  */
125  void FilterModalEnergy::v_Initialise(
126  const Array<OneD, const MultiRegions::ExpListSharedPtr> &pFields,
127  const NekDouble &time)
128  {
129  LibUtilities::CommSharedPtr vComm = pFields[0]->GetComm();
130 
131  if (vComm->GetRank() == 0)
132  {
133  // Open output stream
134  if(m_isHomogeneous1D)
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  */
157  void FilterModalEnergy::v_Update(
158  const Array<OneD, const MultiRegions::ExpListSharedPtr> &pFields,
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);
179  Array<OneD, NekDouble> tmp;
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  {
198  SpatialDomains::MeshGraphSharedPtr graphShrPtr =
199  SpatialDomains::MeshGraph::Read(m_session);
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  */
331  void FilterModalEnergy::v_Finalise(
332  const Array<OneD, const MultiRegions::ExpListSharedPtr> &pFields,
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  */
345  NekDouble FilterModalEnergy::L2Error(
346  const Array<OneD, const MultiRegions::ExpListSharedPtr> &pFields,
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  */
369  void FilterModalEnergy::SetUpBaseFields(
370  SpatialDomains::MeshGraphSharedPtr &graphShrPtr)
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  {
385  m_projectionType = MultiRegions::eGalerkin;
386  }
387  else if ((ProjectStr == "MixedCGDG") ||
388  (ProjectStr == "Mixed_CG_Discontinuous"))
389  {
390  m_projectionType = MultiRegions::eMixed_CG_Discontinuous;
391  }
392  else if(ProjectStr == "DisContinuous")
393  {
394  m_projectionType = MultiRegions::eDiscontinuous;
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;
405  m_projectionType = MultiRegions::eGalerkin;
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 
453  if (m_projectionType == MultiRegions::eGalerkin ||
454  m_projectionType == MultiRegions::eMixed_CG_Discontinuous)
455  {
456  switch (m_expdim)
457  {
458  case 1:
459  {
460  for(i = 0; i < m_base.num_elements(); i++)
461  {
462  m_base[i] = MemoryManager<MultiRegions::ContField1D>
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,
476  LibUtilities::eFourierSingleModeSpaced);
477  const LibUtilities::BasisKey BkeyZ(
478  LibUtilities::eFourier,
479  m_npointsZ, PkeyZ);
480 
481  for (i = 0 ; i < m_base.num_elements(); i++)
482  {
483  m_base[i] = MemoryManager<MultiRegions::
484  ContField3DHomogeneous1D>::
485  AllocateSharedPtr(
486  m_session, BkeyZ, m_LhomZ,
487  m_useFFT, m_homogen_dealiasing,
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,
499  LibUtilities::eFourierSingleModeSpaced);
500  const LibUtilities::BasisKey BkeyZ(
501  LibUtilities::eFourierHalfModeRe,
502  m_npointsZ,PkeyZ);
503 
504  for (i = 0 ; i < m_base.num_elements(); i++)
505  {
506  m_base[i] = MemoryManager<MultiRegions::
507  ContField3DHomogeneous1D>::
508  AllocateSharedPtr(
509  m_session, BkeyZ, m_LhomZ,
510  m_useFFT, m_homogen_dealiasing,
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,
521  LibUtilities::eFourierEvenlySpaced);
522  const LibUtilities::BasisKey BkeyZ(
523  LibUtilities::eFourier, m_npointsZ, PkeyZ);
524 
525  for (i = 0 ; i < m_base.num_elements(); i++)
526  {
527  m_base[i] = MemoryManager<MultiRegions::
528  ContField3DHomogeneous1D>::
529  AllocateSharedPtr(
530  m_session, BkeyZ, m_LhomZ,
531  m_useFFT, m_homogen_dealiasing,
532  graphShrPtr,
533  m_session->GetVariable(i));
534 
535  m_base[i]->SetWaveSpace(false);
536  }
537  }
538  }
539  else
540  {
541  i = 0;
542  MultiRegions::ContField2DSharedPtr firstbase =
543  MemoryManager<MultiRegions::ContField2D>::
544  AllocateSharedPtr(
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  {
562  MultiRegions::ContField3DSharedPtr firstbase =
563  MemoryManager<MultiRegions::ContField3D>::
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  {
591  m_base[i] = MemoryManager<MultiRegions::
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  {
602  m_base[i] = MemoryManager<MultiRegions::
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  */
621  void FilterModalEnergy::ImportFldBase(
622  std::string pInfile,
623  SpatialDomains::MeshGraphSharedPtr pGraph)
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  */
658  bool FilterModalEnergy::v_IsTimeDependent()
659  {
660  return true;
661  }
662  }
663 }
Generated on Wed Mar 4 2015 16:33:14 for Nektar++ by   doxygen 1.8.1.2