Nektar++
NavierStokesAdvection.cpp
Go to the documentation of this file.
1///////////////////////////////////////////////////////////////////////////////
2//
3// File: NavierStokesAdvection.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: Evaluation of the Navier Stokes advective term
32//
33///////////////////////////////////////////////////////////////////////////////
34
37
38namespace Nektar
39{
40
43 "True", 0),
45 "False", 1)};
46
49 "Convective", NavierStokesAdvection::create, "Convective");
52 "NonConservative", NavierStokesAdvection::create, "NonConserviative");
53
54/**
55 * Constructor. Creates ...
56 *
57 * \param
58 * \param
59 */
60
62
63{
64}
65
69{
70 m_homogen_dealiasing = pSession->DefinesSolverInfo("dealiasing");
71
72 pSession->MatchSolverInfo("SPECTRALHPDEALIASING", "True",
73 m_specHP_dealiasing, false);
74 pSession->MatchSolverInfo("ModeType", "SingleMode", m_SingleMode, false);
75 pSession->MatchSolverInfo("ModeType", "HalfMode", m_HalfMode, false);
76
77 Advection::v_InitObject(pSession, pFields);
78}
79
81 const int nConvectiveFields,
83 const Array<OneD, Array<OneD, NekDouble>> &advVel,
84 const Array<OneD, Array<OneD, NekDouble>> &inarray,
86 [[maybe_unused]] const NekDouble &time,
87 [[maybe_unused]] const Array<OneD, Array<OneD, NekDouble>> &pFwd,
88 [[maybe_unused]] const Array<OneD, Array<OneD, NekDouble>> &pBwd)
89{
90 int nqtot = fields[0]->GetTotPoints();
91 ASSERTL1(nConvectiveFields == inarray.size(),
92 "Number of convective fields and Inarray are not compatible");
93
94 // use dimension of Velocity vector to dictate dimension of operation
95 int ndim = advVel.size();
96 Array<OneD, Array<OneD, NekDouble>> AdvVel(advVel.size());
97
99
101
102 for (int i = 0; i < ndim; ++i)
103 {
104 if (fields[i]->GetWaveSpace() && !m_SingleMode && !m_HalfMode &&
106 {
107 velocity[i] = Array<OneD, NekDouble>(nqtot, 0.0);
108 fields[i]->HomogeneousBwdTrans(nqtot, advVel[i], velocity[i]);
109 }
110 else
111 {
112 velocity[i] = advVel[i];
113 }
114 }
115
116 int nPointsTot = fields[0]->GetNpoints();
117 Array<OneD, NekDouble> grad0, grad1, grad2, wkSp;
118
119 NekDouble OneDptscale = 1.5; // factor to rescale 1d points in dealiasing
120
122 {
123 // Get number of points to dealias a quadratic non-linearity
124 nPointsTot = fields[0]->Get1DScaledTotPoints(OneDptscale);
125 }
126
127 // interpolate Advection velocity
128 if (m_specHP_dealiasing) // interpolate advection field to higher space.
129 {
130 for (int i = 0; i < ndim; ++i)
131 {
132 AdvVel[i] = Array<OneD, NekDouble>(nPointsTot);
133 // interpolate infield to 3/2 dimension
134 timer.Start();
135 fields[0]->PhysInterp1DScaled(OneDptscale, velocity[i], AdvVel[i]);
136 timer.Stop();
137 timer.AccumulateRegion("Interp1DScaled");
138 }
139 }
140 else
141 {
142 for (int i = 0; i < ndim; ++i)
143 {
144 AdvVel[i] = velocity[i];
145 }
146 }
147
148 wkSp = Array<OneD, NekDouble>(nPointsTot);
149
150 // Evaluate V\cdot Grad(u)
151 switch (ndim)
152 {
153 case 1:
154 grad0 = Array<OneD, NekDouble>(fields[0]->GetNpoints());
155 for (int n = 0; n < nConvectiveFields; ++n)
156 {
157 fields[0]->PhysDeriv(inarray[n], grad0);
158 if (m_specHP_dealiasing) // interpolate gradient field
159 {
160 Array<OneD, NekDouble> Outarray(nPointsTot);
161 fields[0]->PhysInterp1DScaled(OneDptscale, grad0, wkSp);
162 Vmath::Vmul(nPointsTot, wkSp, 1, AdvVel[0], 1, Outarray, 1);
163 // Galerkin project solution back to origianl spac
164 timer.Start();
165 fields[0]->PhysGalerkinProjection1DScaled(
166 OneDptscale, Outarray, outarray[n]);
167 timer.Stop();
168 timer.AccumulateRegion("GalerinProject");
169 }
170 else
171 {
172 Vmath::Vmul(nPointsTot, grad0, 1, AdvVel[0], 1, outarray[n],
173 1);
174 }
175 }
176 break;
177 case 2:
178 grad0 = Array<OneD, NekDouble>(nqtot);
179 grad1 = Array<OneD, NekDouble>(nqtot);
180 for (int n = 0; n < nConvectiveFields; ++n)
181 {
182 fields[0]->PhysDeriv(inarray[n], grad0, grad1);
183
184 if (m_specHP_dealiasing) // interpolate gradient field
185 {
186 Array<OneD, NekDouble> Outarray(nPointsTot);
187 fields[0]->PhysInterp1DScaled(OneDptscale, grad0, wkSp);
188 Vmath::Vmul(nPointsTot, wkSp, 1, AdvVel[0], 1, Outarray, 1);
189 timer.Start();
190 fields[0]->PhysInterp1DScaled(OneDptscale, grad1, wkSp);
191 timer.Stop();
192 timer.AccumulateRegion("Interp1DScaled");
193 Vmath::Vvtvp(nPointsTot, wkSp, 1, AdvVel[1], 1, Outarray, 1,
194 Outarray, 1);
195 // Galerkin project solution back to original space
196 timer.Start();
197 fields[0]->PhysGalerkinProjection1DScaled(
198 OneDptscale, Outarray, outarray[n]);
199 timer.Stop();
200 timer.AccumulateRegion("GalerinProject");
201 }
202 else
203 {
204 Vmath::Vmul(nPointsTot, grad0, 1, AdvVel[0], 1, outarray[n],
205 1);
206 Vmath::Vvtvp(nPointsTot, grad1, 1, AdvVel[1], 1,
207 outarray[n], 1, outarray[n], 1);
208 }
209 }
210 break;
211 case 3:
212 if (m_homogen_dealiasing == true && m_specHP_dealiasing == true)
213 {
216 ndim * nConvectiveFields);
217 Array<OneD, Array<OneD, NekDouble>> Outarray(nConvectiveFields);
218 for (int i = 0; i < ndim; i++)
219 {
220 grad[i] = Array<OneD, NekDouble>(nqtot);
221 }
222 for (int i = 0; i < ndim * nConvectiveFields; i++)
223 {
224 gradScaled[i] = Array<OneD, NekDouble>(nPointsTot);
225 }
226 for (int i = 0; i < nConvectiveFields; i++)
227 {
228 Outarray[i] = Array<OneD, NekDouble>(nPointsTot);
229 }
230
231 for (int n = 0; n < nConvectiveFields; n++)
232 {
233 fields[0]->PhysDeriv(inarray[n], grad[0], grad[1], grad[2]);
234 for (int i = 0; i < ndim; i++)
235 {
236 timer.Start();
237 fields[0]->PhysInterp1DScaled(OneDptscale, grad[i],
238 gradScaled[n * ndim + i]);
239 timer.Stop();
240 timer.AccumulateRegion("Interp1DScaled");
241 }
242 }
243
244 fields[0]->DealiasedDotProd(nPointsTot, AdvVel, gradScaled,
245 Outarray);
246
247 timer.Start();
248 for (int n = 0; n < nConvectiveFields; n++)
249 {
250 fields[0]->PhysGalerkinProjection1DScaled(
251 OneDptscale, Outarray[n], outarray[n]);
252 }
253 timer.Stop();
254 timer.AccumulateRegion("GalerinProject");
255 }
256 else if (m_homogen_dealiasing == true &&
257 m_specHP_dealiasing == false)
258 {
260 nConvectiveFields);
261 Array<OneD, Array<OneD, NekDouble>> Outarray(nConvectiveFields);
262 for (int i = 0; i < ndim * nConvectiveFields; i++)
263 {
264 grad[i] = Array<OneD, NekDouble>(nPointsTot);
265 }
266 for (int i = 0; i < nConvectiveFields; i++)
267 {
268 Outarray[i] = Array<OneD, NekDouble>(nPointsTot);
269 }
270
271 for (int n = 0; n < nConvectiveFields; n++)
272 {
273 fields[0]->PhysDeriv(inarray[n], grad[n * ndim + 0],
274 grad[n * ndim + 1],
275 grad[n * ndim + 2]);
276 }
277
278 fields[0]->DealiasedDotProd(nPointsTot, AdvVel, grad, outarray);
279 }
280 else
281 {
282 grad0 = Array<OneD, NekDouble>(nqtot);
283 grad1 = Array<OneD, NekDouble>(nqtot);
284 grad2 = Array<OneD, NekDouble>(nqtot);
285 Array<OneD, NekDouble> tmp = grad2;
286 for (int n = 0; n < nConvectiveFields; ++n)
287 {
288 if (fields[0]->GetWaveSpace() == true &&
289 fields[0]->GetExpType() == MultiRegions::e3DH1D)
290 {
291 fields[0]->HomogeneousBwdTrans(nqtot, inarray[n], tmp);
292 fields[0]->PhysDeriv(tmp, grad0, grad1);
293 // Take d/dz derivative using wave space field
294 fields[0]->PhysDeriv(MultiRegions::DirCartesianMap[2],
295 inarray[n], outarray[n]);
296 fields[0]->HomogeneousBwdTrans(nqtot, outarray[n],
297 grad2);
298 }
299 else if (fields[0]->GetWaveSpace() == true &&
300 fields[0]->GetExpType() == MultiRegions::e3DH2D)
301 {
302 fields[0]->HomogeneousBwdTrans(nqtot, inarray[n], tmp);
303 fields[0]->PhysDeriv(tmp, grad0);
304 // Take d/dy derivative using wave space field
305 fields[0]->PhysDeriv(MultiRegions::DirCartesianMap[1],
306 inarray[n], outarray[n]);
307 fields[0]->HomogeneousBwdTrans(nqtot, outarray[n],
308 grad1);
309 // Take d/dz derivative using wave space field
310 fields[0]->PhysDeriv(MultiRegions::DirCartesianMap[2],
311 inarray[n], outarray[n]);
312 fields[0]->HomogeneousBwdTrans(nqtot, outarray[n],
313 grad2);
314 }
315 else
316 {
317 fields[0]->PhysDeriv(inarray[n], grad0, grad1, grad2);
318 }
319 if (m_specHP_dealiasing) // interpolate spectral/hp gradient
320 // field
321 {
322 Array<OneD, NekDouble> Outarray(nPointsTot);
323 timer.Start();
324 fields[0]->PhysInterp1DScaled(OneDptscale, grad0, wkSp);
325 timer.Stop();
326 timer.AccumulateRegion("Interp1DScaled");
327 Vmath::Vmul(nPointsTot, wkSp, 1, AdvVel[0], 1, Outarray,
328 1);
329
330 timer.Start();
331 fields[0]->PhysInterp1DScaled(OneDptscale, grad1, wkSp);
332 timer.Stop();
333 timer.AccumulateRegion("Interp1DScaled");
334 Vmath::Vvtvp(nPointsTot, wkSp, 1, AdvVel[1], 1,
335 Outarray, 1, Outarray, 1);
336
337 timer.Start();
338 fields[0]->PhysInterp1DScaled(OneDptscale, grad2, wkSp);
339 timer.Stop();
340 timer.AccumulateRegion("Interp1DScaled");
341 Vmath::Vvtvp(nPointsTot, wkSp, 1, AdvVel[2], 1,
342 Outarray, 1, Outarray, 1);
343 timer.Start();
344 fields[0]->PhysGalerkinProjection1DScaled(
345 OneDptscale, Outarray, outarray[n]);
346 timer.Stop();
347 timer.AccumulateRegion("GalerinProject");
348 }
349 else
350 {
351 Vmath::Vmul(nPointsTot, grad0, 1, AdvVel[0], 1,
352 outarray[n], 1);
353 Vmath::Vvtvp(nPointsTot, grad1, 1, AdvVel[1], 1,
354 outarray[n], 1, outarray[n], 1);
355 Vmath::Vvtvp(nPointsTot, grad2, 1, AdvVel[2], 1,
356 outarray[n], 1, outarray[n], 1);
357 }
358
359 if (fields[0]->GetWaveSpace() == true)
360 {
361 fields[0]->HomogeneousFwdTrans(nqtot, outarray[n],
362 outarray[n]);
363 }
364 }
365 }
366 break;
367 default:
368 ASSERTL0(false, "dimension unknown");
369 }
370
371 for (int n = 0; n < nConvectiveFields; ++n)
372 {
373 Vmath::Neg(nqtot, outarray[n], 1);
374 }
375}
376
377} // namespace Nektar
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:208
#define ASSERTL1(condition, msg)
Assert Level 1 – Debugging which is used whether in FULLDEBUG or DEBUG compilation mode....
Definition: ErrorUtil.hpp:242
tKey RegisterCreatorFunction(tKey idKey, CreatorFunction classCreator, std::string pDesc="")
Register a class with the factory.
static std::string RegisterEnumValue(std::string pEnum, std::string pString, int pEnumValue)
Registers an enumeration value.
void AccumulateRegion(std::string, int iolevel=0)
Accumulate elapsed time for a region.
Definition: Timer.cpp:70
void v_InitObject(LibUtilities::SessionReaderSharedPtr pSession, Array< OneD, MultiRegions::ExpListSharedPtr > pFields) override
Initialises the advection object.
static std::string className
Name of class.
static SolverUtils::AdvectionSharedPtr create(std::string)
Creates an instance of this class.
void v_Advect(const int nConvectiveFields, const Array< OneD, MultiRegions::ExpListSharedPtr > &fields, const Array< OneD, Array< OneD, NekDouble > > &advVel, const Array< OneD, Array< OneD, NekDouble > > &inarray, Array< OneD, Array< OneD, NekDouble > > &outarray, const NekDouble &time, const Array< OneD, Array< OneD, NekDouble > > &pFwd=NullNekDoubleArrayOfArray, const Array< OneD, Array< OneD, NekDouble > > &pBwd=NullNekDoubleArrayOfArray) override
Advects a vector field.
static std::string navierStokesAdvectionTypeLookupIds[]
An abstract base class encapsulating the concept of advection of a vector field.
Definition: Advection.h:81
virtual SOLVER_UTILS_EXPORT void v_InitObject(LibUtilities::SessionReaderSharedPtr pSession, Array< OneD, MultiRegions::ExpListSharedPtr > pFields)
Initialises the advection object.
Definition: Advection.cpp:264
std::shared_ptr< SessionReader > SessionReaderSharedPtr
MultiRegions::Direction const DirCartesianMap[]
Definition: ExpList.h:87
AdvectionFactory & GetAdvectionFactory()
Gets the factory for initialising advection objects.
Definition: Advection.cpp:43
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 Neg(int n, T *x, const int incx)
Negate x = -x.
Definition: Vmath.hpp:292
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