Nektar++
Loading...
Searching...
No Matches
SkewSymmetricAdvection.cpp
Go to the documentation of this file.
1///////////////////////////////////////////////////////////////////////////////
2//
3// File: SkewSymmetricAdvection.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
35#include <IncNavierStokesSolver/AdvectionTerms/SkewSymmetricAdvection.h>
36
37namespace Nektar
38{
39std::string SkewSymmetricAdvection::className =
40 SolverUtils::GetAdvectionFactory().RegisterCreatorFunction(
41 "SkewSymmetric", SkewSymmetricAdvection::create, "Skew Symmetric");
42
43/**
44 *
45 */
50
51/**
52 *
53 */
54void SkewSymmetricAdvection::v_InitObject(
55 LibUtilities::SessionReaderSharedPtr pSession,
56 Array<OneD, MultiRegions::ExpListSharedPtr> pFields)
57{
58 Advection::v_InitObject(pSession, pFields);
59
60 m_homogen_dealiasing = pSession->DefinesSolverInfo("dealiasing");
61 pSession->MatchSolverInfo("ModeType", "SingleMode", m_SingleMode, false);
62 pSession->MatchSolverInfo("ModeType", "HalfMode", m_HalfMode, false);
63}
64
65/**
66 *
67 */
68void SkewSymmetricAdvection::v_Advect(
69 const int nConvectiveFields,
70 const Array<OneD, MultiRegions::ExpListSharedPtr> &fields,
71 const Array<OneD, Array<OneD, NekDouble>> &advVel,
72 const Array<OneD, Array<OneD, NekDouble>> &inarray,
73 Array<OneD, Array<OneD, NekDouble>> &outarray,
74 [[maybe_unused]] const NekDouble &time,
75 [[maybe_unused]] const Array<OneD, Array<OneD, NekDouble>> &pFwd,
76 [[maybe_unused]] const Array<OneD, Array<OneD, NekDouble>> &pBwd)
77{
78 // use dimension of Velocity vector to dictate dimension of operation
79 int ndim = advVel.size();
80 int nqtot = fields[0]->GetTotPoints();
81 ASSERTL1(nConvectiveFields == inarray.size(),
82 "Number of convective fields and Inarray are not compatible");
83
84 Array<OneD, Array<OneD, NekDouble>> velocity(ndim);
85 for (int i = 0; i < ndim; ++i)
86 {
87 if (fields[i]->GetWaveSpace() && !m_SingleMode && !m_HalfMode)
88 {
89 velocity[i] = Array<OneD, NekDouble>(nqtot, 0.0);
90 fields[i]->HomogeneousBwdTrans(nqtot, advVel[i], velocity[i]);
91 }
92 else
93 {
94 velocity[i] = advVel[i];
95 }
96 }
97
98 for (int n = 0; n < nConvectiveFields; ++n)
99 {
100 // ToDo: here we should add a check that V has right dimension
101
102 int nPointsTot = fields[0]->GetNpoints();
103 Array<OneD, NekDouble> gradV0, gradV1, gradV2, tmp, Up;
104
105 gradV0 = Array<OneD, NekDouble>(nPointsTot);
106 tmp = Array<OneD, NekDouble>(nPointsTot);
107
108 // Evaluate V\cdot Grad(u)
109 switch (ndim)
110 {
111 case 1:
112 fields[0]->PhysDeriv(inarray[n], gradV0);
113 Vmath::Vmul(nPointsTot, gradV0, 1, velocity[0], 1, outarray[n],
114 1);
115 Vmath::Vmul(nPointsTot, inarray[n], 1, velocity[0], 1, gradV0,
116 1);
117 fields[0]->PhysDeriv(gradV0, tmp);
118 Vmath::Vadd(nPointsTot, tmp, 1, outarray[n], 1, outarray[n], 1);
119 Vmath::Smul(nPointsTot, 0.5, outarray[n], 1, outarray[n], 1);
120 break;
121 case 2:
122 gradV1 = Array<OneD, NekDouble>(nPointsTot);
123 fields[0]->PhysDeriv(inarray[n], gradV0, gradV1);
124 Vmath::Vmul(nPointsTot, gradV0, 1, velocity[0], 1, outarray[n],
125 1);
126 Vmath::Vvtvp(nPointsTot, gradV1, 1, velocity[1], 1, outarray[n],
127 1, outarray[n], 1);
128 Vmath::Vmul(nPointsTot, inarray[n], 1, velocity[0], 1, gradV0,
129 1);
130 Vmath::Vmul(nPointsTot, inarray[n], 1, velocity[1], 1, gradV1,
131 1);
132 fields[0]->PhysDeriv(MultiRegions::DirCartesianMap[0], gradV0,
133 tmp);
134 Vmath::Vadd(nPointsTot, tmp, 1, outarray[n], 1, outarray[n], 1);
135 fields[0]->PhysDeriv(MultiRegions::DirCartesianMap[1], gradV1,
136 tmp);
137 Vmath::Vadd(nPointsTot, tmp, 1, outarray[n], 1, outarray[n], 1);
138 Vmath::Smul(nPointsTot, 0.5, outarray[n], 1, outarray[n], 1);
139 break;
140 case 3:
141 gradV1 = Array<OneD, NekDouble>(nPointsTot);
142 gradV2 = Array<OneD, NekDouble>(nPointsTot);
143
144 fields[0]->PhysDeriv(inarray[n], gradV0, gradV1, gradV2);
145
146 // outarray[n] = 1/2(u*du/dx + v*du/dy + w*du/dz + duu/dx +
147 // duv/dy + duw/dz)
148
149 if (m_homogen_dealiasing == true &&
150 fields[0]->GetWaveSpace() == false)
151 {
152 fields[0]->DealiasedProd(nPointsTot, velocity[0], gradV0,
153 gradV0);
154 fields[0]->DealiasedProd(nPointsTot, velocity[1], gradV1,
155 gradV1);
156 fields[0]->DealiasedProd(nPointsTot, velocity[2], gradV2,
157 gradV2);
158 Vmath::Vadd(nPointsTot, gradV0, 1, gradV1, 1, outarray[n],
159 1);
160 Vmath::Vadd(nPointsTot, gradV2, 1, outarray[n], 1,
161 outarray[n], 1);
162 fields[0]->DealiasedProd(nPointsTot, inarray[n],
163 velocity[0], gradV0);
164 fields[0]->DealiasedProd(nPointsTot, inarray[n],
165 velocity[1], gradV1);
166 fields[0]->DealiasedProd(nPointsTot, inarray[n],
167 velocity[2], gradV2);
168 fields[0]->PhysDeriv(MultiRegions::DirCartesianMap[0],
169 gradV0, tmp);
170 Vmath::Vadd(nPointsTot, tmp, 1, outarray[n], 1, outarray[n],
171 1);
172 fields[0]->PhysDeriv(MultiRegions::DirCartesianMap[1],
173 gradV1, tmp);
174 Vmath::Vadd(nPointsTot, tmp, 1, outarray[n], 1, outarray[n],
175 1);
176 fields[0]->PhysDeriv(MultiRegions::DirCartesianMap[2],
177 gradV2, tmp);
178 Vmath::Vadd(nPointsTot, tmp, 1, outarray[n], 1, outarray[n],
179 1);
180 Vmath::Smul(nPointsTot, 0.5, outarray[n], 1, outarray[n],
181 1);
182 }
183 else if (fields[0]->GetWaveSpace() == true &&
184 m_homogen_dealiasing == false)
185 {
186 Up = Array<OneD, NekDouble>(nPointsTot);
187 // vector reused to avoid even more memory requirements
188 // names may be misleading
189 fields[0]->HomogeneousBwdTrans(nPointsTot, gradV0, tmp);
190 Vmath::Vmul(nPointsTot, tmp, 1, velocity[0], 1, outarray[n],
191 1); // + u*du/dx
192 fields[0]->HomogeneousBwdTrans(nPointsTot, gradV1, tmp);
193 Vmath::Vvtvp(nPointsTot, tmp, 1, velocity[1], 1,
194 outarray[n], 1, outarray[n], 1); // + v*du/dy
195 fields[0]->HomogeneousBwdTrans(nPointsTot, gradV2, tmp);
196 Vmath::Vvtvp(nPointsTot, tmp, 1, velocity[2], 1,
197 outarray[n], 1, outarray[n], 1); // + w*du/dz
198
199 fields[0]->HomogeneousBwdTrans(nPointsTot, inarray[n], Up);
200 Vmath::Vmul(nPointsTot, Up, 1, velocity[0], 1, gradV0, 1);
201 Vmath::Vmul(nPointsTot, Up, 1, velocity[1], 1, gradV1, 1);
202 Vmath::Vmul(nPointsTot, Up, 1, velocity[2], 1, gradV2, 1);
203
204 fields[0]->SetWaveSpace(false);
205 fields[0]->PhysDeriv(MultiRegions::DirCartesianMap[0],
206 gradV0, tmp); // duu/dx
207 Vmath::Vadd(nPointsTot, tmp, 1, outarray[n], 1, outarray[n],
208 1);
209 fields[0]->PhysDeriv(MultiRegions::DirCartesianMap[1],
210 gradV1, tmp); // duv/dy
211 Vmath::Vadd(nPointsTot, tmp, 1, outarray[n], 1, outarray[n],
212 1);
213 fields[0]->PhysDeriv(MultiRegions::DirCartesianMap[2],
214 gradV2, tmp); // duw/dz
215 Vmath::Vadd(nPointsTot, tmp, 1, outarray[n], 1, outarray[n],
216 1);
217 fields[0]->SetWaveSpace(true);
218
219 Vmath::Smul(nPointsTot, 0.5, outarray[n], 1, tmp, 1);
220 fields[0]->HomogeneousFwdTrans(nPointsTot, tmp,
221 outarray[n]);
222 }
223 else if (fields[0]->GetWaveSpace() == false &&
224 m_homogen_dealiasing == false)
225 {
226 Vmath::Vmul(nPointsTot, gradV0, 1, velocity[0], 1,
227 outarray[n], 1);
228 Vmath::Vvtvp(nPointsTot, gradV1, 1, velocity[1], 1,
229 outarray[n], 1, outarray[n], 1);
230 Vmath::Vvtvp(nPointsTot, gradV2, 1, velocity[2], 1,
231 outarray[n], 1, outarray[n], 1);
232 Vmath::Vmul(nPointsTot, inarray[n], 1, velocity[0], 1,
233 gradV0, 1);
234 Vmath::Vmul(nPointsTot, inarray[n], 1, velocity[1], 1,
235 gradV1, 1);
236 Vmath::Vmul(nPointsTot, inarray[n], 1, velocity[2], 1,
237 gradV2, 1);
238 fields[0]->PhysDeriv(MultiRegions::DirCartesianMap[0],
239 gradV0, tmp);
240 Vmath::Vadd(nPointsTot, tmp, 1, outarray[n], 1, outarray[n],
241 1);
242 fields[0]->PhysDeriv(MultiRegions::DirCartesianMap[1],
243 gradV1, tmp);
244 Vmath::Vadd(nPointsTot, tmp, 1, outarray[n], 1, outarray[n],
245 1);
246 fields[0]->PhysDeriv(MultiRegions::DirCartesianMap[2],
247 gradV2, tmp);
248 Vmath::Vadd(nPointsTot, tmp, 1, outarray[n], 1, outarray[n],
249 1);
250 Vmath::Smul(nPointsTot, 0.5, outarray[n], 1, outarray[n],
251 1);
252 }
253 else
254 {
255 ASSERTL0(false,
256 "Dealiasing is not allowed in combination "
257 "with the Skew-Symmetric advection form for "
258 "efficiency reasons.");
259 }
260 break;
261 default:
262 ASSERTL0(false, "dimension unknown");
263 }
264
265 Vmath::Neg(nqtot, outarray[n], 1);
266 }
267}
268
269} // namespace Nektar
ASSERTL0
#define ASSERTL0(condition, msg)
Definition ErrorUtil.hpp:216
ASSERTL1
#define ASSERTL1(condition, msg)
Assert Level 1 – Debugging which is used whether in FULLDEBUG or DEBUG compilation mode....
Definition ErrorUtil.hpp:250
Nektar::LibUtilities::NekFactory::RegisterCreatorFunction
tKey RegisterCreatorFunction(tKey idKey, CreatorFunction classCreator, std::string pDesc="")
Register a class with the factory.
Definition BasicUtils/NekFactory.hpp:196
Nektar::SkewSymmetricAdvection::create
static SolverUtils::AdvectionSharedPtr create(std::string)
Creates an instance of this class.
Definition SkewSymmetricAdvection.h:50
Nektar::SkewSymmetricAdvection::v_Advect
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.
Definition SkewSymmetricAdvection.cpp:68
Nektar::SkewSymmetricAdvection::v_InitObject
void v_InitObject(LibUtilities::SessionReaderSharedPtr pSession, Array< OneD, MultiRegions::ExpListSharedPtr > pFields) override
Initialises the advection object.
Definition SkewSymmetricAdvection.cpp:54
Nektar::SkewSymmetricAdvection::className
static std::string className
Name of class.
Definition SkewSymmetricAdvection.h:56
Nektar::SolverUtils::Advection
An abstract base class encapsulating the concept of advection of a vector field.
Definition Advection.h:81
Nektar::SolverUtils::Advection::v_InitObject
virtual SOLVER_UTILS_EXPORT void v_InitObject(LibUtilities::SessionReaderSharedPtr pSession, Array< OneD, MultiRegions::ExpListSharedPtr > pFields)
Initialises the advection object.
Definition Advection.cpp:264
Nektar::LibUtilities::SessionReaderSharedPtr
std::shared_ptr< SessionReader > SessionReaderSharedPtr
Definition SessionReader.h:119
Nektar::MultiRegions::DirCartesianMap
MultiRegions::Direction const DirCartesianMap[]
Definition ExpList.h:87
Nektar::SolverUtils::GetAdvectionFactory
AdvectionFactory & GetAdvectionFactory()
Gets the factory for initialising advection objects.
Definition Advection.cpp:43
Vmath::Vmul
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
Vmath::Neg
void Neg(int n, T *x, const int incx)
Negate x = -x.
Definition Vmath.hpp:292
Vmath::Vvtvp
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
Vmath::Vadd
void Vadd(int n, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
Add vector z = x+y.
Definition Vmath.hpp:180
Vmath::Smul
void Smul(int n, const T alpha, const T *x, const int incx, T *y, const int incy)
Scalar multiply y = alpha*x.
Definition Vmath.hpp:100
Generated on Sun Nov 9 2025 20:17:01 for Nektar++ by doxygen 1.9.8