Nektar++
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
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 using namespace std;
41 
42 namespace Nektar
43 {
44 namespace SolverUtils
45 {
46 FilterEnergyBase::FilterEnergyBase(
48  const ParamMap &pParams,
49  const bool pConstDensity)
50  : Filter (pSession),
51  m_index (-1),
52  m_homogeneous (false),
53  m_planes (),
54  m_constDensity(pConstDensity)
55 {
56  ParamMap::const_iterator it;
57  std::string outName;
58 
59  // OutputFile
60  it = pParams.find("OutputFile");
61  if (it == pParams.end())
62  {
63  outName = m_session->GetSessionName();
64  }
65  else
66  {
67  ASSERTL0(it->second.length() > 0, "Missing parameter 'OutputFile'.");
68  outName = it->second;
69  }
70  outName += ".eny";
71 
72  m_comm = pSession->GetComm();
73  if (m_comm->GetRank() == 0)
74  {
75  m_outFile.open(outName.c_str());
76  ASSERTL0(m_outFile.good(), "Unable to open: '" + outName + "'");
77  m_outFile.setf(ios::scientific, ios::floatfield);
78  m_outFile << "# Time Kinetic energy "
79  << "Enstrophy"
80  << endl
81  << "# ---------------------------------------------"
82  << "--------------"
83  << endl;
84  }
85  pSession->LoadParameter("LZ", m_homogeneousLength, 0.0);
86 
87  // OutputFrequency
88  it = pParams.find("OutputFrequency");
89  ASSERTL0(it != pParams.end(), "Missing parameter 'OutputFrequency'.");
90  LibUtilities::Equation equ(m_session, it->second);
91  m_outputFrequency = floor(equ.Evaluate());
92 }
93 
95 {
96 
97 }
98 
101  const NekDouble &time)
102 {
103  m_index = -1;
105 
106  ASSERTL0(pFields[0]->GetExpType() != MultiRegions::e3DH2D,
107  "Homogeneous 2D expansion not supported"
108  "for energy filter");
109 
110  if (pFields[0]->GetExpType() == MultiRegions::e3DH1D)
111  {
112  m_homogeneous = true;
113  }
114 
115  // Calculate area/volume of domain.
116  if (m_homogeneous)
117  {
118  m_planes = pFields[0]->GetZIDs();
119  areaField = pFields[0]->GetPlane(0);
120  }
121  else
122  {
123  areaField = pFields[0];
124  }
125 
126  Array<OneD, NekDouble> inarray(areaField->GetNpoints(), 1.0);
127  m_area = areaField->Integral(inarray);
128 
129  if (m_homogeneous)
130  {
132  }
133 
134  // Output values at initial time.
135  v_Update(pFields, time);
136 }
137 
140  const NekDouble &time)
141 {
142  int i, nPoints = pFields[0]->GetNpoints();
143 
144  m_index++;
145 
146  if (m_index % m_outputFrequency > 0)
147  {
148  return;
149  }
150 
151  // Calculate kinetic energy.
152  NekDouble Ek = 0.0;
153  Array<OneD, NekDouble> tmp(nPoints, 0.0);
155  for (i = 0; i < 3; ++i)
156  {
157  u[i] = Array<OneD, NekDouble>(nPoints);
158 
159  v_GetVelocity(pFields, i, u[i]);
160 
161  if (m_homogeneous && pFields[i]->GetWaveSpace())
162  {
163  pFields[i]->HomogeneousBwdTrans(u[i], u[i]);
164  }
165 
166  Vmath::Vvtvp(nPoints, u[i], 1, u[i], 1, tmp, 1, tmp, 1);
167  }
168 
169  if (!m_constDensity)
170  {
171  Vmath::Vmul(nPoints, v_GetDensity(pFields), 1, tmp, 1, tmp, 1);
172  }
173 
174  if (m_homogeneous)
175  {
176  Array<OneD, NekDouble> tmp2(nPoints, 0.0);
177  pFields[0]->HomogeneousFwdTrans(tmp, tmp2);
178  Ek = pFields[0]->GetPlane(0)->Integral(tmp2) * 2.0 * M_PI;
179  }
180  else
181  {
182  Ek = pFields[0]->Integral(tmp);
183  }
184 
185  Ek /= 2.0 * m_area;
186 
187  if (m_comm->GetRank() == 0)
188  {
189  m_outFile << setw(17) << setprecision(8) << time
190  << setw(22) << setprecision(11) << Ek;
191  }
192 
193  bool waveSpace[3] = {
194  pFields[0]->GetWaveSpace(),
195  pFields[1]->GetWaveSpace(),
196  pFields[2]->GetWaveSpace()
197  };
198 
199  if (m_homogeneous)
200  {
201  for (i = 0; i < 3; ++i)
202  {
203  pFields[i]->SetWaveSpace(false);
204  }
205  }
206 
207  // First calculate vorticity field.
208  Array<OneD, NekDouble> tmp2(nPoints), tmp3(nPoints);
209  Vmath::Zero(nPoints, tmp, 1);
210  for (i = 0; i < 3; ++i)
211  {
212  int f1 = (i+2) % 3, c2 = f1;
213  int c1 = (i+1) % 3, f2 = c1;
214  pFields[f1]->PhysDeriv(c1, u[f1], tmp2);
215  pFields[f2]->PhysDeriv(c2, u[f2], tmp3);
216  Vmath::Vsub (nPoints, tmp2, 1, tmp3, 1, tmp2, 1);
217  Vmath::Vvtvp(nPoints, tmp2, 1, tmp2, 1, tmp, 1, tmp, 1);
218  }
219 
220  if (!m_constDensity)
221  {
222  Vmath::Vmul(nPoints, v_GetDensity(pFields), 1, tmp, 1, tmp, 1);
223  }
224 
225  if (m_homogeneous)
226  {
227  for (i = 0; i < 3; ++i)
228  {
229  pFields[i]->SetWaveSpace(waveSpace[i]);
230  }
231  pFields[0]->HomogeneousFwdTrans(tmp, tmp);
232  Ek = pFields[0]->GetPlane(0)->Integral(tmp) * 2 * M_PI;
233  }
234  else
235  {
236  Ek = pFields[0]->Integral(tmp);
237  }
238 
239  Ek /= 2.0*m_area;
240 
241  if (m_comm->GetRank() == 0)
242  {
243  m_outFile << setw(22) << setprecision(11) << Ek << endl;
244  }
245 }
246 
249  const NekDouble &time)
250 {
251  m_outFile.close();
252 }
253 
255 {
256  return true;
257 }
258 
261  const int i,
262  Array<OneD, NekDouble> &velocity)
263 {
264  ASSERTL0(false, "Needs to implemented by subclass");
265 }
266 
269 {
270  ASSERTL0(false, "Needs to implemented by subclass");
271  return Array<OneD, NekDouble>();
272 }
273 }
274 }
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
STL namespace.
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:51
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.
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
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)