Advects a vector field.
95 int nqtot = fields[0]->GetTotPoints();
96 ASSERTL1(nConvectiveFields == inarray.size(),
97 "Number of convective fields and Inarray are not compatible");
100 int ndim = advVel.size();
101 Array<OneD, Array<OneD, NekDouble>> AdvVel(advVel.size());
103 Array<OneD, Array<OneD, NekDouble>> velocity(ndim);
105 LibUtilities::Timer timer;
107 for (
int i = 0; i < ndim; ++i)
112 velocity[i] = Array<OneD, NekDouble>(nqtot, 0.0);
113 fields[i]->HomogeneousBwdTrans(advVel[i], velocity[i]);
117 velocity[i] = advVel[i];
121 int nPointsTot = fields[0]->GetNpoints();
122 Array<OneD, NekDouble> grad0, grad1, grad2, wkSp;
129 nPointsTot = fields[0]->Get1DScaledTotPoints(OneDptscale);
135 for (
int i = 0; i < ndim; ++i)
137 AdvVel[i] = Array<OneD, NekDouble>(nPointsTot);
140 fields[0]->PhysInterp1DScaled(OneDptscale, velocity[i], AdvVel[i]);
142 timer.AccumulateRegion(
"Interp1DScaled");
147 for (
int i = 0; i < ndim; ++i)
149 AdvVel[i] = velocity[i];
153 wkSp = Array<OneD, NekDouble>(nPointsTot);
159 grad0 = Array<OneD, NekDouble>(fields[0]->GetNpoints());
160 for (
int n = 0; n < nConvectiveFields; ++n)
162 fields[0]->PhysDeriv(inarray[n], grad0);
165 Array<OneD, NekDouble> Outarray(nPointsTot);
166 fields[0]->PhysInterp1DScaled(OneDptscale, grad0, wkSp);
167 Vmath::Vmul(nPointsTot, wkSp, 1, AdvVel[0], 1, Outarray, 1);
170 fields[0]->PhysGalerkinProjection1DScaled(
171 OneDptscale, Outarray, outarray[n]);
173 timer.AccumulateRegion(
"GalerinProject");
177 Vmath::Vmul(nPointsTot, grad0, 1, AdvVel[0], 1, outarray[n],
183 grad0 = Array<OneD, NekDouble>(fields[0]->GetNpoints());
184 grad1 = Array<OneD, NekDouble>(fields[0]->GetNpoints());
185 for (
int n = 0; n < nConvectiveFields; ++n)
187 fields[0]->PhysDeriv(inarray[n], grad0, grad1);
191 Array<OneD, NekDouble> Outarray(nPointsTot);
192 fields[0]->PhysInterp1DScaled(OneDptscale, grad0, wkSp);
193 Vmath::Vmul(nPointsTot, wkSp, 1, AdvVel[0], 1, Outarray, 1);
195 fields[0]->PhysInterp1DScaled(OneDptscale, grad1, wkSp);
197 timer.AccumulateRegion(
"Interp1DScaled");
198 Vmath::Vvtvp(nPointsTot, wkSp, 1, AdvVel[1], 1, Outarray, 1,
202 fields[0]->PhysGalerkinProjection1DScaled(
203 OneDptscale, Outarray, outarray[n]);
205 timer.AccumulateRegion(
"GalerinProject");
209 Vmath::Vmul(nPointsTot, grad0, 1, AdvVel[0], 1, outarray[n],
212 outarray[n], 1, outarray[n], 1);
219 Array<OneD, Array<OneD, NekDouble>> grad(ndim);
220 Array<OneD, Array<OneD, NekDouble>> gradScaled(
221 ndim * nConvectiveFields);
222 Array<OneD, Array<OneD, NekDouble>> Outarray(nConvectiveFields);
223 for (
int i = 0; i < ndim; i++)
225 grad[i] = Array<OneD, NekDouble>(fields[0]->GetNpoints());
227 for (
int i = 0; i < ndim * nConvectiveFields; i++)
229 gradScaled[i] = Array<OneD, NekDouble>(nPointsTot);
231 for (
int i = 0; i < nConvectiveFields; i++)
233 Outarray[i] = Array<OneD, NekDouble>(nPointsTot);
236 for (
int n = 0; n < nConvectiveFields; n++)
238 fields[0]->PhysDeriv(inarray[n], grad[0], grad[1], grad[2]);
239 for (
int i = 0; i < ndim; i++)
242 fields[0]->PhysInterp1DScaled(OneDptscale, grad[i],
243 gradScaled[n * ndim + i]);
245 timer.AccumulateRegion(
"Interp1DScaled");
249 fields[0]->DealiasedDotProd(AdvVel, gradScaled, Outarray);
252 for (
int n = 0; n < nConvectiveFields; n++)
254 fields[0]->PhysGalerkinProjection1DScaled(
255 OneDptscale, Outarray[n], outarray[n]);
258 timer.AccumulateRegion(
"GalerinProject");
263 Array<OneD, Array<OneD, NekDouble>> grad(ndim *
265 Array<OneD, Array<OneD, NekDouble>> Outarray(nConvectiveFields);
266 for (
int i = 0; i < ndim * nConvectiveFields; i++)
268 grad[i] = Array<OneD, NekDouble>(nPointsTot);
270 for (
int i = 0; i < nConvectiveFields; i++)
272 Outarray[i] = Array<OneD, NekDouble>(nPointsTot);
275 for (
int n = 0; n < nConvectiveFields; n++)
277 fields[0]->PhysDeriv(inarray[n], grad[n * ndim + 0],
282 fields[0]->DealiasedDotProd(AdvVel, grad, outarray);
286 grad0 = Array<OneD, NekDouble>(fields[0]->GetNpoints());
287 grad1 = Array<OneD, NekDouble>(fields[0]->GetNpoints());
288 grad2 = Array<OneD, NekDouble>(fields[0]->GetNpoints());
289 Array<OneD, NekDouble> tmp = grad2;
290 for (
int n = 0; n < nConvectiveFields; ++n)
292 if (fields[0]->GetWaveSpace() ==
true &&
295 fields[0]->HomogeneousBwdTrans(inarray[n], tmp);
296 fields[0]->PhysDeriv(tmp, grad0, grad1);
299 inarray[n], outarray[n]);
300 fields[0]->HomogeneousBwdTrans(outarray[n], grad2);
302 else if (fields[0]->GetWaveSpace() ==
true &&
305 fields[0]->HomogeneousBwdTrans(inarray[n], tmp);
306 fields[0]->PhysDeriv(tmp, grad0);
309 inarray[n], outarray[n]);
310 fields[0]->HomogeneousBwdTrans(outarray[n], grad1);
313 inarray[n], outarray[n]);
314 fields[0]->HomogeneousBwdTrans(outarray[n], grad2);
318 fields[0]->PhysDeriv(inarray[n], grad0, grad1, grad2);
323 Array<OneD, NekDouble> Outarray(nPointsTot);
325 fields[0]->PhysInterp1DScaled(OneDptscale, grad0, wkSp);
327 timer.AccumulateRegion(
"Interp1DScaled");
328 Vmath::Vmul(nPointsTot, wkSp, 1, AdvVel[0], 1, Outarray,
332 fields[0]->PhysInterp1DScaled(OneDptscale, grad1, wkSp);
334 timer.AccumulateRegion(
"Interp1DScaled");
336 Outarray, 1, Outarray, 1);
339 fields[0]->PhysInterp1DScaled(OneDptscale, grad2, wkSp);
341 timer.AccumulateRegion(
"Interp1DScaled");
343 Outarray, 1, Outarray, 1);
345 fields[0]->PhysGalerkinProjection1DScaled(
346 OneDptscale, Outarray, outarray[n]);
348 timer.AccumulateRegion(
"GalerinProject");
355 outarray[n], 1, outarray[n], 1);
357 outarray[n], 1, outarray[n], 1);
360 if (fields[0]->GetWaveSpace() ==
true)
362 fields[0]->HomogeneousFwdTrans(outarray[n],
369 ASSERTL0(
false,
"dimension unknown");
372 for (
int n = 0; n < nConvectiveFields; ++n)
#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