Advects a vector field.
86 boost::ignore_unused(time, pFwd, pBwd);
89 int ndim = advVel.size();
90 int nqtot = fields[0]->GetTotPoints();
91 ASSERTL1(nConvectiveFields == inarray.size(),
92 "Number of convective fields and Inarray are not compatible");
94 Array<OneD, Array<OneD, NekDouble>> velocity(ndim);
95 for (
int i = 0; i < ndim; ++i)
99 velocity[i] = Array<OneD, NekDouble>(nqtot, 0.0);
100 fields[i]->HomogeneousBwdTrans(nqtot, advVel[i], velocity[i]);
104 velocity[i] = advVel[i];
108 for (
int n = 0; n < nConvectiveFields; ++n)
112 int nPointsTot = fields[0]->GetNpoints();
113 Array<OneD, NekDouble> gradV0, gradV1, gradV2, tmp, Up;
115 gradV0 = Array<OneD, NekDouble>(nPointsTot);
116 tmp = Array<OneD, NekDouble>(nPointsTot);
122 fields[0]->PhysDeriv(inarray[n], gradV0);
123 Vmath::Vmul(nPointsTot, gradV0, 1, velocity[0], 1, outarray[n],
125 Vmath::Vmul(nPointsTot, inarray[n], 1, velocity[0], 1, gradV0,
127 fields[0]->PhysDeriv(gradV0, tmp);
128 Vmath::Vadd(nPointsTot, tmp, 1, outarray[n], 1, outarray[n], 1);
129 Vmath::Smul(nPointsTot, 0.5, outarray[n], 1, outarray[n], 1);
132 gradV1 = Array<OneD, NekDouble>(nPointsTot);
133 fields[0]->PhysDeriv(inarray[n], gradV0, gradV1);
134 Vmath::Vmul(nPointsTot, gradV0, 1, velocity[0], 1, outarray[n],
136 Vmath::Vvtvp(nPointsTot, gradV1, 1, velocity[1], 1, outarray[n],
138 Vmath::Vmul(nPointsTot, inarray[n], 1, velocity[0], 1, gradV0,
140 Vmath::Vmul(nPointsTot, inarray[n], 1, velocity[1], 1, gradV1,
144 Vmath::Vadd(nPointsTot, tmp, 1, outarray[n], 1, outarray[n], 1);
147 Vmath::Vadd(nPointsTot, tmp, 1, outarray[n], 1, outarray[n], 1);
148 Vmath::Smul(nPointsTot, 0.5, outarray[n], 1, outarray[n], 1);
151 gradV1 = Array<OneD, NekDouble>(nPointsTot);
152 gradV2 = Array<OneD, NekDouble>(nPointsTot);
154 fields[0]->PhysDeriv(inarray[n], gradV0, gradV1, gradV2);
160 fields[0]->GetWaveSpace() ==
false)
162 fields[0]->DealiasedProd(nPointsTot, velocity[0], gradV0,
164 fields[0]->DealiasedProd(nPointsTot, velocity[1], gradV1,
166 fields[0]->DealiasedProd(nPointsTot, velocity[2], gradV2,
168 Vmath::Vadd(nPointsTot, gradV0, 1, gradV1, 1, outarray[n],
172 fields[0]->DealiasedProd(nPointsTot, inarray[n],
173 velocity[0], gradV0);
174 fields[0]->DealiasedProd(nPointsTot, inarray[n],
175 velocity[1], gradV1);
176 fields[0]->DealiasedProd(nPointsTot, inarray[n],
177 velocity[2], gradV2);
180 Vmath::Vadd(nPointsTot, tmp, 1, outarray[n], 1, outarray[n],
184 Vmath::Vadd(nPointsTot, tmp, 1, outarray[n], 1, outarray[n],
188 Vmath::Vadd(nPointsTot, tmp, 1, outarray[n], 1, outarray[n],
190 Vmath::Smul(nPointsTot, 0.5, outarray[n], 1, outarray[n],
193 else if (fields[0]->GetWaveSpace() ==
true &&
196 Up = Array<OneD, NekDouble>(nPointsTot);
199 fields[0]->HomogeneousBwdTrans(nPointsTot, gradV0, tmp);
200 Vmath::Vmul(nPointsTot, tmp, 1, velocity[0], 1, outarray[n],
202 fields[0]->HomogeneousBwdTrans(nPointsTot, gradV1, tmp);
204 outarray[n], 1, outarray[n], 1);
205 fields[0]->HomogeneousBwdTrans(nPointsTot, gradV2, tmp);
207 outarray[n], 1, outarray[n], 1);
209 fields[0]->HomogeneousBwdTrans(nPointsTot, inarray[n], Up);
210 Vmath::Vmul(nPointsTot, Up, 1, velocity[0], 1, gradV0, 1);
211 Vmath::Vmul(nPointsTot, Up, 1, velocity[1], 1, gradV1, 1);
212 Vmath::Vmul(nPointsTot, Up, 1, velocity[2], 1, gradV2, 1);
214 fields[0]->SetWaveSpace(
false);
217 Vmath::Vadd(nPointsTot, tmp, 1, outarray[n], 1, outarray[n],
221 Vmath::Vadd(nPointsTot, tmp, 1, outarray[n], 1, outarray[n],
225 Vmath::Vadd(nPointsTot, tmp, 1, outarray[n], 1, outarray[n],
227 fields[0]->SetWaveSpace(
true);
229 Vmath::Smul(nPointsTot, 0.5, outarray[n], 1, tmp, 1);
230 fields[0]->HomogeneousFwdTrans(nPointsTot, tmp,
233 else if (fields[0]->GetWaveSpace() ==
false &&
239 outarray[n], 1, outarray[n], 1);
241 outarray[n], 1, outarray[n], 1);
242 Vmath::Vmul(nPointsTot, inarray[n], 1, velocity[0], 1,
244 Vmath::Vmul(nPointsTot, inarray[n], 1, velocity[1], 1,
246 Vmath::Vmul(nPointsTot, inarray[n], 1, velocity[2], 1,
250 Vmath::Vadd(nPointsTot, tmp, 1, outarray[n], 1, outarray[n],
254 Vmath::Vadd(nPointsTot, tmp, 1, outarray[n], 1, outarray[n],
258 Vmath::Vadd(nPointsTot, tmp, 1, outarray[n], 1, outarray[n],
260 Vmath::Smul(nPointsTot, 0.5, outarray[n], 1, outarray[n],
266 "Dealiasing is not allowed in combination "
267 "with the Skew-Symmetric advection form for "
268 "efficiency reasons.");
272 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.