Advects a vector field.
86{
87
88 int ndim = advVel.size();
89 int nqtot = fields[0]->GetTotPoints();
90 ASSERTL1(nConvectiveFields == inarray.size(),
91 "Number of convective fields and Inarray are not compatible");
92
93 Array<OneD, Array<OneD, NekDouble>> velocity(ndim);
94 for (int i = 0; i < ndim; ++i)
95 {
97 {
98 velocity[i] = Array<OneD, NekDouble>(nqtot, 0.0);
99 fields[i]->HomogeneousBwdTrans(nqtot, advVel[i], velocity[i]);
100 }
101 else
102 {
103 velocity[i] = advVel[i];
104 }
105 }
106
107 for (int n = 0; n < nConvectiveFields; ++n)
108 {
109
110
111 int nPointsTot = fields[0]->GetNpoints();
112 Array<OneD, NekDouble> gradV0, gradV1, gradV2, tmp, Up;
113
114 gradV0 = Array<OneD, NekDouble>(nPointsTot);
115 tmp = Array<OneD, NekDouble>(nPointsTot);
116
117
118 switch (ndim)
119 {
120 case 1:
121 fields[0]->PhysDeriv(inarray[n], gradV0);
122 Vmath::Vmul(nPointsTot, gradV0, 1, velocity[0], 1, outarray[n],
123 1);
124 Vmath::Vmul(nPointsTot, inarray[n], 1, velocity[0], 1, gradV0,
125 1);
126 fields[0]->PhysDeriv(gradV0, tmp);
127 Vmath::Vadd(nPointsTot, tmp, 1, outarray[n], 1, outarray[n], 1);
128 Vmath::Smul(nPointsTot, 0.5, outarray[n], 1, outarray[n], 1);
129 break;
130 case 2:
131 gradV1 = Array<OneD, NekDouble>(nPointsTot);
132 fields[0]->PhysDeriv(inarray[n], gradV0, gradV1);
133 Vmath::Vmul(nPointsTot, gradV0, 1, velocity[0], 1, outarray[n],
134 1);
135 Vmath::Vvtvp(nPointsTot, gradV1, 1, velocity[1], 1, outarray[n],
136 1, outarray[n], 1);
137 Vmath::Vmul(nPointsTot, inarray[n], 1, velocity[0], 1, gradV0,
138 1);
139 Vmath::Vmul(nPointsTot, inarray[n], 1, velocity[1], 1, gradV1,
140 1);
142 tmp);
143 Vmath::Vadd(nPointsTot, tmp, 1, outarray[n], 1, outarray[n], 1);
145 tmp);
146 Vmath::Vadd(nPointsTot, tmp, 1, outarray[n], 1, outarray[n], 1);
147 Vmath::Smul(nPointsTot, 0.5, outarray[n], 1, outarray[n], 1);
148 break;
149 case 3:
150 gradV1 = Array<OneD, NekDouble>(nPointsTot);
151 gradV2 = Array<OneD, NekDouble>(nPointsTot);
152
153 fields[0]->PhysDeriv(inarray[n], gradV0, gradV1, gradV2);
154
155
156
157
159 fields[0]->GetWaveSpace() == false)
160 {
161 fields[0]->DealiasedProd(nPointsTot, velocity[0], gradV0,
162 gradV0);
163 fields[0]->DealiasedProd(nPointsTot, velocity[1], gradV1,
164 gradV1);
165 fields[0]->DealiasedProd(nPointsTot, velocity[2], gradV2,
166 gradV2);
167 Vmath::Vadd(nPointsTot, gradV0, 1, gradV1, 1, outarray[n],
168 1);
170 outarray[n], 1);
171 fields[0]->DealiasedProd(nPointsTot, inarray[n],
172 velocity[0], gradV0);
173 fields[0]->DealiasedProd(nPointsTot, inarray[n],
174 velocity[1], gradV1);
175 fields[0]->DealiasedProd(nPointsTot, inarray[n],
176 velocity[2], gradV2);
178 gradV0, tmp);
179 Vmath::Vadd(nPointsTot, tmp, 1, outarray[n], 1, outarray[n],
180 1);
182 gradV1, tmp);
183 Vmath::Vadd(nPointsTot, tmp, 1, outarray[n], 1, outarray[n],
184 1);
186 gradV2, tmp);
187 Vmath::Vadd(nPointsTot, tmp, 1, outarray[n], 1, outarray[n],
188 1);
189 Vmath::Smul(nPointsTot, 0.5, outarray[n], 1, outarray[n],
190 1);
191 }
192 else if (fields[0]->GetWaveSpace() == true &&
194 {
195 Up = Array<OneD, NekDouble>(nPointsTot);
196
197
198 fields[0]->HomogeneousBwdTrans(nPointsTot, gradV0, tmp);
199 Vmath::Vmul(nPointsTot, tmp, 1, velocity[0], 1, outarray[n],
200 1);
201 fields[0]->HomogeneousBwdTrans(nPointsTot, gradV1, tmp);
203 outarray[n], 1, outarray[n], 1);
204 fields[0]->HomogeneousBwdTrans(nPointsTot, gradV2, tmp);
206 outarray[n], 1, outarray[n], 1);
207
208 fields[0]->HomogeneousBwdTrans(nPointsTot, inarray[n], Up);
209 Vmath::Vmul(nPointsTot, Up, 1, velocity[0], 1, gradV0, 1);
210 Vmath::Vmul(nPointsTot, Up, 1, velocity[1], 1, gradV1, 1);
211 Vmath::Vmul(nPointsTot, Up, 1, velocity[2], 1, gradV2, 1);
212
213 fields[0]->SetWaveSpace(false);
215 gradV0, tmp);
216 Vmath::Vadd(nPointsTot, tmp, 1, outarray[n], 1, outarray[n],
217 1);
219 gradV1, tmp);
220 Vmath::Vadd(nPointsTot, tmp, 1, outarray[n], 1, outarray[n],
221 1);
223 gradV2, tmp);
224 Vmath::Vadd(nPointsTot, tmp, 1, outarray[n], 1, outarray[n],
225 1);
226 fields[0]->SetWaveSpace(true);
227
228 Vmath::Smul(nPointsTot, 0.5, outarray[n], 1, tmp, 1);
229 fields[0]->HomogeneousFwdTrans(nPointsTot, tmp,
230 outarray[n]);
231 }
232 else if (fields[0]->GetWaveSpace() == false &&
234 {
236 outarray[n], 1);
238 outarray[n], 1, outarray[n], 1);
240 outarray[n], 1, outarray[n], 1);
241 Vmath::Vmul(nPointsTot, inarray[n], 1, velocity[0], 1,
242 gradV0, 1);
243 Vmath::Vmul(nPointsTot, inarray[n], 1, velocity[1], 1,
244 gradV1, 1);
245 Vmath::Vmul(nPointsTot, inarray[n], 1, velocity[2], 1,
246 gradV2, 1);
248 gradV0, tmp);
249 Vmath::Vadd(nPointsTot, tmp, 1, outarray[n], 1, outarray[n],
250 1);
252 gradV1, tmp);
253 Vmath::Vadd(nPointsTot, tmp, 1, outarray[n], 1, outarray[n],
254 1);
256 gradV2, tmp);
257 Vmath::Vadd(nPointsTot, tmp, 1, outarray[n], 1, outarray[n],
258 1);
259 Vmath::Smul(nPointsTot, 0.5, outarray[n], 1, outarray[n],
260 1);
261 }
262 else
263 {
265 "Dealiasing is not allowed in combination "
266 "with the Skew-Symmetric advection form for "
267 "efficiency reasons.");
268 }
269 break;
270 default:
271 ASSERTL0(
false,
"dimension unknown");
272 }
273
275 }
276}
#define ASSERTL0(condition, msg)
#define ASSERTL1(condition, msg)
Assert Level 1 – Debugging which is used whether in FULLDEBUG or DEBUG compilation mode....
bool m_homogen_dealiasing
MultiRegions::Direction const DirCartesianMap[]
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.
void Neg(int n, T *x, const int incx)
Negate x = -x.
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
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.
void Smul(int n, const T alpha, const T *x, const int incx, T *y, const int incy)
Scalar multiply y = alpha*x.