Advects a vector field.
95 boost::ignore_unused(time, pFwd, pBwd);
97 int nqtot = fields[0]->GetTotPoints();
98 ASSERTL1(nConvectiveFields == inarray.size(),
99 "Number of convective fields and Inarray are not compatible");
102 int ndim = advVel.size();
103 Array<OneD, Array<OneD, NekDouble>> AdvVel(advVel.size());
105 Array<OneD, Array<OneD, NekDouble>> velocity(ndim);
107 LibUtilities::Timer timer;
109 for (
int i = 0; i < ndim; ++i)
114 velocity[i] = Array<OneD, NekDouble>(nqtot, 0.0);
115 fields[i]->HomogeneousBwdTrans(nqtot, advVel[i], velocity[i]);
119 velocity[i] = advVel[i];
123 int nPointsTot = fields[0]->GetNpoints();
124 Array<OneD, NekDouble> grad0, grad1, grad2, wkSp;
131 nPointsTot = fields[0]->Get1DScaledTotPoints(OneDptscale);
137 for (
int i = 0; i < ndim; ++i)
139 AdvVel[i] = Array<OneD, NekDouble>(nPointsTot);
142 fields[0]->PhysInterp1DScaled(OneDptscale, velocity[i], AdvVel[i]);
144 timer.AccumulateRegion(
"Interp1DScaled");
149 for (
int i = 0; i < ndim; ++i)
151 AdvVel[i] = velocity[i];
155 wkSp = Array<OneD, NekDouble>(nPointsTot);
161 grad0 = Array<OneD, NekDouble>(fields[0]->GetNpoints());
162 for (
int n = 0; n < nConvectiveFields; ++n)
164 fields[0]->PhysDeriv(inarray[n], grad0);
167 Array<OneD, NekDouble> Outarray(nPointsTot);
168 fields[0]->PhysInterp1DScaled(OneDptscale, grad0, wkSp);
169 Vmath::Vmul(nPointsTot, wkSp, 1, AdvVel[0], 1, Outarray, 1);
172 fields[0]->PhysGalerkinProjection1DScaled(
173 OneDptscale, Outarray, outarray[n]);
175 timer.AccumulateRegion(
"GalerinProject");
179 Vmath::Vmul(nPointsTot, grad0, 1, AdvVel[0], 1, outarray[n],
185 grad0 = Array<OneD, NekDouble>(nqtot);
186 grad1 = Array<OneD, NekDouble>(nqtot);
187 for (
int n = 0; n < nConvectiveFields; ++n)
189 fields[0]->PhysDeriv(inarray[n], grad0, grad1);
193 Array<OneD, NekDouble> Outarray(nPointsTot);
194 fields[0]->PhysInterp1DScaled(OneDptscale, grad0, wkSp);
195 Vmath::Vmul(nPointsTot, wkSp, 1, AdvVel[0], 1, Outarray, 1);
197 fields[0]->PhysInterp1DScaled(OneDptscale, grad1, wkSp);
199 timer.AccumulateRegion(
"Interp1DScaled");
200 Vmath::Vvtvp(nPointsTot, wkSp, 1, AdvVel[1], 1, Outarray, 1,
204 fields[0]->PhysGalerkinProjection1DScaled(
205 OneDptscale, Outarray, outarray[n]);
207 timer.AccumulateRegion(
"GalerinProject");
211 Vmath::Vmul(nPointsTot, grad0, 1, AdvVel[0], 1, outarray[n],
214 outarray[n], 1, outarray[n], 1);
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++)
227 grad[i] = Array<OneD, NekDouble>(nqtot);
229 for (
int i = 0; i < ndim * nConvectiveFields; i++)
231 gradScaled[i] = Array<OneD, NekDouble>(nPointsTot);
233 for (
int i = 0; i < nConvectiveFields; i++)
235 Outarray[i] = Array<OneD, NekDouble>(nPointsTot);
238 for (
int n = 0; n < nConvectiveFields; n++)
240 fields[0]->PhysDeriv(inarray[n], grad[0], grad[1], grad[2]);
241 for (
int i = 0; i < ndim; i++)
244 fields[0]->PhysInterp1DScaled(OneDptscale, grad[i],
245 gradScaled[n * ndim + i]);
247 timer.AccumulateRegion(
"Interp1DScaled");
251 fields[0]->DealiasedDotProd(nPointsTot, AdvVel, gradScaled,
255 for (
int n = 0; n < nConvectiveFields; n++)
257 fields[0]->PhysGalerkinProjection1DScaled(
258 OneDptscale, Outarray[n], outarray[n]);
261 timer.AccumulateRegion(
"GalerinProject");
266 Array<OneD, Array<OneD, NekDouble>> grad(ndim *
268 Array<OneD, Array<OneD, NekDouble>> Outarray(nConvectiveFields);
269 for (
int i = 0; i < ndim * nConvectiveFields; i++)
271 grad[i] = Array<OneD, NekDouble>(nPointsTot);
273 for (
int i = 0; i < nConvectiveFields; i++)
275 Outarray[i] = Array<OneD, NekDouble>(nPointsTot);
278 for (
int n = 0; n < nConvectiveFields; n++)
280 fields[0]->PhysDeriv(inarray[n], grad[n * ndim + 0],
285 fields[0]->DealiasedDotProd(nPointsTot, AdvVel, grad, outarray);
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)
295 if (fields[0]->GetWaveSpace() ==
true &&
298 fields[0]->HomogeneousBwdTrans(nqtot, inarray[n], tmp);
299 fields[0]->PhysDeriv(tmp, grad0, grad1);
302 inarray[n], outarray[n]);
303 fields[0]->HomogeneousBwdTrans(nqtot, outarray[n],
306 else if (fields[0]->GetWaveSpace() ==
true &&
309 fields[0]->HomogeneousBwdTrans(nqtot, inarray[n], tmp);
310 fields[0]->PhysDeriv(tmp, grad0);
313 inarray[n], outarray[n]);
314 fields[0]->HomogeneousBwdTrans(nqtot, outarray[n],
318 inarray[n], outarray[n]);
319 fields[0]->HomogeneousBwdTrans(nqtot, outarray[n],
324 fields[0]->PhysDeriv(inarray[n], grad0, grad1, grad2);
329 Array<OneD, NekDouble> Outarray(nPointsTot);
331 fields[0]->PhysInterp1DScaled(OneDptscale, grad0, wkSp);
333 timer.AccumulateRegion(
"Interp1DScaled");
334 Vmath::Vmul(nPointsTot, wkSp, 1, AdvVel[0], 1, Outarray,
338 fields[0]->PhysInterp1DScaled(OneDptscale, grad1, wkSp);
340 timer.AccumulateRegion(
"Interp1DScaled");
342 Outarray, 1, Outarray, 1);
345 fields[0]->PhysInterp1DScaled(OneDptscale, grad2, wkSp);
347 timer.AccumulateRegion(
"Interp1DScaled");
349 Outarray, 1, Outarray, 1);
351 fields[0]->PhysGalerkinProjection1DScaled(
352 OneDptscale, Outarray, outarray[n]);
354 timer.AccumulateRegion(
"GalerinProject");
361 outarray[n], 1, outarray[n], 1);
363 outarray[n], 1, outarray[n], 1);
366 if (fields[0]->GetWaveSpace() ==
true)
368 fields[0]->HomogeneousFwdTrans(nqtot, outarray[n],
375 ASSERTL0(
false,
"dimension unknown");
378 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