Nektar++
AlternateSkewAdvection.cpp
Go to the documentation of this file.
1///////////////////////////////////////////////////////////////////////////////
2//
3// File: AlternateSkewAdvection.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
36
37namespace Nektar
38{
39
42 "AlternateSkew", AlternateSkewAdvection::create,
43 "Alternating Skew Symmetric");
44
45/**
46 * Constructor. Creates ...
47 *
48 * \param
49 * \param
50 */
52{
53}
54
57 [[maybe_unused]] Array<OneD, MultiRegions::ExpListSharedPtr> fields)
58{
59 pSession->MatchSolverInfo("ModeType", "SingleMode", m_SingleMode, false);
60 pSession->MatchSolverInfo("ModeType", "HalfMode", m_HalfMode, false);
61}
62
64 const int nConvectiveFields,
66 const Array<OneD, Array<OneD, NekDouble>> &advVel,
67 const Array<OneD, Array<OneD, NekDouble>> &inarray,
69 [[maybe_unused]] const NekDouble &time,
70 [[maybe_unused]] const Array<OneD, Array<OneD, NekDouble>> &pFwd,
71 [[maybe_unused]] const Array<OneD, Array<OneD, NekDouble>> &pBwd)
72{
73 // use dimension of Velocity vector to dictate dimension of operation
74 int ndim = advVel.size();
75 int nPointsTot = fields[0]->GetNpoints();
77 for (int i = 0; i < ndim; ++i)
78 {
79 if (fields[i]->GetWaveSpace() && !m_SingleMode && !m_HalfMode)
80 {
81 velocity[i] = Array<OneD, NekDouble>(nPointsTot, 0.0);
82 fields[i]->HomogeneousBwdTrans(nPointsTot, advVel[i], velocity[i]);
83 }
84 else
85 {
86 velocity[i] = advVel[i];
87 }
88 }
89 for (int n = 0; n < nConvectiveFields; ++n)
90 {
91 // ToDo: here we should add a check that V has right dimension
92 Array<OneD, NekDouble> gradV0, gradV1, gradV2, tmp, Up;
93
94 gradV0 = Array<OneD, NekDouble>(nPointsTot);
95 tmp = Array<OneD, NekDouble>(nPointsTot);
96
97 // Evaluate V\cdot Grad(u)
98 switch (ndim)
99 {
100 case 1:
101 if (m_advectioncalls % 2 == 0)
102 {
103 fields[0]->PhysDeriv(inarray[n], gradV0);
104 Vmath::Vmul(nPointsTot, gradV0, 1, velocity[0], 1,
105 outarray[n], 1);
106 }
107 else
108 {
109 Vmath::Vmul(nPointsTot, inarray[n], 1, velocity[0], 1,
110 gradV0, 1);
111 fields[0]->PhysDeriv(gradV0, outarray[n]);
112 }
113 Vmath::Smul(nPointsTot, 0.5, outarray[n], 1, outarray[n],
114 1); // must be mult by 0.5????
115 break;
116 case 2:
117 gradV1 = Array<OneD, NekDouble>(nPointsTot);
118 if (m_advectioncalls % 2 == 0)
119 {
120 fields[0]->PhysDeriv(inarray[n], gradV0, gradV1);
121 Vmath::Vmul(nPointsTot, gradV0, 1, velocity[0], 1,
122 outarray[n], 1);
123 Vmath::Vvtvp(nPointsTot, gradV1, 1, velocity[1], 1,
124 outarray[n], 1, outarray[n], 1);
125 }
126 else
127 {
128 Vmath::Vmul(nPointsTot, inarray[n], 1, velocity[0], 1,
129 gradV0, 1);
130 Vmath::Vmul(nPointsTot, inarray[n], 1, velocity[1], 1,
131 gradV1, 1);
132 fields[0]->PhysDeriv(MultiRegions::DirCartesianMap[0],
133 gradV0, outarray[n]);
134 fields[0]->PhysDeriv(MultiRegions::DirCartesianMap[1],
135 gradV1, tmp);
136 Vmath::Vadd(nPointsTot, tmp, 1, outarray[n], 1, outarray[n],
137 1);
138 }
139 Vmath::Smul(nPointsTot, 1.0, outarray[n], 1, outarray[n],
140 1); // must be mult by 0.5????
141 break;
142 case 3:
143 gradV1 = Array<OneD, NekDouble>(nPointsTot);
144 gradV2 = Array<OneD, NekDouble>(nPointsTot);
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 (fields[0]->GetWaveSpace() == true)
150 {
151 if (m_advectioncalls % 2 == 0)
152 {
153 // vector reused to avoid even more memory requirements
154 // names may be misleading
155 fields[0]->PhysDeriv(inarray[n], gradV0, gradV1,
156 gradV2);
157 fields[0]->HomogeneousBwdTrans(nPointsTot, gradV0, tmp);
158 Vmath::Vmul(nPointsTot, tmp, 1, velocity[0], 1,
159 outarray[n], 1); // + u*du/dx
160 fields[0]->HomogeneousBwdTrans(nPointsTot, gradV1, tmp);
161 Vmath::Vvtvp(nPointsTot, tmp, 1, velocity[1], 1,
162 outarray[n], 1, outarray[n],
163 1); // + v*du/dy
164 fields[0]->HomogeneousBwdTrans(nPointsTot, gradV2, tmp);
165 Vmath::Vvtvp(nPointsTot, tmp, 1, velocity[2], 1,
166 outarray[n], 1, outarray[n],
167 1); // + w*du/dz
168 }
169 else
170 {
171 Up = Array<OneD, NekDouble>(nPointsTot);
172 fields[0]->HomogeneousBwdTrans(nPointsTot, inarray[n],
173 Up);
174 Vmath::Vmul(nPointsTot, Up, 1, velocity[0], 1, gradV0,
175 1);
176 Vmath::Vmul(nPointsTot, Up, 1, velocity[1], 1, gradV1,
177 1);
178 Vmath::Vmul(nPointsTot, Up, 1, velocity[2], 1, gradV2,
179 1);
180
181 fields[0]->SetWaveSpace(false);
182 fields[0]->PhysDeriv(MultiRegions::DirCartesianMap[0],
183 gradV0, outarray[n]); // duu/dx
184 fields[0]->PhysDeriv(MultiRegions::DirCartesianMap[1],
185 gradV1, tmp); // duv/dy
186 Vmath::Vadd(nPointsTot, tmp, 1, outarray[n], 1,
187 outarray[n], 1);
188 fields[0]->PhysDeriv(MultiRegions::DirCartesianMap[2],
189 gradV2, tmp); // duw/dz
190 Vmath::Vadd(nPointsTot, tmp, 1, outarray[n], 1,
191 outarray[n], 1);
192 fields[0]->SetWaveSpace(true);
193 }
194
195 Vmath::Smul(nPointsTot, 1.0, outarray[n], 1, tmp,
196 1); // must be mult by 0.5????
197 fields[0]->HomogeneousFwdTrans(nPointsTot, tmp,
198 outarray[n]);
199 }
200 else
201 {
202 if (m_advectioncalls % 2 == 0)
203 {
204 fields[0]->PhysDeriv(inarray[n], gradV0, gradV1,
205 gradV2);
206 Vmath::Vmul(nPointsTot, gradV0, 1, velocity[0], 1,
207 outarray[n], 1);
208 Vmath::Vvtvp(nPointsTot, gradV1, 1, velocity[1], 1,
209 outarray[n], 1, outarray[n], 1);
210 Vmath::Vvtvp(nPointsTot, gradV2, 1, velocity[2], 1,
211 outarray[n], 1, outarray[n], 1);
212 }
213 else
214 {
215 Vmath::Vmul(nPointsTot, inarray[n], 1, velocity[0], 1,
216 gradV0, 1);
217 Vmath::Vmul(nPointsTot, inarray[n], 1, velocity[1], 1,
218 gradV1, 1);
219 Vmath::Vmul(nPointsTot, inarray[n], 1, velocity[2], 1,
220 gradV2, 1);
221 fields[0]->PhysDeriv(MultiRegions::DirCartesianMap[0],
222 gradV0, outarray[n]);
223 fields[0]->PhysDeriv(MultiRegions::DirCartesianMap[1],
224 gradV1, tmp);
225 Vmath::Vadd(nPointsTot, tmp, 1, outarray[n], 1,
226 outarray[n], 1);
227 fields[0]->PhysDeriv(MultiRegions::DirCartesianMap[2],
228 gradV2, tmp);
229 Vmath::Vadd(nPointsTot, tmp, 1, outarray[n], 1,
230 outarray[n], 1);
231 }
232 Vmath::Smul(nPointsTot, 1.0, outarray[n], 1, outarray[n],
233 1); // must be mult by 0.5????
234 }
235 break;
236 default:
237 ASSERTL0(false, "dimension unknown");
238 }
239 }
240}
241
242} // namespace Nektar
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:208
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.
void v_InitObject(LibUtilities::SessionReaderSharedPtr pSession, Array< OneD, MultiRegions::ExpListSharedPtr > pFields) override
Initialises the advection object.
static std::string className
Name of class.
tKey RegisterCreatorFunction(tKey idKey, CreatorFunction classCreator, std::string pDesc="")
Register a class with the factory.
An abstract base class encapsulating the concept of advection of a vector field.
Definition: Advection.h:81
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 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 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
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