Advects a vector field.
94{
95 boost::ignore_unused(time, pFwd, pBwd);
96
97 int nqtot = fields[0]->GetTotPoints();
98 ASSERTL1(nConvectiveFields == inarray.size(),
99 "Number of convective fields and Inarray are not compatible");
100
101
102 int ndim = advVel.size();
103 Array<OneD, Array<OneD, NekDouble>> AdvVel(advVel.size());
104
105 Array<OneD, Array<OneD, NekDouble>> velocity(ndim);
106
107 LibUtilities::Timer timer;
108
109 for (int i = 0; i < ndim; ++i)
110 {
113 {
114 velocity[i] = Array<OneD, NekDouble>(nqtot, 0.0);
115 fields[i]->HomogeneousBwdTrans(nqtot, advVel[i], velocity[i]);
116 }
117 else
118 {
119 velocity[i] = advVel[i];
120 }
121 }
122
123 int nPointsTot = fields[0]->GetNpoints();
124 Array<OneD, NekDouble> grad0, grad1, grad2, wkSp;
125
127
129 {
130
131 nPointsTot = fields[0]->Get1DScaledTotPoints(OneDptscale);
132 }
133
134
136 {
137 for (int i = 0; i < ndim; ++i)
138 {
139 AdvVel[i] = Array<OneD, NekDouble>(nPointsTot);
140
141 timer.Start();
142 fields[0]->PhysInterp1DScaled(OneDptscale, velocity[i], AdvVel[i]);
143 timer.Stop();
144 timer.AccumulateRegion("Interp1DScaled");
145 }
146 }
147 else
148 {
149 for (int i = 0; i < ndim; ++i)
150 {
151 AdvVel[i] = velocity[i];
152 }
153 }
154
155 wkSp = Array<OneD, NekDouble>(nPointsTot);
156
157
158 switch (ndim)
159 {
160 case 1:
161 grad0 = Array<OneD, NekDouble>(fields[0]->GetNpoints());
162 for (int n = 0; n < nConvectiveFields; ++n)
163 {
164 fields[0]->PhysDeriv(inarray[n], grad0);
166 {
167 Array<OneD, NekDouble> Outarray(nPointsTot);
168 fields[0]->PhysInterp1DScaled(OneDptscale, grad0, wkSp);
169 Vmath::Vmul(nPointsTot, wkSp, 1, AdvVel[0], 1, Outarray, 1);
170
171 timer.Start();
172 fields[0]->PhysGalerkinProjection1DScaled(
173 OneDptscale, Outarray, outarray[n]);
174 timer.Stop();
175 timer.AccumulateRegion("GalerinProject");
176 }
177 else
178 {
179 Vmath::Vmul(nPointsTot, grad0, 1, AdvVel[0], 1, outarray[n],
180 1);
181 }
182 }
183 break;
184 case 2:
185 grad0 = Array<OneD, NekDouble>(nqtot);
186 grad1 = Array<OneD, NekDouble>(nqtot);
187 for (int n = 0; n < nConvectiveFields; ++n)
188 {
189 fields[0]->PhysDeriv(inarray[n], grad0, grad1);
190
192 {
193 Array<OneD, NekDouble> Outarray(nPointsTot);
194 fields[0]->PhysInterp1DScaled(OneDptscale, grad0, wkSp);
195 Vmath::Vmul(nPointsTot, wkSp, 1, AdvVel[0], 1, Outarray, 1);
196 timer.Start();
197 fields[0]->PhysInterp1DScaled(OneDptscale, grad1, wkSp);
198 timer.Stop();
199 timer.AccumulateRegion("Interp1DScaled");
200 Vmath::Vvtvp(nPointsTot, wkSp, 1, AdvVel[1], 1, Outarray, 1,
201 Outarray, 1);
202
203 timer.Start();
204 fields[0]->PhysGalerkinProjection1DScaled(
205 OneDptscale, Outarray, outarray[n]);
206 timer.Stop();
207 timer.AccumulateRegion("GalerinProject");
208 }
209 else
210 {
211 Vmath::Vmul(nPointsTot, grad0, 1, AdvVel[0], 1, outarray[n],
212 1);
214 outarray[n], 1, outarray[n], 1);
215 }
216 }
217 break;
218 case 3:
220 {
221 Array<OneD, Array<OneD, NekDouble>> grad(ndim);
222 Array<OneD, Array<OneD, NekDouble>> gradScaled(
223 ndim * nConvectiveFields);
224 Array<OneD, Array<OneD, NekDouble>> Outarray(nConvectiveFields);
225 for (int i = 0; i < ndim; i++)
226 {
227 grad[i] = Array<OneD, NekDouble>(nqtot);
228 }
229 for (int i = 0; i < ndim * nConvectiveFields; i++)
230 {
231 gradScaled[i] = Array<OneD, NekDouble>(nPointsTot);
232 }
233 for (int i = 0; i < nConvectiveFields; i++)
234 {
235 Outarray[i] = Array<OneD, NekDouble>(nPointsTot);
236 }
237
238 for (int n = 0; n < nConvectiveFields; n++)
239 {
240 fields[0]->PhysDeriv(inarray[n], grad[0], grad[1], grad[2]);
241 for (int i = 0; i < ndim; i++)
242 {
243 timer.Start();
244 fields[0]->PhysInterp1DScaled(OneDptscale, grad[i],
245 gradScaled[n * ndim + i]);
246 timer.Stop();
247 timer.AccumulateRegion("Interp1DScaled");
248 }
249 }
250
251 fields[0]->DealiasedDotProd(nPointsTot, AdvVel, gradScaled,
252 Outarray);
253
254 timer.Start();
255 for (int n = 0; n < nConvectiveFields; n++)
256 {
257 fields[0]->PhysGalerkinProjection1DScaled(
258 OneDptscale, Outarray[n], outarray[n]);
259 }
260 timer.Stop();
261 timer.AccumulateRegion("GalerinProject");
262 }
265 {
266 Array<OneD, Array<OneD, NekDouble>> grad(ndim *
267 nConvectiveFields);
268 Array<OneD, Array<OneD, NekDouble>> Outarray(nConvectiveFields);
269 for (int i = 0; i < ndim * nConvectiveFields; i++)
270 {
271 grad[i] = Array<OneD, NekDouble>(nPointsTot);
272 }
273 for (int i = 0; i < nConvectiveFields; i++)
274 {
275 Outarray[i] = Array<OneD, NekDouble>(nPointsTot);
276 }
277
278 for (int n = 0; n < nConvectiveFields; n++)
279 {
280 fields[0]->PhysDeriv(inarray[n], grad[n * ndim + 0],
281 grad[n * ndim + 1],
282 grad[n * ndim + 2]);
283 }
284
285 fields[0]->DealiasedDotProd(nPointsTot, AdvVel, grad, outarray);
286 }
287 else
288 {
289 grad0 = Array<OneD, NekDouble>(nqtot);
290 grad1 = Array<OneD, NekDouble>(nqtot);
291 grad2 = Array<OneD, NekDouble>(nqtot);
292 Array<OneD, NekDouble> tmp = grad2;
293 for (int n = 0; n < nConvectiveFields; ++n)
294 {
295 if (fields[0]->GetWaveSpace() == true &&
297 {
298 fields[0]->HomogeneousBwdTrans(nqtot, inarray[n], tmp);
299 fields[0]->PhysDeriv(tmp, grad0, grad1);
300
302 inarray[n], outarray[n]);
303 fields[0]->HomogeneousBwdTrans(nqtot, outarray[n],
304 grad2);
305 }
306 else if (fields[0]->GetWaveSpace() == true &&
308 {
309 fields[0]->HomogeneousBwdTrans(nqtot, inarray[n], tmp);
310 fields[0]->PhysDeriv(tmp, grad0);
311
313 inarray[n], outarray[n]);
314 fields[0]->HomogeneousBwdTrans(nqtot, outarray[n],
315 grad1);
316
318 inarray[n], outarray[n]);
319 fields[0]->HomogeneousBwdTrans(nqtot, outarray[n],
320 grad2);
321 }
322 else
323 {
324 fields[0]->PhysDeriv(inarray[n], grad0, grad1, grad2);
325 }
327
328 {
329 Array<OneD, NekDouble> Outarray(nPointsTot);
330 timer.Start();
331 fields[0]->PhysInterp1DScaled(OneDptscale, grad0, wkSp);
332 timer.Stop();
333 timer.AccumulateRegion("Interp1DScaled");
334 Vmath::Vmul(nPointsTot, wkSp, 1, AdvVel[0], 1, Outarray,
335 1);
336
337 timer.Start();
338 fields[0]->PhysInterp1DScaled(OneDptscale, grad1, wkSp);
339 timer.Stop();
340 timer.AccumulateRegion("Interp1DScaled");
342 Outarray, 1, Outarray, 1);
343
344 timer.Start();
345 fields[0]->PhysInterp1DScaled(OneDptscale, grad2, wkSp);
346 timer.Stop();
347 timer.AccumulateRegion("Interp1DScaled");
349 Outarray, 1, Outarray, 1);
350 timer.Start();
351 fields[0]->PhysGalerkinProjection1DScaled(
352 OneDptscale, Outarray, outarray[n]);
353 timer.Stop();
354 timer.AccumulateRegion("GalerinProject");
355 }
356 else
357 {
359 outarray[n], 1);
361 outarray[n], 1, outarray[n], 1);
363 outarray[n], 1, outarray[n], 1);
364 }
365
366 if (fields[0]->GetWaveSpace() == true)
367 {
368 fields[0]->HomogeneousFwdTrans(nqtot, outarray[n],
369 outarray[n]);
370 }
371 }
372 }
373 break;
374 default:
375 ASSERTL0(
false,
"dimension unknown");
376 }
377
378 for (int n = 0; n < nConvectiveFields; ++n)
379 {
381 }
382}
#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