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