Advects a vector field.
95 int nqtot = fields[0]->GetTotPoints();
96 ASSERTL1(nConvectiveFields == inarray.num_elements(),
"Number of convective fields and Inarray are not compatible");
99 int ndim = advVel.num_elements();
100 Array<OneD, Array<OneD, NekDouble> > AdvVel (advVel.num_elements());
102 Array<OneD, Array<OneD, NekDouble> > velocity(ndim);
103 for(
int i = 0; i < ndim; ++i)
108 velocity[i] = Array<OneD, NekDouble>(nqtot,0.0);
109 fields[i]->HomogeneousBwdTrans(advVel[i],velocity[i]);
113 velocity[i] = advVel[i];
117 int nPointsTot = fields[0]->GetNpoints();
118 Array<OneD, NekDouble> grad0,grad1,grad2,wkSp;
125 nPointsTot = fields[0]->Get1DScaledTotPoints(OneDptscale);
131 for(
int i = 0; i < ndim; ++i)
133 AdvVel[i] = Array<OneD, NekDouble> (nPointsTot);
135 fields[0]->PhysInterp1DScaled(OneDptscale,velocity[i],AdvVel[i]);
140 for(
int i = 0; i < ndim; ++i)
142 AdvVel[i] = velocity[i];
146 wkSp = Array<OneD, NekDouble> (nPointsTot);
152 grad0 = Array<OneD, NekDouble> (fields[0]->GetNpoints());
153 for(
int n = 0; n < nConvectiveFields; ++n)
155 fields[0]->PhysDeriv(inarray[n],grad0);
158 Array<OneD, NekDouble> Outarray(nPointsTot);
159 fields[0]->PhysInterp1DScaled(OneDptscale,grad0,wkSp);
160 Vmath::Vmul (nPointsTot,wkSp,1,AdvVel[0],1,Outarray,1);
162 fields[0]->PhysGalerkinProjection1DScaled(OneDptscale,Outarray,outarray[n]);
166 Vmath::Vmul(nPointsTot,grad0,1,AdvVel[0],1,outarray[n],1);
171 grad0 = Array<OneD, NekDouble> (fields[0]->GetNpoints());
172 grad1 = Array<OneD, NekDouble> (fields[0]->GetNpoints());
173 for(
int n = 0; n < nConvectiveFields; ++n)
175 fields[0]->PhysDeriv(inarray[n],grad0,grad1);
179 Array<OneD, NekDouble> Outarray(nPointsTot);
180 fields[0]->PhysInterp1DScaled(OneDptscale,grad0,wkSp);
181 Vmath::Vmul (nPointsTot,wkSp,1,AdvVel[0],1,Outarray,1);
182 fields[0]->PhysInterp1DScaled(OneDptscale,grad1,wkSp);
183 Vmath::Vvtvp(nPointsTot,wkSp,1,AdvVel[1],1,Outarray,1,Outarray,1);
185 fields[0]->PhysGalerkinProjection1DScaled(OneDptscale,Outarray,outarray[n]);
189 Vmath::Vmul (nPointsTot,grad0,1,AdvVel[0],1,outarray[n],1);
190 Vmath::Vvtvp(nPointsTot,grad1,1,AdvVel[1],1,outarray[n],1,outarray[n],1);
197 Array<OneD, Array<OneD, NekDouble> > grad (ndim);
198 Array<OneD, Array<OneD, NekDouble> > gradScaled (ndim*nConvectiveFields);
199 Array<OneD, Array<OneD, NekDouble> > Outarray (nConvectiveFields);
200 for (
int i = 0; i < ndim; i++)
202 grad[i] = Array<OneD, NekDouble>(fields[0]->GetNpoints());
204 for (
int i = 0; i < ndim*nConvectiveFields; i++)
206 gradScaled[i] = Array<OneD, NekDouble>(nPointsTot);
208 for (
int i = 0; i < nConvectiveFields; i++)
210 Outarray[i] = Array<OneD, NekDouble>(nPointsTot);
213 for (
int n = 0; n < nConvectiveFields; n++)
215 fields[0]->PhysDeriv(inarray[n],grad[0],grad[1],grad[2]);
216 for (
int i = 0; i < ndim; i++)
218 fields[0]->PhysInterp1DScaled(OneDptscale,grad[i],
219 gradScaled[n*ndim+i]);
223 fields[0]->DealiasedDotProd(AdvVel,gradScaled,Outarray,
m_CoeffState);
225 for (
int n = 0; n < nConvectiveFields; n++)
227 fields[0]->PhysGalerkinProjection1DScaled(OneDptscale,
228 Outarray[n],outarray[n]);
233 Array<OneD, Array<OneD, NekDouble> > grad (ndim*nConvectiveFields);
234 Array<OneD, Array<OneD, NekDouble> > Outarray (nConvectiveFields);
235 for (
int i = 0; i < ndim*nConvectiveFields; i++)
237 grad[i] = Array<OneD, NekDouble>(nPointsTot);
239 for (
int i = 0; i < nConvectiveFields; i++)
241 Outarray[i] = Array<OneD, NekDouble>(nPointsTot);
244 for (
int n = 0; n < nConvectiveFields; n++)
246 fields[0]->PhysDeriv(inarray[n],grad[n*ndim+0],
251 fields[0]->DealiasedDotProd(AdvVel,grad,outarray,
m_CoeffState);
255 grad0 = Array<OneD, NekDouble> (fields[0]->GetNpoints());
256 grad1 = Array<OneD, NekDouble> (fields[0]->GetNpoints());
257 grad2 = Array<OneD, NekDouble> (fields[0]->GetNpoints());
258 for(
int n = 0; n < nConvectiveFields; ++n)
260 if (fields[0]->GetWaveSpace() ==
true &&
266 fields[0]->PhysDeriv(velocity[n],grad0,grad1);
270 fields[0]->HomogeneousBwdTrans(inarray[n],wkSp);
271 fields[0]->PhysDeriv(wkSp,grad0,grad1);
277 fields[0]->HomogeneousBwdTrans(outarray[n],grad2);
279 else if (fields[0]->GetWaveSpace() ==
true &&
285 fields[0]->PhysDeriv(velocity[n],grad0);
289 fields[0]->HomogeneousBwdTrans(inarray[n],wkSp);
290 fields[0]->PhysDeriv(wkSp,grad0);
295 fields[0]->HomogeneousBwdTrans(outarray[n],grad1);
299 fields[0]->HomogeneousBwdTrans(outarray[n],grad2);
303 fields[0]->PhysDeriv(inarray[n],grad0,grad1,grad2);
307 Array<OneD, NekDouble> Outarray(nPointsTot);
308 fields[0]->PhysInterp1DScaled(OneDptscale,grad0,wkSp);
309 Vmath::Vmul(nPointsTot,wkSp,1,AdvVel[0],1,Outarray,1);
311 fields[0]->PhysInterp1DScaled(OneDptscale,grad1,wkSp);
315 fields[0]->PhysInterp1DScaled(OneDptscale,grad2,wkSp);
318 fields[0]->PhysGalerkinProjection1DScaled(OneDptscale,
319 Outarray,outarray[n]);
323 Vmath::Vmul(nPointsTot,grad0,1,AdvVel[0],1,outarray[n],1);
324 Vmath::Vvtvp(nPointsTot,grad1,1,AdvVel[1],1,outarray[n],1,
326 Vmath::Vvtvp(nPointsTot,grad2,1,AdvVel[2],1,outarray[n],1,
330 if(fields[0]->GetWaveSpace() ==
true)
332 fields[0]->HomogeneousFwdTrans(outarray[n],outarray[n]);
338 ASSERTL0(
false,
"dimension unknown");
341 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....
MultiRegions::CoeffState m_CoeffState
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