Advects a vector field.
72{
73
74 int ndim = advVel.size();
75 int nPointsTot = fields[0]->GetNpoints();
76 Array<OneD, Array<OneD, NekDouble>> velocity(ndim);
77 for (int i = 0; i < ndim; ++i)
78 {
80 {
81 velocity[i] = Array<OneD, NekDouble>(nPointsTot, 0.0);
82 fields[i]->HomogeneousBwdTrans(nPointsTot, advVel[i], velocity[i]);
83 }
84 else
85 {
86 velocity[i] = advVel[i];
87 }
88 }
89 for (int n = 0; n < nConvectiveFields; ++n)
90 {
91
92 Array<OneD, NekDouble> gradV0, gradV1, gradV2, tmp, Up;
93
94 gradV0 = Array<OneD, NekDouble>(nPointsTot);
95 tmp = Array<OneD, NekDouble>(nPointsTot);
96
97
98 switch (ndim)
99 {
100 case 1:
102 {
103 fields[0]->PhysDeriv(inarray[n], gradV0);
105 outarray[n], 1);
106 }
107 else
108 {
109 Vmath::Vmul(nPointsTot, inarray[n], 1, velocity[0], 1,
110 gradV0, 1);
111 fields[0]->PhysDeriv(gradV0, outarray[n]);
112 }
113 Vmath::Smul(nPointsTot, 0.5, outarray[n], 1, outarray[n],
114 1);
115 break;
116 case 2:
117 gradV1 = Array<OneD, NekDouble>(nPointsTot);
119 {
120 fields[0]->PhysDeriv(inarray[n], gradV0, gradV1);
122 outarray[n], 1);
124 outarray[n], 1, outarray[n], 1);
125 }
126 else
127 {
128 Vmath::Vmul(nPointsTot, inarray[n], 1, velocity[0], 1,
129 gradV0, 1);
130 Vmath::Vmul(nPointsTot, inarray[n], 1, velocity[1], 1,
131 gradV1, 1);
133 gradV0, outarray[n]);
135 gradV1, tmp);
136 Vmath::Vadd(nPointsTot, tmp, 1, outarray[n], 1, outarray[n],
137 1);
138 }
139 Vmath::Smul(nPointsTot, 1.0, outarray[n], 1, outarray[n],
140 1);
141 break;
142 case 3:
143 gradV1 = Array<OneD, NekDouble>(nPointsTot);
144 gradV2 = Array<OneD, NekDouble>(nPointsTot);
145
146
147
148
149 if (fields[0]->GetWaveSpace() == true)
150 {
152 {
153
154
155 fields[0]->PhysDeriv(inarray[n], gradV0, gradV1,
156 gradV2);
157 fields[0]->HomogeneousBwdTrans(nPointsTot, gradV0, tmp);
159 outarray[n], 1);
160 fields[0]->HomogeneousBwdTrans(nPointsTot, gradV1, tmp);
162 outarray[n], 1, outarray[n],
163 1);
164 fields[0]->HomogeneousBwdTrans(nPointsTot, gradV2, tmp);
166 outarray[n], 1, outarray[n],
167 1);
168 }
169 else
170 {
171 Up = Array<OneD, NekDouble>(nPointsTot);
172 fields[0]->HomogeneousBwdTrans(nPointsTot, inarray[n],
173 Up);
174 Vmath::Vmul(nPointsTot, Up, 1, velocity[0], 1, gradV0,
175 1);
176 Vmath::Vmul(nPointsTot, Up, 1, velocity[1], 1, gradV1,
177 1);
178 Vmath::Vmul(nPointsTot, Up, 1, velocity[2], 1, gradV2,
179 1);
180
181 fields[0]->SetWaveSpace(false);
183 gradV0, outarray[n]);
185 gradV1, tmp);
187 outarray[n], 1);
189 gradV2, tmp);
191 outarray[n], 1);
192 fields[0]->SetWaveSpace(true);
193 }
194
196 1);
197 fields[0]->HomogeneousFwdTrans(nPointsTot, tmp,
198 outarray[n]);
199 }
200 else
201 {
203 {
204 fields[0]->PhysDeriv(inarray[n], gradV0, gradV1,
205 gradV2);
207 outarray[n], 1);
209 outarray[n], 1, outarray[n], 1);
211 outarray[n], 1, outarray[n], 1);
212 }
213 else
214 {
215 Vmath::Vmul(nPointsTot, inarray[n], 1, velocity[0], 1,
216 gradV0, 1);
217 Vmath::Vmul(nPointsTot, inarray[n], 1, velocity[1], 1,
218 gradV1, 1);
219 Vmath::Vmul(nPointsTot, inarray[n], 1, velocity[2], 1,
220 gradV2, 1);
222 gradV0, outarray[n]);
224 gradV1, tmp);
226 outarray[n], 1);
228 gradV2, tmp);
230 outarray[n], 1);
231 }
232 Vmath::Smul(nPointsTot, 1.0, outarray[n], 1, outarray[n],
233 1);
234 }
235 break;
236 default:
237 ASSERTL0(
false,
"dimension unknown");
238 }
239 }
240}
#define ASSERTL0(condition, msg)
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 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.