Advects a vector field.
77{
78
79 int ndim = advVel.size();
80 int nqtot = fields[0]->GetTotPoints();
81 ASSERTL1(nConvectiveFields == inarray.size(),
82 "Number of convective fields and Inarray are not compatible");
83
84 Array<OneD, Array<OneD, NekDouble>> velocity(ndim);
85 for (int i = 0; i < ndim; ++i)
86 {
88 {
89 velocity[i] = Array<OneD, NekDouble>(nqtot, 0.0);
90 fields[i]->HomogeneousBwdTrans(nqtot, advVel[i], velocity[i]);
91 }
92 else
93 {
94 velocity[i] = advVel[i];
95 }
96 }
97
98 for (int n = 0; n < nConvectiveFields; ++n)
99 {
100
101
102 int nPointsTot = fields[0]->GetNpoints();
103 Array<OneD, NekDouble> gradV0, gradV1, gradV2, tmp, Up;
104
105 gradV0 = Array<OneD, NekDouble>(nPointsTot);
106 tmp = Array<OneD, NekDouble>(nPointsTot);
107
108
109 switch (ndim)
110 {
111 case 1:
112 fields[0]->PhysDeriv(inarray[n], gradV0);
113 Vmath::Vmul(nPointsTot, gradV0, 1, velocity[0], 1, outarray[n],
114 1);
115 Vmath::Vmul(nPointsTot, inarray[n], 1, velocity[0], 1, gradV0,
116 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);
120 break;
121 case 2:
122 gradV1 = Array<OneD, NekDouble>(nPointsTot);
123 fields[0]->PhysDeriv(inarray[n], gradV0, gradV1);
124 Vmath::Vmul(nPointsTot, gradV0, 1, velocity[0], 1, outarray[n],
125 1);
126 Vmath::Vvtvp(nPointsTot, gradV1, 1, velocity[1], 1, outarray[n],
127 1, outarray[n], 1);
128 Vmath::Vmul(nPointsTot, inarray[n], 1, velocity[0], 1, gradV0,
129 1);
130 Vmath::Vmul(nPointsTot, inarray[n], 1, velocity[1], 1, gradV1,
131 1);
133 tmp);
134 Vmath::Vadd(nPointsTot, tmp, 1, outarray[n], 1, outarray[n], 1);
136 tmp);
137 Vmath::Vadd(nPointsTot, tmp, 1, outarray[n], 1, outarray[n], 1);
138 Vmath::Smul(nPointsTot, 0.5, outarray[n], 1, outarray[n], 1);
139 break;
140 case 3:
141 gradV1 = Array<OneD, NekDouble>(nPointsTot);
142 gradV2 = Array<OneD, NekDouble>(nPointsTot);
143
144 fields[0]->PhysDeriv(inarray[n], gradV0, gradV1, gradV2);
145
146
147
148
150 fields[0]->GetWaveSpace() == false)
151 {
152 fields[0]->DealiasedProd(nPointsTot, velocity[0], gradV0,
153 gradV0);
154 fields[0]->DealiasedProd(nPointsTot, velocity[1], gradV1,
155 gradV1);
156 fields[0]->DealiasedProd(nPointsTot, velocity[2], gradV2,
157 gradV2);
158 Vmath::Vadd(nPointsTot, gradV0, 1, gradV1, 1, outarray[n],
159 1);
161 outarray[n], 1);
162 fields[0]->DealiasedProd(nPointsTot, inarray[n],
163 velocity[0], gradV0);
164 fields[0]->DealiasedProd(nPointsTot, inarray[n],
165 velocity[1], gradV1);
166 fields[0]->DealiasedProd(nPointsTot, inarray[n],
167 velocity[2], gradV2);
169 gradV0, tmp);
170 Vmath::Vadd(nPointsTot, tmp, 1, outarray[n], 1, outarray[n],
171 1);
173 gradV1, tmp);
174 Vmath::Vadd(nPointsTot, tmp, 1, outarray[n], 1, outarray[n],
175 1);
177 gradV2, tmp);
178 Vmath::Vadd(nPointsTot, tmp, 1, outarray[n], 1, outarray[n],
179 1);
180 Vmath::Smul(nPointsTot, 0.5, outarray[n], 1, outarray[n],
181 1);
182 }
183 else if (fields[0]->GetWaveSpace() == true &&
185 {
186 Up = Array<OneD, NekDouble>(nPointsTot);
187
188
189 fields[0]->HomogeneousBwdTrans(nPointsTot, gradV0, tmp);
190 Vmath::Vmul(nPointsTot, tmp, 1, velocity[0], 1, outarray[n],
191 1);
192 fields[0]->HomogeneousBwdTrans(nPointsTot, gradV1, tmp);
194 outarray[n], 1, outarray[n], 1);
195 fields[0]->HomogeneousBwdTrans(nPointsTot, gradV2, tmp);
197 outarray[n], 1, outarray[n], 1);
198
199 fields[0]->HomogeneousBwdTrans(nPointsTot, inarray[n], Up);
200 Vmath::Vmul(nPointsTot, Up, 1, velocity[0], 1, gradV0, 1);
201 Vmath::Vmul(nPointsTot, Up, 1, velocity[1], 1, gradV1, 1);
202 Vmath::Vmul(nPointsTot, Up, 1, velocity[2], 1, gradV2, 1);
203
204 fields[0]->SetWaveSpace(false);
206 gradV0, tmp);
207 Vmath::Vadd(nPointsTot, tmp, 1, outarray[n], 1, outarray[n],
208 1);
210 gradV1, tmp);
211 Vmath::Vadd(nPointsTot, tmp, 1, outarray[n], 1, outarray[n],
212 1);
214 gradV2, tmp);
215 Vmath::Vadd(nPointsTot, tmp, 1, outarray[n], 1, outarray[n],
216 1);
217 fields[0]->SetWaveSpace(true);
218
219 Vmath::Smul(nPointsTot, 0.5, outarray[n], 1, tmp, 1);
220 fields[0]->HomogeneousFwdTrans(nPointsTot, tmp,
221 outarray[n]);
222 }
223 else if (fields[0]->GetWaveSpace() == false &&
225 {
227 outarray[n], 1);
229 outarray[n], 1, outarray[n], 1);
231 outarray[n], 1, outarray[n], 1);
232 Vmath::Vmul(nPointsTot, inarray[n], 1, velocity[0], 1,
233 gradV0, 1);
234 Vmath::Vmul(nPointsTot, inarray[n], 1, velocity[1], 1,
235 gradV1, 1);
236 Vmath::Vmul(nPointsTot, inarray[n], 1, velocity[2], 1,
237 gradV2, 1);
239 gradV0, tmp);
240 Vmath::Vadd(nPointsTot, tmp, 1, outarray[n], 1, outarray[n],
241 1);
243 gradV1, tmp);
244 Vmath::Vadd(nPointsTot, tmp, 1, outarray[n], 1, outarray[n],
245 1);
247 gradV2, tmp);
248 Vmath::Vadd(nPointsTot, tmp, 1, outarray[n], 1, outarray[n],
249 1);
250 Vmath::Smul(nPointsTot, 0.5, outarray[n], 1, outarray[n],
251 1);
252 }
253 else
254 {
256 "Dealiasing is not allowed in combination "
257 "with the Skew-Symmetric advection form for "
258 "efficiency reasons.");
259 }
260 break;
261 default:
262 ASSERTL0(
false,
"dimension unknown");
263 }
264
266 }
267}
#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.