76 int nvariables = physarray.num_elements();
82 int nPointsTot =
m_fields[0]->GetTotPoints();
83 int nPointsTot_plane =
m_fields[0]->GetPlane(0)->GetTotPoints();
84 n_planes = nPointsTot/nPointsTot_plane;
87 int id2, id2_plane, e_max;
91 for(
int e = 0; e < e_max; ++e)
94 GetBndCondExpansions()[
m_bcRegion]->GetExp(e)->GetTotPoints();
96 GetBndCondExpansions()[
m_bcRegion]->GetPhys_Offset(e);
101 int m_offset_plane =
m_offset/n_planes;
103 int e_max_plane = e_max/n_planes;
104 int nTracePts_plane =
m_fields[0]->GetTrace()->GetNpoints();
106 int planeID = floor((e + 0.5 )/ e_max_plane );
107 e_plane = e - e_max_plane*planeID;
109 id2_plane =
m_fields[0]->GetTrace()->GetPhys_Offset(
111 GetBndCondCoeffsToGlobalCoeffsMap(
112 m_offset_plane + e_plane));
113 id2 = id2_plane + planeID*nTracePts_plane;
118 GetTrace()->GetPhys_Offset(
m_fields[0]->GetTraceMap()->
119 GetBndCondTraceToGlobalTraceMap(
m_offset+e));
122 Array<OneD,NekDouble> x0(npoints, 0.0);
123 Array<OneD,NekDouble> x1(npoints, 0.0);
124 Array<OneD,NekDouble> x2(npoints, 0.0);
127 GetExp(e)->GetCoords(x0, x1, x2);
130 NekDouble c, k, phi, r, J, VV, pp, sint,
P, ss;
144 NekDouble gamma_1_2 = (gamma - 1.0) / 2.0;
147 for (
int j = 0; j < npoints; j++)
150 while ((abs(errV) > toll) || (abs(errTheta) > toll))
154 c = sqrt(1.0 - gamma_1_2 * VV);
158 J = 1.0 / c + 1.0 / (3.0 * c * c * c) +
159 1.0 / (5.0 * c * c * c * c * c) -
160 0.5 * log((1.0 + c) / (1.0 - c));
162 r = pow(c, 1.0 / gamma_1_2);
163 xi = 1.0 / (2.0 * r) * (1.0 / VV - 2.0 * pp) + J / 2.0;
164 yi = phi / (r * V) * sqrt(1.0 - VV * pp);
165 par1 = 25.0 - 5.0 * VV;
171 J11 = 39062.5 / pow(par1, 3.5) *
172 (1.0 / VV - 2.0 / VV * ss) * V + 1562.5 /
173 pow(par1, 2.5) * (-2.0 / (VV * V) + 4.0 /
174 (VV * V) * ss) + 12.5 / pow(par1, 1.5) * V +
175 312.5 / pow(par1, 2.5) * V + 7812.5 /
176 pow(par1, 3.5) * V - 0.25 *
177 (-1.0 / pow(par1, 0.5) * V / (1.0 - 0.2 *
178 pow(par1, 0.5)) - (1.0 + 0.2 * pow(par1, 0.5)) /
179 pow((1.0 - 0.2 * pow(par1, 0.5)), 2.0) /
180 pow(par1, 0.5) * V) / (1.0 + 0.2 * pow(par1, 0.5)) *
181 (1.0 - 0.2 * pow(par1, 0.5));
183 J12 = -6250.0 / pow(par1, 2.5) / VV * sint * cos(theta);
184 J21 = -6250.0 / (VV * V) * sint / pow(par1, 2.5) *
185 pow((1.0 - ss), 0.5) + 78125.0 / V * sint /
186 pow(par1, 3.5) * pow((1.0 - ss), 0.5);
189 if (abs(x1[j]) < toll && abs(cos(theta)) < toll)
191 J22 = -39062.5 / pow(par1, 3.5) / V + 3125 /
192 pow(par1, 2.5) / (VV * V) + 12.5 /
193 pow(par1, 1.5) * V + 312.5 / pow(par1, 2.5) *
194 V + 7812.5 / pow(par1, 3.5) * V - 0.25 *
195 (-1.0 / pow(par1, 0.5) * V / (1.0 - 0.2 *
196 pow(par1, 0.5)) - (1.0 + 0.2 * pow(par1, 0.5)) /
197 pow((1.0 - 0.2 * pow(par1, 0.5)), 2.0) /
198 pow(par1, 0.5) * V) / (1.0 + 0.2 *
199 pow(par1, 0.5)) * (1.0 - 0.2 * pow(par1, 0.5));
202 dV = -1.0 / J22 * Fx;
208 J22 = 3125.0 / VV * cos(theta) / pow(par1, 2.5) *
209 pow((1.0 - ss), 0.5) - 3125.0 / VV * ss /
210 pow(par1, 2.5) / pow((1.0 - ss), 0.5) *
213 det = -1.0 / (J11 * J22 - J12 * J21);
216 dV = det * ( J22 * Fx - J12 * Fy);
217 dtheta = det * (-J21 * Fx + J11 * Fy);
221 theta = theta + dtheta;
224 errTheta = abs(dtheta);
227 c = sqrt(1.0 - gamma_1_2 * VV);
231 !(
m_session->DefinesFunction(
"InitialConditions") &&
232 m_session->GetFunctionType(
"InitialConditions", 0) ==
235 Fwd[0][kk] = pow(c, 1.0 / gamma_1_2) *
236 exp(-1.0 + time /timeramp);
238 Fwd[1][kk] = Fwd[0][kk] * V * cos(theta) *
239 exp(-1 + time / timeramp);
241 Fwd[2][kk] = Fwd[0][kk] * V * sin(theta) *
242 exp(-1 + time / timeramp);
246 Fwd[0][kk] = pow(c, 1.0 / gamma_1_2);
247 Fwd[1][kk] = Fwd[0][kk] * V * cos(theta);
248 Fwd[2][kk] = Fwd[0][kk] * V * sin(theta);
251 P = (c * c) * Fwd[0][kk] / gamma;
252 Fwd[3][kk] = P / (gamma - 1.0) + 0.5 *
253 (Fwd[1][kk] * Fwd[1][kk] / Fwd[0][kk] +
254 Fwd[2][kk] * Fwd[2][kk] / Fwd[0][kk]);
259 V = kExt * sin(theta);
262 for (
int i = 0; i < nvariables; ++i)
266 UpdatePhys())[id1],1);
LibUtilities::SessionReaderSharedPtr m_session
Session reader.
Array< OneD, MultiRegions::ExpListSharedPtr > m_fields
Array of fields.
NekDouble m_gamma
Parameters of the flow.
void Vcopy(int n, const T *x, const int incx, T *y, const int incy)
int m_bcRegion
Id of the boundary region.