Advects a vector field.
92 int nqtot = fields[0]->GetTotPoints();
93 ASSERTL1(nConvectiveFields == inarray.num_elements(),
"Number of convective fields and Inarray are not compatible");
95 Array<OneD, NekDouble > Deriv = Array<OneD, NekDouble> (nqtot*nConvectiveFields);
97 for(
int n = 0; n < nConvectiveFields; ++n)
100 int ndim = advVel.num_elements();
104 int nPointsTot = fields[0]->GetNpoints();
105 Array<OneD, NekDouble> gradV0,gradV1,gradV2, tmp, Up;
107 gradV0 = Array<OneD, NekDouble> (nPointsTot);
108 tmp = Array<OneD, NekDouble> (nPointsTot);
114 fields[0]->PhysDeriv(inarray[n],gradV0);
115 Vmath::Vmul(nPointsTot,gradV0,1,advVel[0],1,outarray[n],1);
116 Vmath::Vmul(nPointsTot,inarray[n],1,advVel[0],1,gradV0,1);
117 fields[0]->PhysDeriv(gradV0,tmp);
118 Vmath::Vadd(nPointsTot,tmp,1,outarray[n],1,outarray[n],1);
119 Vmath::Smul(nPointsTot,0.5,outarray[n],1,outarray[n],1);
122 gradV1 = Array<OneD, NekDouble> (nPointsTot);
123 fields[0]->PhysDeriv(inarray[n],gradV0,gradV1);
124 Vmath::Vmul (nPointsTot,gradV0,1,advVel[0],1,outarray[n],1);
125 Vmath::Vvtvp(nPointsTot,gradV1,1,advVel[1],1,outarray[n],1,outarray[n],1);
126 Vmath::Vmul(nPointsTot,inarray[n],1,advVel[0],1,gradV0,1);
127 Vmath::Vmul(nPointsTot,inarray[n],1,advVel[1],1,gradV1,1);
129 Vmath::Vadd(nPointsTot,tmp,1,outarray[n],1,outarray[n],1);
131 Vmath::Vadd(nPointsTot,tmp,1,outarray[n],1,outarray[n],1);
132 Vmath::Smul(nPointsTot,0.5,outarray[n],1,outarray[n],1);
135 gradV1 = Array<OneD, NekDouble> (nPointsTot);
136 gradV2 = Array<OneD, NekDouble> (nPointsTot);
138 fields[0]->PhysDeriv(inarray[n],gradV0,gradV1,gradV2);
144 fields[0]->DealiasedProd(advVel[0],gradV0,gradV0,
m_CoeffState);
145 fields[0]->DealiasedProd(advVel[1],gradV1,gradV1,
m_CoeffState);
146 fields[0]->DealiasedProd(advVel[2],gradV2,gradV2,
m_CoeffState);
147 Vmath::Vadd(nPointsTot,gradV0,1,gradV1,1,outarray[n],1);
148 Vmath::Vadd(nPointsTot,gradV2,1,outarray[n],1,outarray[n],1);
149 fields[0]->DealiasedProd(inarray[n],advVel[0],gradV0,
m_CoeffState);
150 fields[0]->DealiasedProd(inarray[n],advVel[1],gradV1,
m_CoeffState);
151 fields[0]->DealiasedProd(inarray[n],advVel[2],gradV2,
m_CoeffState);
153 Vmath::Vadd(nPointsTot,tmp,1,outarray[n],1,outarray[n],1);
155 Vmath::Vadd(nPointsTot,tmp,1,outarray[n],1,outarray[n],1);
157 Vmath::Vadd(nPointsTot,tmp,1,outarray[n],1,outarray[n],1);
158 Vmath::Smul(nPointsTot,0.5,outarray[n],1,outarray[n],1);
162 Up = Array<OneD, NekDouble> (nPointsTot);
165 fields[0]->HomogeneousBwdTrans(gradV0,tmp);
166 Vmath::Vmul(nPointsTot,tmp,1,advVel[0],1,outarray[n],1);
167 fields[0]->HomogeneousBwdTrans(gradV1,tmp);
168 Vmath::Vvtvp(nPointsTot,tmp,1,advVel[1],1,outarray[n],1,outarray[n],1);
169 fields[0]->HomogeneousBwdTrans(gradV2,tmp);
170 Vmath::Vvtvp(nPointsTot,tmp,1,advVel[2],1,outarray[n],1,outarray[n],1);
172 fields[0]->HomogeneousBwdTrans(inarray[n],Up);
177 fields[0]->SetWaveSpace(
false);
179 Vmath::Vadd(nPointsTot,tmp,1,outarray[n],1,outarray[n],1);
181 Vmath::Vadd(nPointsTot,tmp,1,outarray[n],1,outarray[n],1);
183 Vmath::Vadd(nPointsTot,tmp,1,outarray[n],1,outarray[n],1);
184 fields[0]->SetWaveSpace(
true);
187 fields[0]->HomogeneousFwdTrans(tmp,outarray[n]);
191 Vmath::Vmul(nPointsTot,gradV0,1,advVel[0],1,outarray[n],1);
192 Vmath::Vvtvp(nPointsTot,gradV1,1,advVel[1],1,outarray[n],1,outarray[n],1);
193 Vmath::Vvtvp(nPointsTot,gradV2,1,advVel[2],1,outarray[n],1,outarray[n],1);
194 Vmath::Vmul(nPointsTot,inarray[n],1,advVel[0],1,gradV0,1);
195 Vmath::Vmul(nPointsTot,inarray[n],1,advVel[1],1,gradV1,1);
196 Vmath::Vmul(nPointsTot,inarray[n],1,advVel[2],1,gradV2,1);
198 Vmath::Vadd(nPointsTot,tmp,1,outarray[n],1,outarray[n],1);
200 Vmath::Vadd(nPointsTot,tmp,1,outarray[n],1,outarray[n],1);
202 Vmath::Vadd(nPointsTot,tmp,1,outarray[n],1,outarray[n],1);
203 Vmath::Smul(nPointsTot,0.5,outarray[n],1,outarray[n],1);
207 ASSERTL0(
false,
"Dealiasing is not allowed in combination "
208 "with the Skew-Symmetric advection form for "
209 "efficiency reasons.");
213 ASSERTL0(
false,
"dimension unknown");
#define ASSERTL0(condition, msg)
MultiRegions::CoeffState m_CoeffState
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 Smul(int n, const T alpha, const T *x, const int incx, T *y, const int incy)
Scalar multiply y = alpha*y.
void Neg(int n, T *x, const int incx)
Negate x = -x.
MultiRegions::Direction const DirCartesianMap[]
bool m_homogen_dealiasing
#define ASSERTL1(condition, msg)
Assert Level 1 – Debugging which is used whether in FULLDEBUG or DEBUG compilation mode...
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 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.