Nektar++
FilterEnergyBase.cpp
Go to the documentation of this file.
1 ///////////////////////////////////////////////////////////////////////////////
2 //
3 // File FilterEnergyBase.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 kinetic energy and enstrophy.
33 //
34 ///////////////////////////////////////////////////////////////////////////////
35 
36 #include <iomanip>
37 
39 
40 namespace Nektar
41 {
42 namespace SolverUtils
43 {
46  const ParamMap &pParams,
47  const bool pConstDensity)
48  : Filter (pSession),
49  m_index (-1),
50  m_homogeneous (false),
51  m_planes (),
52  m_constDensity(pConstDensity)
53 {
54  ParamMap::const_iterator it;
55  std::string outName;
56 
57  // OutputFile
58  it = pParams.find("OutputFile");
59  if (it == pParams.end())
60  {
61  outName = m_session->GetSessionName();
62  }
63  else
64  {
65  ASSERTL0(it->second.length() > 0, "Missing parameter 'OutputFile'.");
66  outName = it->second;
67  }
68  outName += ".eny";
69 
70  m_comm = pSession->GetComm();
71  if (m_comm->GetRank() == 0)
72  {
73  m_outFile.open(outName.c_str());
74  ASSERTL0(m_outFile.good(), "Unable to open: '" + outName + "'");
75  m_outFile.setf(ios::scientific, ios::floatfield);
76  m_outFile << "# Time Kinetic energy "
77  << "Enstrophy"
78  << endl
79  << "# ---------------------------------------------"
80  << "--------------"
81  << endl;
82  }
83  pSession->LoadParameter("LZ", m_homogeneousLength, 0.0);
84 
85  // OutputFrequency
86  it = pParams.find("OutputFrequency");
87  ASSERTL0(it != pParams.end(), "Missing parameter 'OutputFrequency'.");
88  LibUtilities::Equation equ(m_session, it->second);
89  m_outputFrequency = floor(equ.Evaluate());
90 }
91 
93 {
94 
95 }
96 
99  const NekDouble &time)
100 {
101  m_index = -1;
103 
104  ASSERTL0(pFields[0]->GetExpType() != MultiRegions::e3DH2D,
105  "Homogeneous 2D expansion not supported"
106  "for energy filter");
107 
108  if (pFields[0]->GetExpType() == MultiRegions::e3DH1D)
109  {
110  m_homogeneous = true;
111  }
112 
113  // Calculate area/volume of domain.
114  if (m_homogeneous)
115  {
116  m_planes = pFields[0]->GetZIDs();
117  areaField = pFields[0]->GetPlane(0);
118  }
119  else
120  {
121  areaField = pFields[0];
122  }
123 
124  Array<OneD, NekDouble> inarray(areaField->GetNpoints(), 1.0);
125  m_area = areaField->Integral(inarray);
126 
127  if (m_homogeneous)
128  {
130  }
131 
132  // Output values at initial time.
133  v_Update(pFields, time);
134 }
135 
138  const NekDouble &time)
139 {
140  int i, nPoints = pFields[0]->GetNpoints();
141 
142  m_index++;
143 
144  if (m_index % m_outputFrequency > 0)
145  {
146  return;
147  }
148 
149  // Calculate kinetic energy.
150  NekDouble Ek = 0.0;
151  Array<OneD, NekDouble> tmp(nPoints, 0.0);
153  for (i = 0; i < 3; ++i)
154  {
155  u[i] = Array<OneD, NekDouble>(nPoints);
156 
157  v_GetVelocity(pFields, i, u[i]);
158 
159  if (m_homogeneous && pFields[i]->GetWaveSpace())
160  {
161  pFields[i]->HomogeneousBwdTrans(u[i], u[i]);
162  }
163 
164  Vmath::Vvtvp(nPoints, u[i], 1, u[i], 1, tmp, 1, tmp, 1);
165  }
166 
167  if (!m_constDensity)
168  {
169  Vmath::Vmul(nPoints, v_GetDensity(pFields), 1, tmp, 1, tmp, 1);
170  }
171 
172  if (m_homogeneous)
173  {
174  Array<OneD, NekDouble> tmp2(nPoints, 0.0);
175  pFields[0]->HomogeneousFwdTrans(tmp, tmp2);
176  Ek = pFields[0]->GetPlane(0)->Integral(tmp2) * 2.0 * M_PI;
177  }
178  else
179  {
180  Ek = pFields[0]->Integral(tmp);
181  }
182 
183  Ek /= 2.0 * m_area;
184 
185  if (m_comm->GetRank() == 0)
186  {
187  m_outFile << setw(17) << setprecision(8) << time
188  << setw(22) << setprecision(11) << Ek;
189  }
190 
191  bool waveSpace[3] = {
192  pFields[0]->GetWaveSpace(),
193  pFields[1]->GetWaveSpace(),
194  pFields[2]->GetWaveSpace()
195  };
196 
197  if (m_homogeneous)
198  {
199  for (i = 0; i < 3; ++i)
200  {
201  pFields[i]->SetWaveSpace(false);
202  }
203  }
204 
205  // First calculate vorticity field.
206  Array<OneD, NekDouble> tmp2(nPoints), tmp3(nPoints);
207  Vmath::Zero(nPoints, tmp, 1);
208  for (i = 0; i < 3; ++i)
209  {
210  int f1 = (i+2) % 3, c2 = f1;
211  int c1 = (i+1) % 3, f2 = c1;
212  pFields[f1]->PhysDeriv(c1, u[f1], tmp2);
213  pFields[f2]->PhysDeriv(c2, u[f2], tmp3);
214  Vmath::Vsub (nPoints, tmp2, 1, tmp3, 1, tmp2, 1);
215  Vmath::Vvtvp(nPoints, tmp2, 1, tmp2, 1, tmp, 1, tmp, 1);
216  }
217 
218  if (!m_constDensity)
219  {
220  Vmath::Vmul(nPoints, v_GetDensity(pFields), 1, tmp, 1, tmp, 1);
221  }
222 
223  if (m_homogeneous)
224  {
225  for (i = 0; i < 3; ++i)
226  {
227  pFields[i]->SetWaveSpace(waveSpace[i]);
228  }
229  pFields[0]->HomogeneousFwdTrans(tmp, tmp);
230  Ek = pFields[0]->GetPlane(0)->Integral(tmp) * 2 * M_PI;
231  }
232  else
233  {
234  Ek = pFields[0]->Integral(tmp);
235  }
236 
237  Ek /= 2.0*m_area;
238 
239  if (m_comm->GetRank() == 0)
240  {
241  m_outFile << setw(22) << setprecision(11) << Ek << endl;
242  }
243 }
244 
247  const NekDouble &time)
248 {
249  m_outFile.close();
250 }
251 
253 {
254  return true;
255 }
256 
259  const int i,
260  Array<OneD, NekDouble> &velocity)
261 {
262  ASSERTL0(false, "Needs to implemented by subclass");
263 }
264 
267 {
268  ASSERTL0(false, "Needs to implemented by subclass");
269  return Array<OneD, NekDouble>();
270 }
271 }
272 }
virtual SOLVER_UTILS_EXPORT bool v_IsTimeDependent()
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:161
virtual SOLVER_UTILS_EXPORT Array< OneD, NekDouble > v_GetDensity(const Array< OneD, const MultiRegions::ExpListSharedPtr > &pFld)
void Vvtvp(int n, const T *w, const int incw, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
vvtvp (vector times vector plus vector): z = w*x + y
Definition: Vmath.cpp:428
virtual SOLVER_UTILS_EXPORT void v_Finalise(const Array< OneD, const MultiRegions::ExpListSharedPtr > &pField, const NekDouble &time)
Array< OneD, unsigned int > m_planes
LibUtilities::CommSharedPtr m_comm
boost::shared_ptr< SessionReader > SessionReaderSharedPtr
Definition: MeshPartition.h:50
virtual SOLVER_UTILS_EXPORT void v_Update(const Array< OneD, const MultiRegions::ExpListSharedPtr > &pField, const NekDouble &time)
boost::shared_ptr< ExpList > ExpListSharedPtr
Shared pointer to an ExpList object.
Definition: ExpList.h:1340
virtual SOLVER_UTILS_EXPORT void v_Initialise(const Array< OneD, const MultiRegions::ExpListSharedPtr > &pField, const NekDouble &time)
double NekDouble
std::map< std::string, std::string > ParamMap
Definition: Filter.h:67
LibUtilities::SessionReaderSharedPtr m_session
Definition: Filter.h:84
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
SOLVER_UTILS_EXPORT FilterEnergyBase(const LibUtilities::SessionReaderSharedPtr &pSession, const ParamMap &pParams, const bool pConstDensity=true)
void Zero(int n, T *x, const int incx)
Zero vector.
Definition: Vmath.cpp:359
void Vmul(int n, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
Multiply vector z = x*y.
Definition: Vmath.cpp:169
virtual SOLVER_UTILS_EXPORT void v_GetVelocity(const Array< OneD, const MultiRegions::ExpListSharedPtr > &pField, const int i, Array< OneD, NekDouble > &velocity)