Advects a vector field.
87 int ndim = advVel.size();
88 int nqtot = fields[0]->GetTotPoints();
89 ASSERTL1(nConvectiveFields == inarray.size(),
90 "Number of convective fields and Inarray are not compatible");
92 Array<OneD, Array<OneD, NekDouble>> velocity(ndim);
93 for (
int i = 0; i < ndim; ++i)
97 velocity[i] = Array<OneD, NekDouble>(nqtot, 0.0);
98 fields[i]->HomogeneousBwdTrans(advVel[i], velocity[i]);
102 velocity[i] = advVel[i];
106 for (
int n = 0; n < nConvectiveFields; ++n)
110 int nPointsTot = fields[0]->GetNpoints();
111 Array<OneD, NekDouble> gradV0, gradV1, gradV2, tmp, Up;
113 gradV0 = Array<OneD, NekDouble>(nPointsTot);
114 tmp = Array<OneD, NekDouble>(nPointsTot);
120 fields[0]->PhysDeriv(inarray[n], gradV0);
121 Vmath::Vmul(nPointsTot, gradV0, 1, velocity[0], 1, outarray[n],
123 Vmath::Vmul(nPointsTot, inarray[n], 1, velocity[0], 1, gradV0,
125 fields[0]->PhysDeriv(gradV0, tmp);
126 Vmath::Vadd(nPointsTot, tmp, 1, outarray[n], 1, outarray[n], 1);
127 Vmath::Smul(nPointsTot, 0.5, outarray[n], 1, outarray[n], 1);
130 gradV1 = Array<OneD, NekDouble>(nPointsTot);
131 fields[0]->PhysDeriv(inarray[n], gradV0, gradV1);
132 Vmath::Vmul(nPointsTot, gradV0, 1, velocity[0], 1, outarray[n],
134 Vmath::Vvtvp(nPointsTot, gradV1, 1, velocity[1], 1, outarray[n],
136 Vmath::Vmul(nPointsTot, inarray[n], 1, velocity[0], 1, gradV0,
138 Vmath::Vmul(nPointsTot, inarray[n], 1, velocity[1], 1, gradV1,
142 Vmath::Vadd(nPointsTot, tmp, 1, outarray[n], 1, outarray[n], 1);
145 Vmath::Vadd(nPointsTot, tmp, 1, outarray[n], 1, outarray[n], 1);
146 Vmath::Smul(nPointsTot, 0.5, outarray[n], 1, outarray[n], 1);
149 gradV1 = Array<OneD, NekDouble>(nPointsTot);
150 gradV2 = Array<OneD, NekDouble>(nPointsTot);
152 fields[0]->PhysDeriv(inarray[n], gradV0, gradV1, gradV2);
158 fields[0]->GetWaveSpace() ==
false)
160 fields[0]->DealiasedProd(velocity[0], gradV0, gradV0);
161 fields[0]->DealiasedProd(velocity[1], gradV1, gradV1);
162 fields[0]->DealiasedProd(velocity[2], gradV2, gradV2);
163 Vmath::Vadd(nPointsTot, gradV0, 1, gradV1, 1, outarray[n],
167 fields[0]->DealiasedProd(inarray[n], velocity[0], gradV0);
168 fields[0]->DealiasedProd(inarray[n], velocity[1], gradV1);
169 fields[0]->DealiasedProd(inarray[n], velocity[2], gradV2);
172 Vmath::Vadd(nPointsTot, tmp, 1, outarray[n], 1, outarray[n],
176 Vmath::Vadd(nPointsTot, tmp, 1, outarray[n], 1, outarray[n],
180 Vmath::Vadd(nPointsTot, tmp, 1, outarray[n], 1, outarray[n],
182 Vmath::Smul(nPointsTot, 0.5, outarray[n], 1, outarray[n],
185 else if (fields[0]->GetWaveSpace() ==
true &&
188 Up = Array<OneD, NekDouble>(nPointsTot);
191 fields[0]->HomogeneousBwdTrans(gradV0, tmp);
192 Vmath::Vmul(nPointsTot, tmp, 1, velocity[0], 1, outarray[n],
194 fields[0]->HomogeneousBwdTrans(gradV1, tmp);
196 outarray[n], 1, outarray[n], 1);
197 fields[0]->HomogeneousBwdTrans(gradV2, tmp);
199 outarray[n], 1, outarray[n], 1);
201 fields[0]->HomogeneousBwdTrans(inarray[n], Up);
202 Vmath::Vmul(nPointsTot, Up, 1, velocity[0], 1, gradV0, 1);
203 Vmath::Vmul(nPointsTot, Up, 1, velocity[1], 1, gradV1, 1);
204 Vmath::Vmul(nPointsTot, Up, 1, velocity[2], 1, gradV2, 1);
206 fields[0]->SetWaveSpace(
false);
209 Vmath::Vadd(nPointsTot, tmp, 1, outarray[n], 1, outarray[n],
213 Vmath::Vadd(nPointsTot, tmp, 1, outarray[n], 1, outarray[n],
217 Vmath::Vadd(nPointsTot, tmp, 1, outarray[n], 1, outarray[n],
219 fields[0]->SetWaveSpace(
true);
221 Vmath::Smul(nPointsTot, 0.5, outarray[n], 1, tmp, 1);
222 fields[0]->HomogeneousFwdTrans(tmp, outarray[n]);
224 else if (fields[0]->GetWaveSpace() ==
false &&
230 outarray[n], 1, outarray[n], 1);
232 outarray[n], 1, outarray[n], 1);
233 Vmath::Vmul(nPointsTot, inarray[n], 1, velocity[0], 1,
235 Vmath::Vmul(nPointsTot, inarray[n], 1, velocity[1], 1,
237 Vmath::Vmul(nPointsTot, inarray[n], 1, velocity[2], 1,
241 Vmath::Vadd(nPointsTot, tmp, 1, outarray[n], 1, outarray[n],
245 Vmath::Vadd(nPointsTot, tmp, 1, outarray[n], 1, outarray[n],
249 Vmath::Vadd(nPointsTot, tmp, 1, outarray[n], 1, outarray[n],
251 Vmath::Smul(nPointsTot, 0.5, outarray[n], 1, outarray[n],
257 "Dealiasing is not allowed in combination "
258 "with the Skew-Symmetric advection form for "
259 "efficiency reasons.");
263 ASSERTL0(
false,
"dimension unknown");
#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
void Vadd(int n, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
Add vector z = x+y.
void Smul(int n, const T alpha, const T *x, const int incx, T *y, const int incy)
Scalar multiply y = alpha*x.