Nektar++
FilterEnergy.cpp
Go to the documentation of this file.
1///////////////////////////////////////////////////////////////////////////////
2//
3// File: FilterEnergy.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// Permission is hereby granted, free of charge, to any person obtaining a
14// copy of this software and associated documentation files (the "Software"),
15// to deal in the Software without restriction, including without limitation
16// the rights to use, copy, modify, merge, publish, distribute, sublicense,
17// and/or sell copies of the Software, and to permit persons to whom the
18// Software is furnished to do so, subject to the following conditions:
19//
20// The above copyright notice and this permission notice shall be included
21// in all copies or substantial portions of the Software.
22//
23// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
24// OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
25// FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
26// THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
27// LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
28// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
29// DEALINGS IN THE SOFTWARE.
30//
31// Description: Output kinetic energy and enstrophy.
32//
33///////////////////////////////////////////////////////////////////////////////
34
35#include <iomanip>
36
39
40using namespace std;
41
42namespace Nektar::SolverUtils
43{
44std::string FilterEnergy::className =
46 "Energy", FilterEnergy::create);
47
49 const std::shared_ptr<EquationSystem> &pEquation,
50 const ParamMap &pParams)
51 : Filter(pSession, pEquation), m_index(-1), m_homogeneous(false), m_planes()
52{
53 std::string outName;
54 std::string ext = ".eny";
55 outName = Filter::SetupOutput(ext, pParams);
56
57 m_comm = pSession->GetComm();
58 if (m_comm->GetRank() == 0)
59 {
60 m_outFile.open(outName.c_str());
61 ASSERTL0(m_outFile.good(), "Unable to open: '" + outName + "'");
62 m_outFile.setf(ios::scientific, ios::floatfield);
63 m_outFile << "# Time Kinetic energy "
64 << "Enstrophy" << endl
65 << "# ---------------------------------------------"
66 << "--------------" << endl;
67 }
68 pSession->LoadParameter("LZ", m_homogeneousLength, 0.0);
69
70 // OutputFrequency
71 auto it = pParams.find("OutputFrequency");
72 ASSERTL0(it != pParams.end(), "Missing parameter 'OutputFrequency'.");
73 LibUtilities::Equation equ(m_session->GetInterpreter(), it->second);
74 m_outputFrequency = round(equ.Evaluate());
75}
76
78{
79}
80
83 [[maybe_unused]] const NekDouble &time)
84{
85 m_index = -1;
87
88 ASSERTL0(pFields[0]->GetExpType() != MultiRegions::e1D,
89 "1D expansion not supported for energy filter");
90
91 ASSERTL0(pFields[0]->GetExpType() != MultiRegions::e2D,
92 "2D expansion not supported for energy filter");
93
94 ASSERTL0(pFields[0]->GetExpType() != MultiRegions::e3DH2D,
95 "Homogeneous 2D expansion not supported for energy filter");
96
97 if (pFields[0]->GetExpType() == MultiRegions::e3DH1D)
98 {
99 m_homogeneous = true;
100 }
101
102 // Calculate area/volume of domain.
103 if (m_homogeneous)
104 {
105 m_planes = pFields[0]->GetZIDs();
106 areaField = pFields[0]->GetPlane(0);
107 }
108 else
109 {
110 areaField = pFields[0];
111 }
112
113 Array<OneD, NekDouble> inarray(areaField->GetNpoints(), 1.0);
114 m_area = areaField->Integral(inarray);
115
116 if (m_homogeneous)
117 {
119 }
121 {
122 v_Update(pFields, time);
123 }
124}
125
128 const NekDouble &time)
129{
130 int i, nPoints = pFields[0]->GetNpoints();
131
132 m_index++;
133
134 if (m_index % m_outputFrequency > 0)
135 {
136 return;
137 }
138
139 // Lock equation system pointer
140 auto equ = m_equ.lock();
141 ASSERTL0(equ, "Weak pointer expired");
142
143 auto fluidEqu = std::dynamic_pointer_cast<FluidInterface>(equ);
144 ASSERTL0(fluidEqu, "Energy filter is incompatible with this solver.");
145
146 // Store physical values in an array
147 Array<OneD, Array<OneD, NekDouble>> physfields(pFields.size());
148 for (i = 0; i < pFields.size(); ++i)
149 {
150 physfields[i] = pFields[i]->GetPhys();
151 }
152
153 // Calculate kinetic energy.
154 NekDouble Ek = 0.0;
155 Array<OneD, NekDouble> tmp(nPoints, 0.0);
158 for (i = 0; i < 3; ++i)
159 {
160 u[i] = Array<OneD, NekDouble>(nPoints);
161 }
162 fluidEqu->GetVelocity(physfields, u);
163
164 for (i = 0; i < 3; ++i)
165 {
166 if (m_homogeneous && pFields[i]->GetWaveSpace())
167 {
168 pFields[i]->HomogeneousBwdTrans(nPoints, u[i], u[i]);
169 }
170
171 Vmath::Vvtvp(nPoints, u[i], 1, u[i], 1, tmp, 1, tmp, 1);
172 }
173
174 if (!fluidEqu->HasConstantDensity())
175 {
176 density = Array<OneD, NekDouble>(nPoints);
177 fluidEqu->GetDensity(physfields, density);
178 Vmath::Vmul(nPoints, density, 1, tmp, 1, tmp, 1);
179 }
180
181 if (m_homogeneous)
182 {
183 Array<OneD, NekDouble> tmp2(nPoints, 0.0);
184 pFields[0]->HomogeneousFwdTrans(nPoints, tmp, tmp2);
185 Ek = pFields[0]->GetPlane(0)->Integral(tmp2) * m_homogeneousLength;
186 }
187 else
188 {
189 Ek = pFields[0]->Integral(tmp);
190 }
191
192 Ek /= 2.0 * m_area;
193
194 if (m_comm->GetRank() == 0)
195 {
196 m_outFile << setw(17) << setprecision(8) << time << setw(22)
197 << setprecision(11) << Ek;
198 }
199
200 bool waveSpace[3] = {pFields[0]->GetWaveSpace(), pFields[1]->GetWaveSpace(),
201 pFields[2]->GetWaveSpace()};
202
203 if (m_homogeneous)
204 {
205 for (i = 0; i < 3; ++i)
206 {
207 pFields[i]->SetWaveSpace(false);
208 }
209 }
210
211 // First calculate vorticity field.
212 Array<OneD, NekDouble> tmp2(nPoints), tmp3(nPoints);
213 Vmath::Zero(nPoints, tmp, 1);
214 for (i = 0; i < 3; ++i)
215 {
216 int f1 = (i + 2) % 3, c2 = f1;
217 int c1 = (i + 1) % 3, f2 = c1;
218 pFields[f1]->PhysDeriv(c1, u[f1], tmp2);
219 pFields[f2]->PhysDeriv(c2, u[f2], tmp3);
220 Vmath::Vsub(nPoints, tmp2, 1, tmp3, 1, tmp2, 1);
221 Vmath::Vvtvp(nPoints, tmp2, 1, tmp2, 1, tmp, 1, tmp, 1);
222 }
223
224 if (!fluidEqu->HasConstantDensity())
225 {
226 Vmath::Vmul(nPoints, density, 1, tmp, 1, tmp, 1);
227 }
228
229 if (m_homogeneous)
230 {
231 for (i = 0; i < 3; ++i)
232 {
233 pFields[i]->SetWaveSpace(waveSpace[i]);
234 }
235 pFields[0]->HomogeneousFwdTrans(nPoints, tmp, tmp);
236 Ek = pFields[0]->GetPlane(0)->Integral(tmp) * m_homogeneousLength;
237 }
238 else
239 {
240 Ek = pFields[0]->Integral(tmp);
241 }
242
243 Ek /= 2.0 * m_area;
244
245 if (m_comm->GetRank() == 0)
246 {
247 m_outFile << setw(22) << setprecision(11) << Ek << endl;
248 }
249}
250
253 &pFields,
254 [[maybe_unused]] const NekDouble &time)
255{
256 m_outFile.close();
257}
258
260{
261 return true;
262}
263
264} // namespace Nektar::SolverUtils
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:208
tKey RegisterCreatorFunction(tKey idKey, CreatorFunction classCreator, std::string pDesc="")
Register a class with the factory.
SOLVER_UTILS_EXPORT void v_Initialise(const Array< OneD, const MultiRegions::ExpListSharedPtr > &pField, const NekDouble &time) override
SOLVER_UTILS_EXPORT void v_Update(const Array< OneD, const MultiRegions::ExpListSharedPtr > &pField, const NekDouble &time) override
static std::string className
Name of the class.
Definition: FilterEnergy.h:58
SOLVER_UTILS_EXPORT void v_Finalise(const Array< OneD, const MultiRegions::ExpListSharedPtr > &pField, const NekDouble &time) override
LibUtilities::CommSharedPtr m_comm
Definition: FilterEnergy.h:85
SOLVER_UTILS_EXPORT FilterEnergy(const LibUtilities::SessionReaderSharedPtr &pSession, const std::shared_ptr< EquationSystem > &pEquation, const ParamMap &pParams)
Array< OneD, unsigned int > m_planes
Definition: FilterEnergy.h:86
SOLVER_UTILS_EXPORT ~FilterEnergy() override
static SolverUtils::FilterSharedPtr create(const LibUtilities::SessionReaderSharedPtr &pSession, const std::shared_ptr< SolverUtils::EquationSystem > &pEquation, const ParamMap &pParams)
Creates an instance of this class.
Definition: FilterEnergy.h:46
SOLVER_UTILS_EXPORT bool v_IsTimeDependent() override
SOLVER_UTILS_EXPORT std::string SetupOutput(const std::string ext, const ParamMap &pParams)
Definition: Filter.h:139
LibUtilities::SessionReaderSharedPtr m_session
Definition: Filter.h:93
const std::weak_ptr< EquationSystem > m_equ
Definition: Filter.h:94
std::map< std::string, std::string > ParamMap
Definition: Filter.h:66
std::shared_ptr< SessionReader > SessionReaderSharedPtr
std::shared_ptr< ExpList > ExpListSharedPtr
Shared pointer to an ExpList object.
FilterFactory & GetFilterFactory()
double NekDouble
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.hpp:72
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.hpp:366
void Zero(int n, T *x, const int incx)
Zero vector.
Definition: Vmath.hpp:273
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.hpp:220
STL namespace.