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().RegisterCreatorFunction("ModalEnergy", FilterModalEnergy::create);
46 
47  /**
48  *
49  */
50  FilterModalEnergy::FilterModalEnergy(
51  const LibUtilities::SessionReaderSharedPtr &pSession,
52  const std::map<std::string, std::string> &pParams) :
53  Filter(pSession)
54  {
55  if (pParams.find("OutputFile") == pParams.end())
56  {
57  m_outputFile = m_session->GetSessionName();
58  }
59  else
60  {
61  ASSERTL0(!(pParams.find("OutputFile")->second.empty()),
62  "Missing parameter 'OutputFile'.");
63  m_outputFile = pParams.find("OutputFile")->second;
64  }
65  if (!(m_outputFile.length() >= 4
66  && m_outputFile.substr(m_outputFile.length() - 4) == ".mdl"))
67  {
68  m_outputFile += ".mdl";
69  }
70 
71  if (pParams.find("OutputFrequency") == pParams.end())
72  {
73  m_outputFrequency = 1;
74  }
75  else
76  {
77  m_outputFrequency =
78  atoi(pParams.find("OutputFrequency")->second.c_str());
79  }
80 
81 
82  m_session->MatchSolverInfo("Homogeneous", "1D",
83  m_isHomogeneous1D, false);
84  m_session->MatchSolverInfo("Homogeneous", "2D",
85  m_isHomogeneous2D, false);
86  m_session->MatchSolverInfo("CalculatePerturbationEnergy", "True",
87  m_PertEnergy, false);
88  m_session->LoadParameter ("NumQuadPointsError",
89  m_NumQuadPointsError, 0);
90 
91  if(m_isHomogeneous1D || m_isHomogeneous2D)
92  {
93  m_session->LoadParameter("LZ", m_LhomZ);
94 
95  if (pParams.find("OutputPlane") == pParams.end())
96  {
97  m_outputPlane = 0;
98  }
99  else
100  {
101  m_outputPlane =
102  atoi(pParams.find("OutputPlane")->second.c_str());
103  }
104  }
105 
106  m_fld = MemoryManager<LibUtilities::FieldIO>::AllocateSharedPtr(pSession->GetComm());
107 
108  }
109 
110  /**
111  *
112  */
113  FilterModalEnergy::~FilterModalEnergy()
114  {
115 
116  }
117 
118  /**
119  *
120  */
121  void FilterModalEnergy::v_Initialise(
122  const Array<OneD, const MultiRegions::ExpListSharedPtr> &pFields,
123  const NekDouble &time)
124  {
125  LibUtilities::CommSharedPtr vComm = pFields[0]->GetComm();
126 
127  if (vComm->GetRank() == 0)
128  {
129  // Open output stream
130  if(m_isHomogeneous1D)
131  {
132  m_outputStream.open(m_outputFile.c_str());
133  m_outputStream << "# Time, Fourier Mode, Energy ";
134  m_outputStream << endl;
135  }
136  else
137  {
138  m_outputStream.open(m_outputFile.c_str());
139  m_outputStream << "# Time, Energy ";
140  m_outputStream << endl;
141  }
142 
143  }
144 
145  v_Update(pFields, time);
146  }
147 
148 
149  /**
150  *
151  */
152  void FilterModalEnergy::v_Update(
153  const Array<OneD, const MultiRegions::ExpListSharedPtr> &pFields,
154  const NekDouble &time)
155  {
156  // Only output every m_outputFrequency.
157  if ((m_index++) % m_outputFrequency)
158  {
159  return;
160  }
161 
162  LibUtilities::CommSharedPtr vComm = pFields[0]->GetComm();
163 
164 
165  if(m_isHomogeneous1D)
166  {
167  int colrank = vComm->GetColumnComm()->GetRank();
168  int nproc = vComm->GetColumnComm()->GetSize();
169  m_npointsZ = (m_session->GetParameter("HomModesZ"));
170  int locsize = m_npointsZ/nproc/2;
171 
172  Array<OneD, NekDouble> energy (locsize,0.0);
173  Array<OneD, NekDouble> energy_tmp(locsize,0.0);
174  Array<OneD, NekDouble> tmp;
175 
176  // Calculate the energy of the perturbation for stability
177  // analysis
178  if(m_PertEnergy)
179  {
180  SpatialDomains::MeshGraphSharedPtr graphShrPtr = SpatialDomains::MeshGraph::Read(m_session);
181  SetUpBaseFields(graphShrPtr);
182  string file = m_session->GetFunctionFilename("BaseFlow", 0);
183  ImportFldBase(file,graphShrPtr);
184 
185  for(int i = 0; i < pFields.num_elements()-1; ++i)
186  {
187  Vmath::Vsub(pFields[i]->GetNcoeffs(),
188  pFields[i]->GetCoeffs(),1,
189  m_base [i]->GetCoeffs(),1,
190  pFields[i]->UpdateCoeffs(),1);
191 
192  energy_tmp = pFields[i]->HomogeneousEnergy();
193  Vmath::Vadd(locsize,energy_tmp,1,energy,1,energy,1);
194 
195  Vmath::Vadd(pFields[i]->GetNcoeffs(),
196  pFields[i]->GetCoeffs(),1,
197  m_base[i]->GetCoeffs(),1,
198  pFields[i]->UpdateCoeffs(),1);
199  }
200  }
201  else
202  {
203  for(int i = 0; i < pFields.num_elements()-1; ++i)
204  {
205  energy_tmp =pFields[i]->HomogeneousEnergy();
206  Vmath::Vadd(locsize,energy_tmp,1,energy,1,energy,1);
207  }
208  }
209 
210  // Send to root process.
211  if (colrank == 0)
212  {
213  int j, m = 0;
214 
215  for (j = 0; j < energy.num_elements(); ++j, ++m)
216  {
217  m_outputStream << setw(10) << time
218  << setw(5) << m
219  << setw(18) << energy[j] << endl;
220  }
221 
222  for (int i = 1; i < nproc; ++i)
223  {
224  vComm->GetColumnComm()->Recv(i, energy);
225 
226  for (j = 0; j < energy.num_elements(); ++j, ++m)
227  {
228  m_outputStream << setw(10) << time
229  << setw(5) << m
230  << setw(18) << energy[j] << endl;
231  }
232  }
233  }
234  else
235  {
236  vComm->GetColumnComm()->Send(0, energy);
237  }
238  }
239  else if(m_isHomogeneous2D)
240  {
241  ASSERTL0(false,"3D Homogeneous 2D energy dumping not implemented yet");
242  }
243  else
244  {
245  NekDouble energy = 0.0;
246  for(int i = 0; i < pFields.num_elements()-1; ++i)
247  {
248  pFields[i]->SetPhysState(true);
249  NekDouble norm = L2Error(pFields,i, time);
250  energy += norm*norm;
251  }
252  m_outputStream << setprecision(6) << time;
253  m_outputStream.width(25);
254  m_outputStream << setprecision(8) << 0.5*energy;
255  m_outputStream << endl;
256  }
257  }
258 
259  /**
260  *
261  */
262  void FilterModalEnergy::v_Finalise(
263  const Array<OneD, const MultiRegions::ExpListSharedPtr> &pFields,
264  const NekDouble &time)
265  {
266  if (pFields[0]->GetComm()->GetRank() == 0)
267  {
268  m_outputStream.close();
269  }
270  }
271 
272  NekDouble FilterModalEnergy::L2Error(
273  const Array<OneD, const MultiRegions::ExpListSharedPtr> &pFields,
274  unsigned int field,
275  const NekDouble &time)
276  {
277  NekDouble L2error = -1.0;
278  LibUtilities::CommSharedPtr vComm = pFields[0]->GetComm();
279 
280  if(m_NumQuadPointsError == 0)
281  {
282  if(pFields[field]->GetPhysState() == false)
283  {
284  pFields[field]->BwdTrans(pFields[field]->GetCoeffs(),
285  pFields[field]->UpdatePhys());
286  }
287  }
288 
289  L2error = pFields[field]->L2(pFields[field]->GetPhys());
290  return L2error;
291  }
292 
293  void FilterModalEnergy::SetUpBaseFields(
294  SpatialDomains::MeshGraphSharedPtr &graphShrPtr)
295  {
296  int i;
297  int m_expdim = graphShrPtr->GetMeshDimension();
298 
299  //definition of the projection tipe:
300  if(m_session->DefinesSolverInfo("PROJECTION"))
301  {
302  std::string ProjectStr = m_session->GetSolverInfo("PROJECTION");
303 
304  if((ProjectStr == "Continuous") || (ProjectStr == "Galerkin")
305  || (ProjectStr == "CONTINUOUS") || (ProjectStr == "GALERKIN"))
306  {
307  m_projectionType = MultiRegions::eGalerkin;
308  }
309  else if((ProjectStr == "MixedCGDG")||(ProjectStr == "Mixed_CG_Discontinuous"))
310  {
311  m_projectionType = MultiRegions::eMixed_CG_Discontinuous;
312  }
313  else if(ProjectStr == "DisContinuous")
314  {
315  m_projectionType = MultiRegions::eDiscontinuous;
316  }
317  else
318  {
319  ASSERTL0(false,"PROJECTION value not recognised");
320  }
321  }
322  else
323  {
324  cerr << "Projection type not specified in SOLVERINFO,"
325  "defaulting to continuous Galerkin" << endl;
326  m_projectionType = MultiRegions::eGalerkin;
327  }
328 
329  if(m_session->DefinesSolverInfo("ModeType"))
330  {
331  m_session->MatchSolverInfo("ModeType", "SingleMode",
332  m_SingleMode, false);
333  m_session->MatchSolverInfo("ModeType", "HalfMode",
334  m_HalfMode, false);
335  m_session->MatchSolverInfo("ModeType", "MultipleModes",
336  m_MultipleModes, false);
337  }
338 
339  m_session->MatchSolverInfo("USEFFT","FFTW",m_useFFT,false);
340  m_session->MatchSolverInfo("DEALIASING","True",m_homogen_dealiasing,false);
341  if(m_homogen_dealiasing == false)
342  {
343  m_session->MatchSolverInfo("DEALIASING","On",m_homogen_dealiasing,false);
344  }
345 
346  // Stability Analysis flags
347  if(m_session->DefinesSolverInfo("ModeType"))
348  {
349  if(m_SingleMode)
350  {
351  m_npointsZ = 2;
352  }
353  else if(m_HalfMode)
354  {
355  m_npointsZ = 1;
356  }
357  else if(m_MultipleModes)
358  {
359  m_npointsZ = m_session->GetParameter("HomModesZ");
360  }
361  else
362  {
363  ASSERTL0(false, "SolverInfo ModeType not valid");
364  }
365  }
366  else
367  {
368  m_npointsZ = m_session->GetParameter("HomModesZ");
369  }
370 
371  if (m_projectionType == MultiRegions::eGalerkin ||
372  m_projectionType == MultiRegions::eMixed_CG_Discontinuous)
373  {
374  switch (m_expdim)
375  {
376  case 1:
377  {
378  for(i = 0; i < m_base.num_elements(); i++)
379  {
380  m_base[i] = MemoryManager<MultiRegions::ContField1D>
381  ::AllocateSharedPtr(m_session, graphShrPtr,
382  m_session->GetVariable(0));
383  }
384  }
385  break;
386  case 2:
387  {
388  if(m_isHomogeneous1D)
389  {
390  if(m_SingleMode)
391  {
392  const LibUtilities::PointsKey PkeyZ(m_npointsZ,
393  LibUtilities::eFourierSingleModeSpaced);
394  const LibUtilities::BasisKey BkeyZ(
395  LibUtilities::eFourier,
396  m_npointsZ,PkeyZ);
397 
398  for(i = 0 ; i < m_base.num_elements(); i++)
399  {
400  m_base[i] = MemoryManager<MultiRegions::ContField3DHomogeneous1D> ::AllocateSharedPtr(m_session,BkeyZ,m_LhomZ,m_useFFT,m_homogen_dealiasing,graphShrPtr,m_session->GetVariable(i));
401  m_base[i]->SetWaveSpace(true);
402  }
403  }
404  else if(m_HalfMode)
405  {
406  //1 plane field (half mode expansion)
407  const LibUtilities::PointsKey PkeyZ(m_npointsZ,
408  LibUtilities::eFourierSingleModeSpaced);
409  const LibUtilities::BasisKey BkeyZ(
410  LibUtilities::eFourierHalfModeRe,
411  m_npointsZ,PkeyZ);
412 
413  for(i = 0 ; i < m_base.num_elements(); i++)
414  {
415  m_base[i] = MemoryManager<MultiRegions::ContField3DHomogeneous1D>
416  ::AllocateSharedPtr(m_session,BkeyZ,m_LhomZ,m_useFFT,m_homogen_dealiasing,graphShrPtr,m_session->GetVariable(i));
417  m_base[i]->SetWaveSpace(true);
418  }
419  }
420  else
421  {
422  const LibUtilities::PointsKey PkeyZ(m_npointsZ,
423  LibUtilities::eFourierEvenlySpaced);
424  const LibUtilities::BasisKey BkeyZ(
425  LibUtilities::eFourier,m_npointsZ,
426  PkeyZ);
427 
428  for(i = 0 ; i < m_base.num_elements(); i++)
429  {
430  m_base[i] = MemoryManager<MultiRegions::ContField3DHomogeneous1D>
431  ::AllocateSharedPtr(m_session,BkeyZ,m_LhomZ,m_useFFT,m_homogen_dealiasing,graphShrPtr,m_session->GetVariable(i));
432  m_base[i]->SetWaveSpace(false);
433  }
434  }
435  }
436  else
437  {
438  i = 0;
439  MultiRegions::ContField2DSharedPtr firstbase =
440  MemoryManager<MultiRegions::ContField2D>
441  ::AllocateSharedPtr(m_session,graphShrPtr,
442  m_session->GetVariable(i));
443  m_base[0]=firstbase;
444 
445  for(i = 1 ; i < m_base.num_elements(); i++)
446  {
447  m_base[i] = MemoryManager<MultiRegions::ContField2D>
448  ::AllocateSharedPtr(*firstbase,graphShrPtr,
449  m_session->GetVariable(i));
450  }
451  }
452  }
453  break;
454  case 3:
455  {
456  MultiRegions::ContField3DSharedPtr firstbase =
457  MemoryManager<MultiRegions::ContField3D>
458  ::AllocateSharedPtr(m_session, graphShrPtr,
459  m_session->GetVariable(0));
460  m_base[0] = firstbase;
461  for(i = 1 ; i < m_base.num_elements(); i++)
462  {
463  m_base[i] = MemoryManager<MultiRegions::ContField3D>
464  ::AllocateSharedPtr(*firstbase, graphShrPtr,
465  m_session->GetVariable(0));
466  }
467  }
468  break;
469  default:
470  ASSERTL0(false,"Expansion dimension not recognised");
471  break;
472  }
473  }
474  else
475  {
476  switch(m_expdim)
477  {
478  case 1:
479  {
480  // need to use zero for variable as may be more base
481  // flows than variables
482  for(i = 0 ; i < m_base.num_elements(); i++)
483  {
484  m_base[i] = MemoryManager<MultiRegions
485  ::DisContField1D>
486  ::AllocateSharedPtr(m_session, graphShrPtr,
487  m_session->GetVariable(0));
488  }
489  break;
490  }
491  case 2:
492  {
493  for(i = 0 ; i < m_base.num_elements(); i++)
494  {
495  m_base[i] = MemoryManager<MultiRegions
496  ::DisContField2D>::AllocateSharedPtr(m_session,graphShrPtr,
497  m_session->GetVariable(0));
498  }
499  break;
500  }
501  case 3:
502  ASSERTL0(false, "3 D not set up");
503  default:
504  ASSERTL0(false, "Expansion dimension not recognised");
505  break;
506  }
507  }
508  }
509 
510  void FilterModalEnergy::ImportFldBase(
511  std::string pInfile,
512  SpatialDomains::MeshGraphSharedPtr pGraph)
513  {
514  std::vector<LibUtilities::FieldDefinitionsSharedPtr> FieldDef;
515  std::vector<std::vector<NekDouble> > FieldData;
516 
517  //Get Homogeneous
518  m_fld->Import(pInfile,FieldDef,FieldData);
519 
520  int nvar = m_session->GetVariables().size();
521  if(m_session->DefinesSolverInfo("HOMOGENEOUS"))
522  {
523  std::string HomoStr = m_session->GetSolverInfo("HOMOGENEOUS");
524  }
525  // copy FieldData into m_fields
526  for(int j = 0; j < nvar; ++j)
527  {
528  for(int i = 0; i < FieldDef.size(); ++i)
529  {
530  bool flag = FieldDef[i]->m_fields[j] ==
531  m_session->GetVariable(j);
532  ASSERTL1(flag, (std::string("Order of ") + pInfile
533  + std::string(" data and that defined in "
534  "m_boundaryconditions differs")).c_str());
535 
536  m_base[j]->ExtractDataToCoeffs(FieldDef[i], FieldData[i],
537  FieldDef[i]->m_fields[j],
538  m_base[j]->UpdateCoeffs());
539  }
540  }
541  }
542 
543  /**
544  *
545  */
546  bool FilterModalEnergy::v_IsTimeDependent()
547  {
548  return true;
549  }
550  }
551 }
Generated on Sun Mar 15 2015 00:11:31 for Nektar++ by   doxygen 1.8.1.2