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