Nektar++
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
EulerCFE.cpp
Go to the documentation of this file.
1 ///////////////////////////////////////////////////////////////////////////////
2 //
3 // File EulerCFE.cpp
4 //
5 // For more information, please see: http://www.nektar.info
6 //
7 // The MIT License
8 //
9 // Copyright (c) 2006 Division of Applied Mathematics, Brown University (USA),
10 // Department of Aeronautics, Imperial College London (UK), and Scientific
11 // Computing and Imaging Institute, University of Utah (USA).
12 //
13 // License for the specific language governing rights and limitations under
14 // Permission is hereby granted, free of charge, to any person obtaining a
15 // copy of this software and associated documentation files (the "Software"),
16 // to deal in the Software without restriction, including without limitation
17 // the rights to use, copy, modify, merge, publish, distribute, sublicense,
18 // and/or sell copies of the Software, and to permit persons to whom the
19 // Software is furnished to do so, subject to the following conditions:
20 //
21 // The above copyright notice and this permission notice shall be included
22 // in all copies or substantial portions of the Software.
23 //
24 // THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
25 // OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
26 // FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
27 // THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
28 // LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
29 // FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
30 // DEALINGS IN THE SOFTWARE.
31 //
32 // Description: Euler equations in consƒervative variables without artificial
33 // diffusion
34 //
35 ///////////////////////////////////////////////////////////////////////////////
36 
37 #include <iostream>
38 #include <iomanip>
39 #include <boost/algorithm/string.hpp>
40 
41 #include <MultiRegions/AssemblyMap/AssemblyMapDG.h>
42 #include <CompressibleFlowSolver/EquationSystems/EulerCFE.h>
43 
44 namespace Nektar
45 {
46  string EulerCFE::className =
47  SolverUtils::GetEquationSystemFactory().RegisterCreatorFunction(
48  "EulerCFE", EulerCFE::create,
49  "Euler equations in conservative variables without "
50  "artificial diffusion.");
51 
52  EulerCFE::EulerCFE(
53  const LibUtilities::SessionReaderSharedPtr& pSession)
54  : CompressibleFlowSystem(pSession)
55  {
56  }
57 
58  void EulerCFE::v_InitObject()
59  {
60  CompressibleFlowSystem::v_InitObject();
61 
62  if(m_session->DefinesSolverInfo("PROBLEMTYPE"))
63  {
64  int i;
65  std::string ProblemTypeStr =
66  m_session->GetSolverInfo("PROBLEMTYPE");
67  for (i = 0; i < (int) SIZE_ProblemType; ++i)
68  {
69  if (boost::iequals(ProblemTypeMap[i], ProblemTypeStr))
70  {
71  m_problemType = (ProblemType)i;
72  break;
73  }
74  }
75  }
76  else
77  {
78  m_problemType = (ProblemType)0;
79  }
80 
81  if (m_explicitAdvection)
82  {
83  m_ode.DefineOdeRhs (&EulerCFE::DoOdeRhs, this);
84  m_ode.DefineProjection (&EulerCFE::DoOdeProjection, this);
85  }
86  else
87  {
88  ASSERTL0(false, "Implicit CFE not set up.");
89  }
90  }
91 
92  /**
93  * @brief Destructor for EulerCFE class.
94  */
95  EulerCFE::~EulerCFE()
96  {
97  }
98 
99  /**
100  * @brief Print out a summary with some relevant information.
101  */
102  void EulerCFE::v_GenerateSummary(SolverUtils::SummaryList& s)
103  {
104  CompressibleFlowSystem::v_GenerateSummary(s);
105  SolverUtils::AddSummaryItem(s, "Problem Type", ProblemTypeMap[m_problemType]);
106  }
107 
108  /**
109  * @brief Set the initial conditions.
110  */
111  void EulerCFE::v_SetInitialConditions(
112  NekDouble initialtime,
113  bool dumpInitialConditions,
114  const int domain)
115  {
116  switch (m_problemType)
117  {
118  case eIsentropicVortex:
119  {
120  SetInitialIsentropicVortex(initialtime);
121  break;
122  }
123  default:
124  {
125  EquationSystem::v_SetInitialConditions(initialtime, false);
126  break;
127  }
128  }
129 
130  //insert white noise in initial condition
131  NekDouble Noise;
132  int phystot = m_fields[0]->GetTotPoints();
133  Array<OneD, NekDouble> noise(phystot);
134 
135  m_session->LoadParameter("Noise", Noise,0.0);
136  int m_nConvectiveFields = m_fields.num_elements();
137 
138  if(Noise > 0.0)
139  {
140  for(int i = 0; i < m_nConvectiveFields; i++)
141  {
142  Vmath::FillWhiteNoise(phystot,Noise,noise,1,m_comm->GetColumnComm()->GetRank()+1);
143  Vmath::Vadd(phystot,m_fields[i]->GetPhys(),1,noise,1,m_fields[i]->UpdatePhys(),1);
144  m_fields[i]->FwdTrans_IterPerExp(m_fields[i]->GetPhys(),m_fields[i]->UpdateCoeffs());
145  }
146  }
147 
148  if (dumpInitialConditions)
149  {
150  // Dump initial conditions to file
151  Checkpoint_Output(0);
152  }
153  }
154 
155  /**
156  * @brief Compute the right-hand side.
157  */
158  void EulerCFE::DoOdeRhs(
159  const Array<OneD, const Array<OneD, NekDouble> > &inarray,
160  Array<OneD, Array<OneD, NekDouble> > &outarray,
161  const NekDouble time)
162  {
163  int i;
164  int nvariables = inarray.num_elements();
165  int npoints = GetNpoints();
166 
167  Array<OneD, Array<OneD, NekDouble> > advVel(m_spacedim);
168 
169  m_advection->Advect(nvariables, m_fields, advVel, inarray,
170  outarray, time);
171 
172  for (i = 0; i < nvariables; ++i)
173  {
174  Vmath::Neg(npoints, outarray[i], 1);
175  }
176  }
177 
178  /**
179  * @brief Compute the projection and call the method for imposing the
180  * boundary conditions in case of discontinuous projection.
181  */
182  void EulerCFE::DoOdeProjection(
183  const Array<OneD, const Array<OneD, NekDouble> > &inarray,
184  Array<OneD, Array<OneD, NekDouble> > &outarray,
185  const NekDouble time)
186  {
187  int i;
188  int nvariables = inarray.num_elements();
189 
190  switch (m_projectionType)
191  {
192  case MultiRegions::eDiscontinuous:
193  {
194  // Just copy over array
195  int npoints = GetNpoints();
196 
197  for(i = 0; i < nvariables; ++i)
198  {
199  Vmath::Vcopy(npoints, inarray[i], 1, outarray[i], 1);
200  }
201  SetBoundaryConditions(outarray, time);
202  break;
203  }
204  case MultiRegions::eGalerkin:
205  case MultiRegions::eMixed_CG_Discontinuous:
206  {
207  ASSERTL0(false, "No Continuous Galerkin for Euler equations");
208  break;
209  }
210  default:
211  ASSERTL0(false, "Unknown projection scheme");
212  break;
213  }
214  }
215 
216  /**
217  * @brief Set boundary conditions which can be:
218  * a) Wall and Symmerty BCs implemented at CompressibleFlowSystem level
219  * since they are compressible solver specific;
220  * b) Isentropic vortex and Ringleb flow BCs implemented at EulerCFE level
221  * since they are Euler solver specific;
222  * c) Time dependent BCs.
223  *
224  * @param inarray: fields array.
225  * @param time: time.
226  */
227  void EulerCFE::SetBoundaryConditions(
228  Array<OneD, Array<OneD, NekDouble> > &inarray,
229  NekDouble time)
230  {
231  std::string varName;
232  int nvariables = m_fields.num_elements();
233  int cnt = 0;
234 
235  // Loop over Boundary Regions
236  for (int n = 0; n < m_fields[0]->GetBndConditions().num_elements(); ++n)
237  {
238  // Wall Boundary Condition
239  if (m_fields[0]->GetBndConditions()[n]->GetUserDefined() ==
240  SpatialDomains::eWall)
241  {
242  WallBC(n, cnt, inarray);
243  }
244 
245  // Wall Boundary Condition
246  if (m_fields[0]->GetBndConditions()[n]->GetUserDefined() ==
247  SpatialDomains::eWallViscous)
248  {
249  ASSERTL0(false, "WallViscous is a wrong bc for the "
250  "Euler equations");
251  }
252 
253  // Symmetric Boundary Condition
254  if (m_fields[0]->GetBndConditions()[n]->GetUserDefined() ==
255  SpatialDomains::eSymmetry)
256  {
257  SymmetryBC(n, cnt, inarray);
258  }
259 
260  // Riemann invariant characteristic Boundary Condition (CBC)
261  if (m_fields[0]->GetBndConditions()[n]->GetUserDefined() ==
262  SpatialDomains::eRiemannInvariant)
263  {
264  RiemannInvariantBC(n, cnt, inarray);
265  }
266 
267  // Extrapolation of the data at the boundaries
268  if (m_fields[0]->GetBndConditions()[n]->GetUserDefined() ==
269  SpatialDomains::eExtrapOrder0)
270  {
271  ExtrapOrder0BC(n, cnt, inarray);
272  }
273 
274  // Time Dependent Boundary Condition (specified in meshfile)
275  if (m_fields[0]->GetBndConditions()[n]->GetUserDefined() ==
276  SpatialDomains::eTimeDependent)
277  {
278  for (int i = 0; i < nvariables; ++i)
279  {
280  varName = m_session->GetVariable(i);
281  m_fields[i]->EvaluateBoundaryConditions(time, varName);
282  }
283  }
284 
285  // Isentropic Vortex Boundary Condition
286  if (m_fields[0]->GetBndConditions()[n]->GetUserDefined() ==
287  SpatialDomains::eIsentropicVortex)
288  {
289  SetBoundaryIsentropicVortex(n, time, cnt, inarray);
290  }
291 
292  // Ringleb Flow Inflow and outflow Boundary Condition
293  if (m_fields[0]->GetBndConditions()[n]->GetUserDefined() ==
294  SpatialDomains::eRinglebFlow)
295  {
296  SetBoundaryRinglebFlow(n, time, cnt, inarray);
297  }
298 
299  cnt += m_fields[0]->GetBndCondExpansions()[n]->GetExpSize();
300  }
301  }
302 
303  /**
304  * @brief Get the exact solutions for isentropic vortex and Ringleb
305  * flow problems.
306  */
307  void EulerCFE::v_EvaluateExactSolution(
308  unsigned int field,
309  Array<OneD, NekDouble> &outfield,
310  const NekDouble time)
311  {
312  switch(m_problemType)
313  {
314  case eIsentropicVortex:
315  {
316  GetExactIsentropicVortex(field, outfield, time);
317  break;
318  }
319  case eRinglebFlow:
320  {
321  GetExactRinglebFlow(field, outfield);
322  break;
323  }
324  default:
325  {
326  break;
327  }
328  }
329  }
330 
331  void EulerCFE::EvaluateIsentropicVortex(
332  const Array<OneD, NekDouble> &x,
333  const Array<OneD, NekDouble> &y,
334  const Array<OneD, NekDouble> &z,
335  Array<OneD, Array<OneD, NekDouble> > &u,
336  NekDouble time,
337  const int o)
338  {
339  int nq = x.num_elements();
340 
341  // Flow parameters
342  const NekDouble x0 = 5.0;
343  const NekDouble y0 = 0.0;
344  const NekDouble beta = 5.0;
345  const NekDouble u0 = 1.0;
346  const NekDouble v0 = 0.5;
347  const NekDouble gamma = m_gamma;
348  NekDouble r, xbar, ybar, tmp;
349  NekDouble fac = 1.0/(16.0*gamma*M_PI*M_PI);
350 
351  // In 3D zero rhow field.
352  if (m_spacedim == 3)
353  {
354  Vmath::Zero(nq, &u[3][o], 1);
355  }
356 
357  // Fill storage
358  for (int i = 0; i < nq; ++i)
359  {
360  xbar = x[i] - u0*time - x0;
361  ybar = y[i] - v0*time - y0;
362  r = sqrt(xbar*xbar + ybar*ybar);
363  tmp = beta*exp(1-r*r);
364  u[0][i+o] = pow(1.0 - (gamma-1.0)*tmp*tmp*fac, 1.0/(gamma-1.0));
365  u[1][i+o] = u[0][i+o]*(u0 - tmp*ybar/(2*M_PI));
366  u[2][i+o] = u[0][i+o]*(v0 + tmp*xbar/(2*M_PI));
367  u[m_spacedim+1][i+o] = pow(u[0][i+o], gamma)/(gamma-1.0) +
368  0.5*(u[1][i+o]*u[1][i+o] + u[2][i+o]*u[2][i+o]) / u[0][i+o];
369  }
370  }
371 
372  /**
373  * @brief Compute the exact solution for the isentropic vortex problem.
374  */
375  void EulerCFE::GetExactIsentropicVortex(
376  int field,
377  Array<OneD, NekDouble> &outarray,
378  NekDouble time)
379  {
380  int nTotQuadPoints = GetTotPoints();
381  Array<OneD, NekDouble> x(nTotQuadPoints);
382  Array<OneD, NekDouble> y(nTotQuadPoints);
383  Array<OneD, NekDouble> z(nTotQuadPoints);
384  Array<OneD, Array<OneD, NekDouble> > u(m_spacedim+2);
385 
386  m_fields[0]->GetCoords(x, y, z);
387 
388  for (int i = 0; i < m_spacedim + 2; ++i)
389  {
390  u[i] = Array<OneD, NekDouble>(nTotQuadPoints);
391  }
392 
393  EvaluateIsentropicVortex(x, y, z, u, time);
394 
395  Vmath::Vcopy(nTotQuadPoints, u[field], 1, outarray, 1);
396  }
397 
398  /**
399  * @brief Set the initial condition for the isentropic vortex problem.
400  */
401  void EulerCFE::SetInitialIsentropicVortex(NekDouble initialtime)
402  {
403  int nTotQuadPoints = GetTotPoints();
404  Array<OneD, NekDouble> x(nTotQuadPoints);
405  Array<OneD, NekDouble> y(nTotQuadPoints);
406  Array<OneD, NekDouble> z(nTotQuadPoints);
407  Array<OneD, Array<OneD, NekDouble> > u(m_spacedim+2);
408 
409  m_fields[0]->GetCoords(x, y, z);
410 
411  for (int i = 0; i < m_spacedim + 2; ++i)
412  {
413  u[i] = Array<OneD, NekDouble>(nTotQuadPoints);
414  }
415 
416  EvaluateIsentropicVortex(x, y, z, u, initialtime);
417 
418  // Forward transform to fill the coefficient space
419  for(int i = 0; i < m_fields.num_elements(); ++i)
420  {
421  Vmath::Vcopy(nTotQuadPoints, u[i], 1, m_fields[i]->UpdatePhys(), 1);
422  m_fields[i]->SetPhysState(true);
423  m_fields[i]->FwdTrans(m_fields[i]->GetPhys(),
424  m_fields[i]->UpdateCoeffs());
425  }
426  }
427 
428  /**
429  * @brief Set the boundary conditions for the isentropic vortex problem.
430  */
431  void EulerCFE::SetBoundaryIsentropicVortex(
432  int bcRegion,
433  NekDouble time,
434  int cnt,
435  Array<OneD, Array<OneD, NekDouble> > &physarray)
436  {
437  int nvariables = physarray.num_elements();
438  int nTraceNumPoints = GetTraceTotPoints();
439  Array<OneD, Array<OneD, NekDouble> > Fwd(nvariables);
440  const Array<OneD, const int> &bndTraceMap = m_fields[0]->GetTraceBndMap();
441 
442  // Get physical values of the forward trace (from exp to phys)
443  for (int i = 0; i < nvariables; ++i)
444  {
445  Fwd[i] = Array<OneD, NekDouble>(nTraceNumPoints);
446  m_fields[i]->ExtractTracePhys(physarray[i], Fwd[i]);
447  }
448 
449  int id2, e_max;
450  e_max = m_fields[0]->GetBndCondExpansions()[bcRegion]->GetExpSize();
451 
452  for(int e = 0; e < e_max; ++e)
453  {
454  int npoints = m_fields[0]->
455  GetBndCondExpansions()[bcRegion]->GetExp(e)->GetTotPoints();
456  int id1 = m_fields[0]->
457  GetBndCondExpansions()[bcRegion]->GetPhys_Offset(e);
458  id2 = m_fields[0]->GetTrace()->GetPhys_Offset(bndTraceMap[cnt++]);
459 
460  Array<OneD,NekDouble> x(npoints, 0.0);
461  Array<OneD,NekDouble> y(npoints, 0.0);
462  Array<OneD,NekDouble> z(npoints, 0.0);
463 
464  m_fields[0]->GetBndCondExpansions()[bcRegion]->
465  GetExp(e)->GetCoords(x, y, z);
466 
467  EvaluateIsentropicVortex(x, y, z, Fwd, time, id2);
468 
469  for (int i = 0; i < nvariables; ++i)
470  {
471  Vmath::Vcopy(npoints, &Fwd[i][id2], 1,
472  &(m_fields[i]->GetBndCondExpansions()[bcRegion]->
473  UpdatePhys())[id1], 1);
474  }
475  }
476  }
477 
478  /**
479  * @brief Compute the exact solution for the Ringleb flow problem.
480  */
481  void EulerCFE::GetExactRinglebFlow(
482  int field,
483  Array<OneD, NekDouble> &outarray)
484  {
485  int nTotQuadPoints = GetTotPoints();
486 
487  Array<OneD, NekDouble> rho(nTotQuadPoints, 100.0);
488  Array<OneD, NekDouble> rhou(nTotQuadPoints);
489  Array<OneD, NekDouble> rhov(nTotQuadPoints);
490  Array<OneD, NekDouble> E(nTotQuadPoints);
491  Array<OneD, NekDouble> x(nTotQuadPoints);
492  Array<OneD, NekDouble> y(nTotQuadPoints);
493  Array<OneD, NekDouble> z(nTotQuadPoints);
494 
495  m_fields[0]->GetCoords(x, y, z);
496 
497  // Flow parameters
498  NekDouble c, k, phi, r, J, VV, pp, sint, P, ss;
499  NekDouble J11, J12, J21, J22, det;
500  NekDouble Fx, Fy;
501  NekDouble xi, yi;
502  NekDouble dV;
503  NekDouble dtheta;
504  NekDouble par1;
505  NekDouble theta = M_PI / 4.0;
506  NekDouble kExt = 0.7;
507  NekDouble V = kExt * sin(theta);
508  NekDouble toll = 1.0e-8;
509  NekDouble errV = 1.0;
510  NekDouble errTheta = 1.0;
511  NekDouble gamma = m_gamma;
512  NekDouble gamma_1_2 = (gamma - 1.0) / 2.0;
513 
514  for (int i = 0; i < nTotQuadPoints; ++i)
515  {
516  while ((abs(errV) > toll) || (abs(errTheta) > toll))
517  {
518  VV = V * V;
519  sint = sin(theta);
520  c = sqrt(1.0 - gamma_1_2 * VV);
521  k = V / sint;
522  phi = 1.0 / k;
523  pp = phi * phi;
524  J = 1.0 / c + 1.0 / (3.0 * c * c * c) +
525  1.0 / (5.0 * c * c * c * c * c) -
526  0.5 * log((1.0 + c) / (1.0 - c));
527 
528  r = pow(c, 1.0 / gamma_1_2);
529  xi = 1.0 / (2.0 * r) * (1.0 / VV - 2.0 * pp) + J / 2.0;
530  yi = phi / (r * V) * sqrt(1.0 - VV * pp);
531  par1 = 25.0 - 5.0 * VV;
532  ss = sint * sint;
533 
534  Fx = xi - x[i];
535  Fy = yi - y[i];
536 
537  J11 = 39062.5 / pow(par1, 3.5) * (1.0 / VV - 2.0 / VV * ss) *
538  V + 1562.5 / pow(par1, 2.5) * (-2.0 / (VV * V) + 4.0 /
539  (VV * V) * ss) + 12.5 / pow(par1, 1.5) * V + 312.5 /
540  pow(par1, 2.5) * V + 7812.5 / pow(par1, 3.5) * V -
541  0.25 * (-1.0 / pow(par1, 0.5) * V/(1.0 - 0.2 *
542  pow(par1, 0.5)) - (1.0 + 0.2 * pow(par1, 0.5)) /
543  pow((1.0 - 0.2 * pow(par1, 0.5)), 2.0) /
544  pow(par1, 0.5) * V) / (1.0 + 0.2 * pow(par1, 0.5)) *
545  (1.0 - 0.2 * pow(par1, 0.5));
546 
547  J12 = -6250.0 / pow(par1, 2.5) / VV * sint * cos(theta);
548  J21 = -6250.0 / (VV * V) * sint /
549  pow(par1, 2.5) * pow((1.0 - ss), 0.5) +
550  78125.0 / V * sint / pow(par1, 3.5) *
551  pow((1.0 - ss), 0.5);
552 
553  // the matrix is singular when theta = pi/2
554  if(abs(y[i])<toll && abs(cos(theta))<toll)
555  {
556  J22 = -39062.5 / pow(par1, 3.5) / V + 3125 /
557  pow(par1, 2.5) / (VV * V) + 12.5 / pow(par1, 1.5) *
558  V + 312.5 / pow(par1, 2.5) * V + 7812.5 /
559  pow(par1, 3.5) * V - 0.25 * (-1.0 / pow(par1, 0.5) *
560  V / (1.0 - 0.2 * pow(par1, 0.5)) - (1.0 + 0.2 *
561  pow(par1, 0.5)) / pow((1.0 - 0.2 *
562  pow(par1, 0.5)), 2.0) / pow(par1, 0.5) * V) /
563  (1.0 + 0.2 * pow(par1, 0.5)) * (1.0 - 0.2 *
564  pow(par1,0.5));
565 
566  // dV = -dV/dx * Fx
567  dV = -1.0 / J22 * Fx;
568  dtheta = 0.0;
569  theta = M_PI / 2.0;
570  }
571  else
572  {
573  J22 = 3125.0 / VV * cos(theta) / pow(par1, 2.5) *
574  pow((1.0 - ss), 0.5) - 3125.0 / VV * ss /
575  pow(par1, 2.5) / pow((1.0 - ss), 0.5) * cos(theta);
576 
577  det = -1.0 / (J11 * J22 - J12 * J21);
578 
579  // [dV dtheta]' = -[invJ]*[Fx Fy]'
580  dV = det * ( J22 * Fx - J12 * Fy);
581  dtheta = det * (-J21 * Fx + J11 * Fy);
582  }
583 
584  V = V + dV;
585  theta = theta + dtheta;
586 
587  errV = abs(dV);
588  errTheta = abs(dtheta);
589 
590  }
591 
592  c = sqrt(1.0 - gamma_1_2 * VV);
593  r = pow(c, 1.0 / gamma_1_2);
594 
595  rho[i] = r;
596  rhou[i] = rho[i] * V * cos(theta);
597  rhov[i] = rho[i] * V * sin(theta);
598  P = (c * c) * rho[i] / gamma;
599  E[i] = P / (gamma - 1.0) + 0.5 *
600  (rhou[i] * rhou[i] / rho[i] + rhov[i] * rhov[i] / rho[i]);
601 
602  // Resetting the guess value
603  errV = 1.0;
604  errTheta = 1.0;
605  theta = M_PI/4.0;
606  V = kExt*sin(theta);
607  }
608 
609  switch (field)
610  {
611  case 0:
612  outarray = rho;
613  break;
614  case 1:
615  outarray = rhou;
616  break;
617  case 2:
618  outarray = rhov;
619  break;
620  case 3:
621  outarray = E;
622  break;
623  default:
624  ASSERTL0(false, "Error in variable number!");
625  break;
626  }
627  }
628 
629  /**
630  * @brief Set the initial condition for the Ringleb flow problem.
631  */
632  void EulerCFE::SetInitialRinglebFlow(void)
633  {
634  // Get number of different boundaries in the input file
635  int nbnd = m_fields[0]->GetBndConditions().num_elements();
636 
637  // Loop on all the edges of the input file
638  for(int bcRegion=0; bcRegion < nbnd; ++bcRegion)
639  {
640 
641  int npoints = m_fields[0]->
642  GetBndCondExpansions()[bcRegion]->GetNpoints();
643 
644  Array<OneD,NekDouble> x0(npoints, 0.0);
645  Array<OneD,NekDouble> x1(npoints, 0.0);
646  Array<OneD,NekDouble> x2(npoints, 0.0);
647 
648  Array<OneD, NekDouble> rho(npoints, 0.0);
649  Array<OneD, NekDouble> rhou(npoints, 0.0);
650  Array<OneD, NekDouble> rhov(npoints, 0.0);
651  Array<OneD, NekDouble> E(npoints, 0.0);
652 
653  m_fields[0]->GetBndCondExpansions()[bcRegion]->
654  GetCoords(x0, x1, x2);
655 
656  // Flow parameters
657  NekDouble c, k, phi, r, J, VV, pp, sint, P, ss;
658  NekDouble J11, J12, J21, J22, det;
659  NekDouble Fx, Fy;
660  NekDouble xi, yi;
661  NekDouble dV;
662  NekDouble dtheta;
663  NekDouble par1;
664  NekDouble theta = M_PI / 4.0;
665  NekDouble kExt = 0.7;
666  NekDouble V = kExt * sin(theta);
667  NekDouble toll = 1.0e-8;
668  NekDouble errV = 1.0;
669  NekDouble errTheta = 1.0;
670  NekDouble gamma = m_gamma;
671  NekDouble gamma_1_2 = (gamma - 1.0) / 2.0;
672 
673  // Loop on all the points of that edge
674  for (int j = 0; j < npoints; j++)
675  {
676  while ((abs(errV) > toll) || (abs(errTheta) > toll))
677  {
678 
679  VV = V * V;
680  sint = sin(theta);
681  c = sqrt(1.0 - gamma_1_2 * VV);
682  k = V / sint;
683  phi = 1.0 / k;
684  pp = phi * phi;
685  J = 1.0 / c + 1.0 / (3.0 * c * c * c) +
686  1.0 / (5.0 * c * c * c * c * c) -
687  0.5 * log((1.0 + c) / (1.0 - c));
688 
689  r = pow(c, 1.0 / gamma_1_2);
690  xi = 1.0 / (2.0 * r) * (1.0 / VV - 2.0 * pp) + J / 2.0;
691  yi = phi / (r * V) * sqrt(1.0 - VV * pp);
692  par1 = 25.0 - 5.0 * VV;
693  ss = sint * sint;
694 
695  Fx = xi - x0[j];
696  Fy = yi - x1[j];
697 
698  J11 = 39062.5 / pow(par1, 3.5) *
699  (1.0 / VV - 2.0 / VV * ss) * V +
700  1562.5 / pow(par1, 2.5) * (-2.0 /
701  (VV * V) + 4.0 / (VV * V) * ss) +
702  12.5 / pow(par1, 1.5) * V +
703  312.5 / pow(par1, 2.5) * V +
704  7812.5 / pow(par1, 3.5) * V -
705  0.25 * (-1.0 / pow(par1, 0.5) * V /
706  (1.0 - 0.2 * pow(par1, 0.5)) - (1.0 + 0.2 *
707  pow(par1, 0.5)) / pow((1.0 - 0.2 *
708  pow(par1, 0.5)), 2.0) /
709  pow(par1, 0.5) * V) /
710  (1.0 + 0.2 * pow(par1, 0.5)) *
711  (1.0 - 0.2 * pow(par1, 0.5));
712 
713  J12 = -6250.0 / pow(par1, 2.5) / VV * sint * cos(theta);
714  J21 = -6250.0 / (VV * V) * sint / pow(par1, 2.5) *
715  pow((1.0 - ss), 0.5) + 78125.0 / V * sint /
716  pow(par1, 3.5) * pow((1.0 - ss), 0.5);
717 
718  // the matrix is singular when theta = pi/2
719  if (abs(x1[j]) < toll && abs(cos(theta)) < toll)
720  {
721 
722  J22 = -39062.5 / pow(par1, 3.5) / V +
723  3125 / pow(par1, 2.5) / (VV * V) + 12.5 /
724  pow(par1, 1.5) * V + 312.5 / pow(par1, 2.5) * V +
725  7812.5 / pow(par1, 3.5) * V -
726  0.25 * (-1.0 / pow(par1, 0.5) * V /
727  (1.0 - 0.2 * pow(par1, 0.5)) -
728  (1.0 + 0.2 * pow(par1, 0.5)) /
729  pow((1.0 - 0.2* pow(par1, 0.5)), 2.0) /
730  pow(par1, 0.5) *V) /
731  (1.0 + 0.2 * pow(par1, 0.5)) *
732  (1.0 - 0.2 * pow(par1, 0.5));
733 
734  // dV = -dV/dx * Fx
735  dV = -1.0 / J22 * Fx;
736  dtheta = 0.0;
737  theta = M_PI / 2.0;
738  }
739  else
740  {
741 
742  J22 = 3125.0 / VV * cos(theta) / pow(par1, 2.5) *
743  pow((1.0 - ss), 0.5) - 3125.0 / VV * ss /
744  pow(par1, 2.5) / pow((1.0 - ss), 0.5) *
745  cos(theta);
746 
747  det = -1.0 / (J11 * J22 - J12 * J21);
748 
749  // [dV dtheta]' = -[invJ]*[Fx Fy]'
750  dV = det * ( J22 * Fx - J12 * Fy);
751  dtheta = det * (-J21 * Fx + J11 * Fy);
752  }
753 
754  V = V + dV;
755  theta = theta + dtheta;
756 
757  errV = abs(dV);
758  errTheta = abs(dtheta);
759 
760  }
761 
762  c = sqrt(1.0 - gamma_1_2 * VV);
763  rho[j] = pow(c, 1.0 / gamma_1_2) * exp(-1.0);
764  rhou[j] = rho[j] * V * cos(theta) * exp(-1.0);
765  rhov[j] = rho[j] * V * sin(theta) * exp(-1.0);
766  P = (c * c) * rho[j] / gamma;
767  E[j] = P / (gamma - 1.0) + 0.5 *
768  (rhou[j] * rhou[j] /
769  rho[j] + rhov[j] * rhov[j] / rho[j]);
770  errV = 1.0;
771  errTheta = 1.0;
772  theta = M_PI / 4.0;
773  V = kExt * sin(theta);
774 
775  }
776 
777  // Fill the physical space
778  m_fields[0]->GetBndCondExpansions()[bcRegion]->SetPhys(rho);
779  m_fields[1]->GetBndCondExpansions()[bcRegion]->SetPhys(rhou);
780  m_fields[2]->GetBndCondExpansions()[bcRegion]->SetPhys(rhov);
781  m_fields[3]->GetBndCondExpansions()[bcRegion]->SetPhys(E);
782 
783  // Forward transform to fill the coefficients space
784  for(int i = 0; i < m_fields.num_elements(); ++i)
785  {
786  m_fields[i]->GetBndCondExpansions()[bcRegion]->
787  FwdTrans_BndConstrained(
788  m_fields[i]->GetBndCondExpansions()[bcRegion]->
789  GetPhys(),
790  m_fields[i]->GetBndCondExpansions()[bcRegion]->
791  UpdateCoeffs());
792  }
793 
794  }
795 
796  // Calculation of the initial internal values as a weighted average
797  // over the distance from the boundaries
798  int nq = m_fields[0]->GetNpoints();
799 
800  Array<OneD,NekDouble> x0(nq);
801  Array<OneD,NekDouble> x1(nq);
802  Array<OneD,NekDouble> x2(nq);
803 
804  // Get the coordinates
805  // (assuming all fields have the same discretisation)
806  m_fields[0]->GetCoords(x0, x1, x2);
807 
808  for(int j = 0; j < nq; j++)
809  {
810  NekDouble Dist = 0.0;
811  NekDouble rho = 0.0;
812  NekDouble rhou = 0.0;
813  NekDouble rhov = 0.0;
814  NekDouble E = 0.0;
815  NekDouble SumDist = 0.0;
816 
817  // Calculation of all the distances
818  // Loop on all the edges of the input file
819  for (int bcRegion = 0; bcRegion < nbnd; ++bcRegion)
820  {
821  // Get quadrature points on the edge
822  int npoints = m_fields[0]->
823  GetBndCondExpansions()[bcRegion]->GetNpoints();
824 
825  Array<OneD,NekDouble> xb0(npoints, 0.0);
826  Array<OneD,NekDouble> xb1(npoints, 0.0);
827  Array<OneD,NekDouble> xb2(npoints, 0.0);
828 
829  m_fields[0]->GetBndCondExpansions()[bcRegion]->
830  GetCoords(xb0, xb1, xb2);
831 
832  for (int k = 0; k < npoints; k++)
833  {
834  Dist = sqrt((xb0[k] - x0[j]) * (xb0[k] - x0[j]) +
835  (xb1[k] - x1[j]) * (xb1[k] - x1[j]));
836 
837  SumDist += Dist;
838  rho += Dist * (m_fields[0]->
839  GetBndCondExpansions()[bcRegion]->GetPhys())[k];
840  rhou += Dist * (m_fields[1]->
841  GetBndCondExpansions()[bcRegion]->GetPhys())[k];
842  rhov += Dist * (m_fields[2]->
843  GetBndCondExpansions()[bcRegion]->GetPhys())[k];
844  E += Dist * (m_fields[3]->
845  GetBndCondExpansions()[bcRegion]->GetPhys())[k];
846  }
847  }
848 
849  rho = rho / SumDist;
850  rhou = rhou / SumDist;
851  rhov = rhov / SumDist;
852  E = E / SumDist;
853 
854  (m_fields[0]->UpdatePhys())[j] = rho;
855  (m_fields[1]->UpdatePhys())[j] = rhou;
856  (m_fields[2]->UpdatePhys())[j] = rhov;
857  (m_fields[3]->UpdatePhys())[j] = E;
858 
859  }
860 
861  for (int i = 0 ; i < m_fields.num_elements(); i++)
862  {
863  m_fields[i]->SetPhysState(true);
864  m_fields[i]->FwdTrans(m_fields[i]->GetPhys(),
865  m_fields[i]->UpdateCoeffs());
866  }
867 
868  // Dump initial conditions to file
869  std::string outname = m_sessionName + "_initialRingleb.chk";
870  WriteFld(outname);
871  }
872 
873  /**
874  * @brief Set the boundary conditions for the Ringleb flow problem.
875  */
876  void EulerCFE::SetBoundaryRinglebFlow(
877  int bcRegion,
878  NekDouble time,
879  int cnt,
880  Array<OneD, Array<OneD, NekDouble> > &physarray)
881  {
882  int nvariables = physarray.num_elements();
883  int nTraceNumPoints = GetTraceTotPoints();
884  Array<OneD, Array<OneD, NekDouble> > Fwd(nvariables);
885 
886  // For 3DHomogenoeus1D
887  int n_planes = 1;
888  if (m_expdim == 2 && m_HomogeneousType == eHomogeneous1D)
889  {
890  int nPointsTot = m_fields[0]->GetTotPoints();
891  int nPointsTot_plane = m_fields[0]->GetPlane(0)->GetTotPoints();
892  n_planes = nPointsTot/nPointsTot_plane;
893  nTraceNumPoints = nTraceNumPoints * n_planes;
894  }
895 
896  // Get physical values of the forward trace (from exp to phys)
897  for (int i = 0; i < nvariables; ++i)
898  {
899  Fwd[i] = Array<OneD, NekDouble>(nTraceNumPoints);
900  m_fields[i]->ExtractTracePhys(physarray[i], Fwd[i]);
901  }
902 
903  int id2, id2_plane, e_max;
904 
905  e_max = m_fields[0]->GetBndCondExpansions()[bcRegion]->GetExpSize();
906 
907  for(int e = 0; e < e_max; ++e)
908  {
909  int npoints = m_fields[0]->
910  GetBndCondExpansions()[bcRegion]->GetExp(e)->GetTotPoints();
911  int id1 = m_fields[0]->
912  GetBndCondExpansions()[bcRegion]->GetPhys_Offset(e);
913 
914  // For 3DHomogenoeus1D
915  if (m_expdim == 2 && m_HomogeneousType == eHomogeneous1D)
916  {
917  int cnt_plane = cnt/n_planes;
918  int e_plane;
919  int e_max_plane = e_max/n_planes;
920  int nTracePts_plane = GetTraceTotPoints();
921 
922  int planeID = floor((e + 0.5 )/ e_max_plane );
923  e_plane = e - e_max_plane*planeID;
924 
925  id2_plane = m_fields[0]->GetTrace()->GetPhys_Offset(
926  m_fields[0]->GetTraceMap()->
927  GetBndCondCoeffsToGlobalCoeffsMap(
928  cnt_plane + e_plane));
929  id2 = id2_plane + planeID*nTracePts_plane;
930  }
931  else // For general case
932  {
933  id2 = m_fields[0]->
934  GetTrace()->GetPhys_Offset(m_fields[0]->GetTraceMap()->
935  GetBndCondTraceToGlobalTraceMap(cnt++));
936  }
937 
938  Array<OneD,NekDouble> x0(npoints, 0.0);
939  Array<OneD,NekDouble> x1(npoints, 0.0);
940  Array<OneD,NekDouble> x2(npoints, 0.0);
941 
942  m_fields[0]->GetBndCondExpansions()[bcRegion]->
943  GetExp(e)->GetCoords(x0, x1, x2);
944 
945  // Flow parameters
946  NekDouble c, k, phi, r, J, VV, pp, sint, P, ss;
947  NekDouble J11, J12, J21, J22, det;
948  NekDouble Fx, Fy;
949  NekDouble xi, yi;
950  NekDouble dV;
951  NekDouble dtheta;
952  NekDouble par1;
953  NekDouble theta = M_PI / 4.0;
954  NekDouble kExt = 0.7;
955  NekDouble V = kExt * sin(theta);
956  NekDouble toll = 1.0e-8;
957  NekDouble errV = 1.0;
958  NekDouble errTheta = 1.0;
959  NekDouble gamma = m_gamma;
960  NekDouble gamma_1_2 = (gamma - 1.0) / 2.0;
961 
962  // Loop on all the points of that edge
963  for (int j = 0; j < npoints; j++)
964  {
965 
966  while ((abs(errV) > toll) || (abs(errTheta) > toll))
967  {
968  VV = V * V;
969  sint = sin(theta);
970  c = sqrt(1.0 - gamma_1_2 * VV);
971  k = V / sint;
972  phi = 1.0 / k;
973  pp = phi * phi;
974  J = 1.0 / c + 1.0 / (3.0 * c * c * c) +
975  1.0 / (5.0 * c * c * c * c * c) -
976  0.5 * log((1.0 + c) / (1.0 - c));
977 
978  r = pow(c, 1.0 / gamma_1_2);
979  xi = 1.0 / (2.0 * r) * (1.0 / VV - 2.0 * pp) + J / 2.0;
980  yi = phi / (r * V) * sqrt(1.0 - VV * pp);
981  par1 = 25.0 - 5.0 * VV;
982  ss = sint * sint;
983 
984  Fx = xi - x0[j];
985  Fy = yi - x1[j];
986 
987  J11 = 39062.5 / pow(par1, 3.5) *
988  (1.0 / VV - 2.0 / VV * ss) * V + 1562.5 /
989  pow(par1, 2.5) * (-2.0 / (VV * V) + 4.0 /
990  (VV * V) * ss) + 12.5 / pow(par1, 1.5) * V +
991  312.5 / pow(par1, 2.5) * V + 7812.5 /
992  pow(par1, 3.5) * V - 0.25 *
993  (-1.0 / pow(par1, 0.5) * V / (1.0 - 0.2 *
994  pow(par1, 0.5)) - (1.0 + 0.2 * pow(par1, 0.5)) /
995  pow((1.0 - 0.2 * pow(par1, 0.5)), 2.0) /
996  pow(par1, 0.5) * V) / (1.0 + 0.2 * pow(par1, 0.5)) *
997  (1.0 - 0.2 * pow(par1, 0.5));
998 
999  J12 = -6250.0 / pow(par1, 2.5) / VV * sint * cos(theta);
1000  J21 = -6250.0 / (VV * V) * sint / pow(par1, 2.5) *
1001  pow((1.0 - ss), 0.5) + 78125.0 / V * sint /
1002  pow(par1, 3.5) * pow((1.0 - ss), 0.5);
1003 
1004  // the matrix is singular when theta = pi/2
1005  if (abs(x1[j]) < toll && abs(cos(theta)) < toll)
1006  {
1007  J22 = -39062.5 / pow(par1, 3.5) / V + 3125 /
1008  pow(par1, 2.5) / (VV * V) + 12.5 /
1009  pow(par1, 1.5) * V + 312.5 / pow(par1, 2.5) *
1010  V + 7812.5 / pow(par1, 3.5) * V - 0.25 *
1011  (-1.0 / pow(par1, 0.5) * V / (1.0 - 0.2 *
1012  pow(par1, 0.5)) - (1.0 + 0.2 * pow(par1, 0.5)) /
1013  pow((1.0 - 0.2 * pow(par1, 0.5)), 2.0) /
1014  pow(par1, 0.5) * V) / (1.0 + 0.2 *
1015  pow(par1, 0.5)) * (1.0 - 0.2 * pow(par1, 0.5));
1016 
1017  // dV = -dV/dx * Fx
1018  dV = -1.0 / J22 * Fx;
1019  dtheta = 0.0;
1020  theta = M_PI / 2.0;
1021  }
1022  else
1023  {
1024  J22 = 3125.0 / VV * cos(theta) / pow(par1, 2.5) *
1025  pow((1.0 - ss), 0.5) - 3125.0 / VV * ss /
1026  pow(par1, 2.5) / pow((1.0 - ss), 0.5) *
1027  cos(theta);
1028 
1029  det = -1.0 / (J11 * J22 - J12 * J21);
1030 
1031  // [dV dtheta]' = -[invJ]*[Fx Fy]'
1032  dV = det * ( J22 * Fx - J12 * Fy);
1033  dtheta = det * (-J21 * Fx + J11 * Fy);
1034  }
1035 
1036  V = V + dV;
1037  theta = theta + dtheta;
1038 
1039  errV = abs(dV);
1040  errTheta = abs(dtheta);
1041  }
1042 
1043  c = sqrt(1.0 - gamma_1_2 * VV);
1044  int kk = id2 + j;
1045  NekDouble timeramp = 200.0;
1046  std::string restartstr = "RESTART";
1047  if (time<timeramp &&
1048  !(m_session->DefinesFunction("InitialConditions") &&
1049  m_session->GetFunctionType("InitialConditions", 0) ==
1050  LibUtilities::eFunctionTypeFile))
1051  {
1052  Fwd[0][kk] = pow(c, 1.0 / gamma_1_2) *
1053  exp(-1.0 + time /timeramp);
1054 
1055  Fwd[1][kk] = Fwd[0][kk] * V * cos(theta) *
1056  exp(-1 + time / timeramp);
1057 
1058  Fwd[2][kk] = Fwd[0][kk] * V * sin(theta) *
1059  exp(-1 + time / timeramp);
1060  }
1061  else
1062  {
1063  Fwd[0][kk] = pow(c, 1.0 / gamma_1_2);
1064  Fwd[1][kk] = Fwd[0][kk] * V * cos(theta);
1065  Fwd[2][kk] = Fwd[0][kk] * V * sin(theta);
1066  }
1067 
1068  P = (c * c) * Fwd[0][kk] / gamma;
1069  Fwd[3][kk] = P / (gamma - 1.0) + 0.5 *
1070  (Fwd[1][kk] * Fwd[1][kk] / Fwd[0][kk] +
1071  Fwd[2][kk] * Fwd[2][kk] / Fwd[0][kk]);
1072 
1073  errV = 1.0;
1074  errTheta = 1.0;
1075  theta = M_PI / 4.0;
1076  V = kExt * sin(theta);
1077  }
1078 
1079  for (int i = 0; i < nvariables; ++i)
1080  {
1081  Vmath::Vcopy(npoints, &Fwd[i][id2], 1,
1082  &(m_fields[i]->GetBndCondExpansions()[bcRegion]->
1083  UpdatePhys())[id1],1);
1084  }
1085  }
1086  }
1087 }
Generated on Wed Mar 4 2015 16:33:13 for Nektar++ by   doxygen 1.8.1.2