Advects a vector field.
85{
86 boost::ignore_unused(time, pFwd, pBwd);
87
88
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");
93
94 Array<OneD, Array<OneD, NekDouble>> velocity(ndim);
95 for (int i = 0; i < ndim; ++i)
96 {
98 {
99 velocity[i] = Array<OneD, NekDouble>(nqtot, 0.0);
100 fields[i]->HomogeneousBwdTrans(nqtot, advVel[i], velocity[i]);
101 }
102 else
103 {
104 velocity[i] = advVel[i];
105 }
106 }
107
108 for (int n = 0; n < nConvectiveFields; ++n)
109 {
110
111
112 int nPointsTot = fields[0]->GetNpoints();
113 Array<OneD, NekDouble> gradV0, gradV1, gradV2, tmp, Up;
114
115 gradV0 = Array<OneD, NekDouble>(nPointsTot);
116 tmp = Array<OneD, NekDouble>(nPointsTot);
117
118
119 switch (ndim)
120 {
121 case 1:
122 fields[0]->PhysDeriv(inarray[n], gradV0);
123 Vmath::Vmul(nPointsTot, gradV0, 1, velocity[0], 1, outarray[n],
124 1);
125 Vmath::Vmul(nPointsTot, inarray[n], 1, velocity[0], 1, gradV0,
126 1);
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);
130 break;
131 case 2:
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],
135 1);
136 Vmath::Vvtvp(nPointsTot, gradV1, 1, velocity[1], 1, outarray[n],
137 1, outarray[n], 1);
138 Vmath::Vmul(nPointsTot, inarray[n], 1, velocity[0], 1, gradV0,
139 1);
140 Vmath::Vmul(nPointsTot, inarray[n], 1, velocity[1], 1, gradV1,
141 1);
143 tmp);
144 Vmath::Vadd(nPointsTot, tmp, 1, outarray[n], 1, outarray[n], 1);
146 tmp);
147 Vmath::Vadd(nPointsTot, tmp, 1, outarray[n], 1, outarray[n], 1);
148 Vmath::Smul(nPointsTot, 0.5, outarray[n], 1, outarray[n], 1);
149 break;
150 case 3:
151 gradV1 = Array<OneD, NekDouble>(nPointsTot);
152 gradV2 = Array<OneD, NekDouble>(nPointsTot);
153
154 fields[0]->PhysDeriv(inarray[n], gradV0, gradV1, gradV2);
155
156
157
158
160 fields[0]->GetWaveSpace() == false)
161 {
162 fields[0]->DealiasedProd(nPointsTot, velocity[0], gradV0,
163 gradV0);
164 fields[0]->DealiasedProd(nPointsTot, velocity[1], gradV1,
165 gradV1);
166 fields[0]->DealiasedProd(nPointsTot, velocity[2], gradV2,
167 gradV2);
168 Vmath::Vadd(nPointsTot, gradV0, 1, gradV1, 1, outarray[n],
169 1);
171 outarray[n], 1);
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);
179 gradV0, tmp);
180 Vmath::Vadd(nPointsTot, tmp, 1, outarray[n], 1, outarray[n],
181 1);
183 gradV1, tmp);
184 Vmath::Vadd(nPointsTot, tmp, 1, outarray[n], 1, outarray[n],
185 1);
187 gradV2, tmp);
188 Vmath::Vadd(nPointsTot, tmp, 1, outarray[n], 1, outarray[n],
189 1);
190 Vmath::Smul(nPointsTot, 0.5, outarray[n], 1, outarray[n],
191 1);
192 }
193 else if (fields[0]->GetWaveSpace() == true &&
195 {
196 Up = Array<OneD, NekDouble>(nPointsTot);
197
198
199 fields[0]->HomogeneousBwdTrans(nPointsTot, gradV0, tmp);
200 Vmath::Vmul(nPointsTot, tmp, 1, velocity[0], 1, outarray[n],
201 1);
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);
208
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);
213
214 fields[0]->SetWaveSpace(false);
216 gradV0, tmp);
217 Vmath::Vadd(nPointsTot, tmp, 1, outarray[n], 1, outarray[n],
218 1);
220 gradV1, tmp);
221 Vmath::Vadd(nPointsTot, tmp, 1, outarray[n], 1, outarray[n],
222 1);
224 gradV2, tmp);
225 Vmath::Vadd(nPointsTot, tmp, 1, outarray[n], 1, outarray[n],
226 1);
227 fields[0]->SetWaveSpace(true);
228
229 Vmath::Smul(nPointsTot, 0.5, outarray[n], 1, tmp, 1);
230 fields[0]->HomogeneousFwdTrans(nPointsTot, tmp,
231 outarray[n]);
232 }
233 else if (fields[0]->GetWaveSpace() == false &&
235 {
237 outarray[n], 1);
239 outarray[n], 1, outarray[n], 1);
241 outarray[n], 1, outarray[n], 1);
242 Vmath::Vmul(nPointsTot, inarray[n], 1, velocity[0], 1,
243 gradV0, 1);
244 Vmath::Vmul(nPointsTot, inarray[n], 1, velocity[1], 1,
245 gradV1, 1);
246 Vmath::Vmul(nPointsTot, inarray[n], 1, velocity[2], 1,
247 gradV2, 1);
249 gradV0, tmp);
250 Vmath::Vadd(nPointsTot, tmp, 1, outarray[n], 1, outarray[n],
251 1);
253 gradV1, tmp);
254 Vmath::Vadd(nPointsTot, tmp, 1, outarray[n], 1, outarray[n],
255 1);
257 gradV2, tmp);
258 Vmath::Vadd(nPointsTot, tmp, 1, outarray[n], 1, outarray[n],
259 1);
260 Vmath::Smul(nPointsTot, 0.5, outarray[n], 1, outarray[n],
261 1);
262 }
263 else
264 {
266 "Dealiasing is not allowed in combination "
267 "with the Skew-Symmetric advection form for "
268 "efficiency reasons.");
269 }
270 break;
271 default:
272 ASSERTL0(
false,
"dimension unknown");
273 }
274
276 }
277}
#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.