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