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, outarray);
170 
171  for (i = 0; i < nvariables; ++i)
172  {
173  Vmath::Neg(npoints, outarray[i], 1);
174  }
175  }
176 
177  /**
178  * @brief Compute the projection and call the method for imposing the
179  * boundary conditions in case of discontinuous projection.
180  */
181  void EulerCFE::DoOdeProjection(
182  const Array<OneD, const Array<OneD, NekDouble> > &inarray,
183  Array<OneD, Array<OneD, NekDouble> > &outarray,
184  const NekDouble time)
185  {
186  int i;
187  int nvariables = inarray.num_elements();
188 
189  switch (m_projectionType)
190  {
191  case MultiRegions::eDiscontinuous:
192  {
193  // Just copy over array
194  int npoints = GetNpoints();
195 
196  for(i = 0; i < nvariables; ++i)
197  {
198  Vmath::Vcopy(npoints, inarray[i], 1, outarray[i], 1);
199  }
200  SetBoundaryConditions(outarray, time);
201  break;
202  }
203  case MultiRegions::eGalerkin:
204  case MultiRegions::eMixed_CG_Discontinuous:
205  {
206  ASSERTL0(false, "No Continuous Galerkin for Euler equations");
207  break;
208  }
209  default:
210  ASSERTL0(false, "Unknown projection scheme");
211  break;
212  }
213  }
214 
215  /**
216  * @brief Set boundary conditions which can be:
217  * a) Wall and Symmerty BCs implemented at CompressibleFlowSystem level
218  since they are compressible solver specific;
219  * b) Isentropic vortex and Ringleb flow BCs implemented at EulerCFE level
220  * since they are Euler solver specific;
221  * c) Time dependent BCs.
222  *
223  * @param inarray: fields array.
224  * @param time: time.
225  */
226  void EulerCFE::SetBoundaryConditions(
227  Array<OneD, Array<OneD, NekDouble> > &inarray,
228  NekDouble time)
229  {
230  std::string varName;
231  int nvariables = m_fields.num_elements();
232  int cnt = 0;
233 
234  // Loop over Boundary Regions
235  for (int n = 0; n < m_fields[0]->GetBndConditions().num_elements(); ++n)
236  {
237  // Wall Boundary Condition
238  if (m_fields[0]->GetBndConditions()[n]->GetUserDefined() ==
239  SpatialDomains::eWall)
240  {
241  WallBC(n, cnt, inarray);
242  }
243 
244  // Wall Boundary Condition
245  if (m_fields[0]->GetBndConditions()[n]->GetUserDefined() ==
246  SpatialDomains::eWallViscous)
247  {
248  ASSERTL0(false, "WallViscous is a wrong bc for the "
249  "Euler equations");
250  }
251 
252  // Symmetric Boundary Condition
253  if (m_fields[0]->GetBndConditions()[n]->GetUserDefined() ==
254  SpatialDomains::eSymmetry)
255  {
256  SymmetryBC(n, cnt, inarray);
257  }
258 
259  // Riemann invariant characteristic Boundary Condition (CBC)
260  if (m_fields[0]->GetBndConditions()[n]->GetUserDefined() ==
261  SpatialDomains::eRiemannInvariant)
262  {
263  RiemannInvariantBC(n, cnt, inarray);
264  }
265 
266  // Extrapolation of the data at the boundaries
267  if (m_fields[0]->GetBndConditions()[n]->GetUserDefined() ==
268  SpatialDomains::eExtrapOrder0)
269  {
270  ExtrapOrder0BC(n, cnt, inarray);
271  }
272 
273  // Time Dependent Boundary Condition (specified in meshfile)
274  if (m_fields[0]->GetBndConditions()[n]->GetUserDefined() ==
275  SpatialDomains::eTimeDependent)
276  {
277  for (int i = 0; i < nvariables; ++i)
278  {
279  varName = m_session->GetVariable(i);
280  m_fields[i]->EvaluateBoundaryConditions(time, varName);
281  }
282  }
283 
284  // Isentropic Vortex Boundary Condition
285  if (m_fields[0]->GetBndConditions()[n]->GetUserDefined() ==
286  SpatialDomains::eIsentropicVortex)
287  {
288  SetBoundaryIsentropicVortex(n, time, cnt, inarray);
289  }
290 
291  // Ringleb Flow Inflow and outflow Boundary Condition
292  if (m_fields[0]->GetBndConditions()[n]->GetUserDefined() ==
293  SpatialDomains::eRinglebFlow)
294  {
295  SetBoundaryRinglebFlow(n, time, cnt, inarray);
296  }
297 
298  cnt += m_fields[0]->GetBndCondExpansions()[n]->GetExpSize();
299  }
300  }
301 
302  /**
303  * @brief Get the exact solutions for isentropic vortex and Ringleb
304  * flow problems.
305  */
306  void EulerCFE::v_EvaluateExactSolution(
307  unsigned int field,
308  Array<OneD, NekDouble> &outfield,
309  const NekDouble time)
310  {
311  switch(m_problemType)
312  {
313  case eIsentropicVortex:
314  {
315  GetExactIsentropicVortex(field, outfield, time);
316  break;
317  }
318  case eRinglebFlow:
319  {
320  GetExactRinglebFlow(field, outfield);
321  break;
322  }
323  default:
324  {
325  break;
326  }
327  }
328  }
329 
330  void EulerCFE::EvaluateIsentropicVortex(
331  const Array<OneD, NekDouble> &x,
332  const Array<OneD, NekDouble> &y,
333  const Array<OneD, NekDouble> &z,
334  Array<OneD, Array<OneD, NekDouble> > &u,
335  NekDouble time,
336  const int o)
337  {
338  int nq = x.num_elements();
339 
340  // Flow parameters
341  const NekDouble x0 = 5.0;
342  const NekDouble y0 = 0.0;
343  const NekDouble beta = 5.0;
344  const NekDouble u0 = 1.0;
345  const NekDouble v0 = 0.5;
346  const NekDouble gamma = m_gamma;
347  NekDouble r, xbar, ybar, tmp;
348  NekDouble fac = 1.0/(16.0*gamma*M_PI*M_PI);
349 
350  // In 3D zero rhow field.
351  if (m_spacedim == 3)
352  {
353  Vmath::Zero(nq, &u[3][o], 1);
354  }
355 
356  // Fill storage
357  for (int i = 0; i < nq; ++i)
358  {
359  xbar = x[i] - u0*time - x0;
360  ybar = y[i] - v0*time - y0;
361  r = sqrt(xbar*xbar + ybar*ybar);
362  tmp = beta*exp(1-r*r);
363  u[0][i+o] = pow(1.0 - (gamma-1.0)*tmp*tmp*fac, 1.0/(gamma-1.0));
364  u[1][i+o] = u[0][i+o]*(u0 - tmp*ybar/(2*M_PI));
365  u[2][i+o] = u[0][i+o]*(v0 + tmp*xbar/(2*M_PI));
366  u[m_spacedim+1][i+o] = pow(u[0][i+o], gamma)/(gamma-1.0) +
367  0.5*(u[1][i+o]*u[1][i+o] + u[2][i+o]*u[2][i+o]) / u[0][i+o];
368  }
369  }
370 
371  /**
372  * @brief Compute the exact solution for the isentropic vortex problem.
373  */
374  void EulerCFE::GetExactIsentropicVortex(
375  int field,
376  Array<OneD, NekDouble> &outarray,
377  NekDouble time)
378  {
379  int nTotQuadPoints = GetTotPoints();
380  Array<OneD, NekDouble> x(nTotQuadPoints);
381  Array<OneD, NekDouble> y(nTotQuadPoints);
382  Array<OneD, NekDouble> z(nTotQuadPoints);
383  Array<OneD, Array<OneD, NekDouble> > u(m_spacedim+2);
384 
385  m_fields[0]->GetCoords(x, y, z);
386 
387  for (int i = 0; i < m_spacedim + 2; ++i)
388  {
389  u[i] = Array<OneD, NekDouble>(nTotQuadPoints);
390  }
391 
392  EvaluateIsentropicVortex(x, y, z, u, time);
393 
394  Vmath::Vcopy(nTotQuadPoints, u[field], 1, outarray, 1);
395  }
396 
397  /**
398  * @brief Set the initial condition for the isentropic vortex problem.
399  */
400  void EulerCFE::SetInitialIsentropicVortex(NekDouble initialtime)
401  {
402  int nTotQuadPoints = GetTotPoints();
403  Array<OneD, NekDouble> x(nTotQuadPoints);
404  Array<OneD, NekDouble> y(nTotQuadPoints);
405  Array<OneD, NekDouble> z(nTotQuadPoints);
406  Array<OneD, Array<OneD, NekDouble> > u(m_spacedim+2);
407 
408  m_fields[0]->GetCoords(x, y, z);
409 
410  for (int i = 0; i < m_spacedim + 2; ++i)
411  {
412  u[i] = Array<OneD, NekDouble>(nTotQuadPoints);
413  }
414 
415  EvaluateIsentropicVortex(x, y, z, u, initialtime);
416 
417  // Forward transform to fill the coefficient space
418  for(int i = 0; i < m_fields.num_elements(); ++i)
419  {
420  Vmath::Vcopy(nTotQuadPoints, u[i], 1, m_fields[i]->UpdatePhys(), 1);
421  m_fields[i]->SetPhysState(true);
422  m_fields[i]->FwdTrans(m_fields[i]->GetPhys(),
423  m_fields[i]->UpdateCoeffs());
424  }
425  }
426 
427  /**
428  * @brief Set the boundary conditions for the isentropic vortex problem.
429  */
430  void EulerCFE::SetBoundaryIsentropicVortex(
431  int bcRegion,
432  NekDouble time,
433  int cnt,
434  Array<OneD, Array<OneD, NekDouble> > &physarray)
435  {
436  int nvariables = physarray.num_elements();
437  int nTraceNumPoints = GetTraceTotPoints();
438  Array<OneD, Array<OneD, NekDouble> > Fwd(nvariables);
439  const Array<OneD, const int> &bndTraceMap = m_fields[0]->GetTraceBndMap();
440 
441  // Get physical values of the forward trace (from exp to phys)
442  for (int i = 0; i < nvariables; ++i)
443  {
444  Fwd[i] = Array<OneD, NekDouble>(nTraceNumPoints);
445  m_fields[i]->ExtractTracePhys(physarray[i], Fwd[i]);
446  }
447 
448  int id2, e_max;
449  e_max = m_fields[0]->GetBndCondExpansions()[bcRegion]->GetExpSize();
450 
451  for(int e = 0; e < e_max; ++e)
452  {
453  int npoints = m_fields[0]->
454  GetBndCondExpansions()[bcRegion]->GetExp(e)->GetTotPoints();
455  int id1 = m_fields[0]->
456  GetBndCondExpansions()[bcRegion]->GetPhys_Offset(e);
457  id2 = m_fields[0]->GetTrace()->GetPhys_Offset(bndTraceMap[cnt++]);
458 
459  Array<OneD,NekDouble> x(npoints, 0.0);
460  Array<OneD,NekDouble> y(npoints, 0.0);
461  Array<OneD,NekDouble> z(npoints, 0.0);
462 
463  m_fields[0]->GetBndCondExpansions()[bcRegion]->
464  GetExp(e)->GetCoords(x, y, z);
465 
466  EvaluateIsentropicVortex(x, y, z, Fwd, time, id2);
467 
468  for (int i = 0; i < nvariables; ++i)
469  {
470  Vmath::Vcopy(npoints, &Fwd[i][id2], 1,
471  &(m_fields[i]->GetBndCondExpansions()[bcRegion]->
472  UpdatePhys())[id1], 1);
473  }
474  }
475  }
476 
477  /**
478  * @brief Compute the exact solution for the Ringleb flow problem.
479  */
480  void EulerCFE::GetExactRinglebFlow(
481  int field,
482  Array<OneD, NekDouble> &outarray)
483  {
484  int nTotQuadPoints = GetTotPoints();
485 
486  Array<OneD, NekDouble> rho(nTotQuadPoints, 100.0);
487  Array<OneD, NekDouble> rhou(nTotQuadPoints);
488  Array<OneD, NekDouble> rhov(nTotQuadPoints);
489  Array<OneD, NekDouble> E(nTotQuadPoints);
490  Array<OneD, NekDouble> x(nTotQuadPoints);
491  Array<OneD, NekDouble> y(nTotQuadPoints);
492  Array<OneD, NekDouble> z(nTotQuadPoints);
493 
494  m_fields[0]->GetCoords(x, y, z);
495 
496  // Flow parameters
497  NekDouble c, k, phi, r, J, VV, pp, sint, P, ss;
498  NekDouble J11, J12, J21, J22, det;
499  NekDouble Fx, Fy;
500  NekDouble xi, yi;
501  NekDouble dV;
502  NekDouble dtheta;
503  NekDouble par1;
504  NekDouble theta = M_PI / 4.0;
505  NekDouble kExt = 0.7;
506  NekDouble V = kExt * sin(theta);
507  NekDouble toll = 1.0e-8;
508  NekDouble errV = 1.0;
509  NekDouble errTheta = 1.0;
510  NekDouble gamma = m_gamma;
511  NekDouble gamma_1_2 = (gamma - 1.0) / 2.0;
512 
513  for (int i = 0; i < nTotQuadPoints; ++i)
514  {
515  while ((abs(errV) > toll) || (abs(errTheta) > toll))
516  {
517  VV = V * V;
518  sint = sin(theta);
519  c = sqrt(1.0 - gamma_1_2 * VV);
520  k = V / sint;
521  phi = 1.0 / k;
522  pp = phi * phi;
523  J = 1.0 / c + 1.0 / (3.0 * c * c * c) +
524  1.0 / (5.0 * c * c * c * c * c) -
525  0.5 * log((1.0 + c) / (1.0 - c));
526 
527  r = pow(c, 1.0 / gamma_1_2);
528  xi = 1.0 / (2.0 * r) * (1.0 / VV - 2.0 * pp) + J / 2.0;
529  yi = phi / (r * V) * sqrt(1.0 - VV * pp);
530  par1 = 25.0 - 5.0 * VV;
531  ss = sint * sint;
532 
533  Fx = xi - x[i];
534  Fy = yi - y[i];
535 
536  J11 = 39062.5 / pow(par1, 3.5) * (1.0 / VV - 2.0 / VV * ss) *
537  V + 1562.5 / pow(par1, 2.5) * (-2.0 / (VV * V) + 4.0 /
538  (VV * V) * ss) + 12.5 / pow(par1, 1.5) * V + 312.5 /
539  pow(par1, 2.5) * V + 7812.5 / pow(par1, 3.5) * V -
540  0.25 * (-1.0 / pow(par1, 0.5) * V/(1.0 - 0.2 *
541  pow(par1, 0.5)) - (1.0 + 0.2 * pow(par1, 0.5)) /
542  pow((1.0 - 0.2 * pow(par1, 0.5)), 2.0) /
543  pow(par1, 0.5) * V) / (1.0 + 0.2 * pow(par1, 0.5)) *
544  (1.0 - 0.2 * pow(par1, 0.5));
545 
546  J12 = -6250.0 / pow(par1, 2.5) / VV * sint * cos(theta);
547  J21 = -6250.0 / (VV * V) * sint /
548  pow(par1, 2.5) * pow((1.0 - ss), 0.5) +
549  78125.0 / V * sint / pow(par1, 3.5) *
550  pow((1.0 - ss), 0.5);
551 
552  // the matrix is singular when theta = pi/2
553  if(abs(y[i])<toll && abs(cos(theta))<toll)
554  {
555  J22 = -39062.5 / pow(par1, 3.5) / V + 3125 /
556  pow(par1, 2.5) / (VV * V) + 12.5 / pow(par1, 1.5) *
557  V + 312.5 / pow(par1, 2.5) * V + 7812.5 /
558  pow(par1, 3.5) * V - 0.25 * (-1.0 / pow(par1, 0.5) *
559  V / (1.0 - 0.2 * pow(par1, 0.5)) - (1.0 + 0.2 *
560  pow(par1, 0.5)) / pow((1.0 - 0.2 *
561  pow(par1, 0.5)), 2.0) / pow(par1, 0.5) * V) /
562  (1.0 + 0.2 * pow(par1, 0.5)) * (1.0 - 0.2 *
563  pow(par1,0.5));
564 
565  // dV = -dV/dx * Fx
566  dV = -1.0 / J22 * Fx;
567  dtheta = 0.0;
568  theta = M_PI / 2.0;
569  }
570  else
571  {
572  J22 = 3125.0 / VV * cos(theta) / pow(par1, 2.5) *
573  pow((1.0 - ss), 0.5) - 3125.0 / VV * ss /
574  pow(par1, 2.5) / pow((1.0 - ss), 0.5) * cos(theta);
575 
576  det = -1.0 / (J11 * J22 - J12 * J21);
577 
578  // [dV dtheta]' = -[invJ]*[Fx Fy]'
579  dV = det * ( J22 * Fx - J12 * Fy);
580  dtheta = det * (-J21 * Fx + J11 * Fy);
581  }
582 
583  V = V + dV;
584  theta = theta + dtheta;
585 
586  errV = abs(dV);
587  errTheta = abs(dtheta);
588 
589  }
590 
591  c = sqrt(1.0 - gamma_1_2 * VV);
592  r = pow(c, 1.0 / gamma_1_2);
593 
594  rho[i] = r;
595  rhou[i] = rho[i] * V * cos(theta);
596  rhov[i] = rho[i] * V * sin(theta);
597  P = (c * c) * rho[i] / gamma;
598  E[i] = P / (gamma - 1.0) + 0.5 *
599  (rhou[i] * rhou[i] / rho[i] + rhov[i] * rhov[i] / rho[i]);
600 
601  // Resetting the guess value
602  errV = 1.0;
603  errTheta = 1.0;
604  theta = M_PI/4.0;
605  V = kExt*sin(theta);
606  }
607 
608  switch (field)
609  {
610  case 0:
611  outarray = rho;
612  break;
613  case 1:
614  outarray = rhou;
615  break;
616  case 2:
617  outarray = rhov;
618  break;
619  case 3:
620  outarray = E;
621  break;
622  default:
623  ASSERTL0(false, "Error in variable number!");
624  break;
625  }
626  }
627 
628  /**
629  * @brief Set the initial condition for the Ringleb flow problem.
630  */
631  void EulerCFE::SetInitialRinglebFlow(void)
632  {
633  // Get number of different boundaries in the input file
634  int nbnd = m_fields[0]->GetBndConditions().num_elements();
635 
636  // Loop on all the edges of the input file
637  for(int bcRegion=0; bcRegion < nbnd; ++bcRegion)
638  {
639 
640  int npoints = m_fields[0]->
641  GetBndCondExpansions()[bcRegion]->GetNpoints();
642 
643  Array<OneD,NekDouble> x0(npoints, 0.0);
644  Array<OneD,NekDouble> x1(npoints, 0.0);
645  Array<OneD,NekDouble> x2(npoints, 0.0);
646 
647  Array<OneD, NekDouble> rho(npoints, 0.0);
648  Array<OneD, NekDouble> rhou(npoints, 0.0);
649  Array<OneD, NekDouble> rhov(npoints, 0.0);
650  Array<OneD, NekDouble> E(npoints, 0.0);
651 
652  m_fields[0]->GetBndCondExpansions()[bcRegion]->
653  GetCoords(x0, x1, x2);
654 
655  // Flow parameters
656  NekDouble c, k, phi, r, J, VV, pp, sint, P, ss;
657  NekDouble J11, J12, J21, J22, det;
658  NekDouble Fx, Fy;
659  NekDouble xi, yi;
660  NekDouble dV;
661  NekDouble dtheta;
662  NekDouble par1;
663  NekDouble theta = M_PI / 4.0;
664  NekDouble kExt = 0.7;
665  NekDouble V = kExt * sin(theta);
666  NekDouble toll = 1.0e-8;
667  NekDouble errV = 1.0;
668  NekDouble errTheta = 1.0;
669  NekDouble gamma = m_gamma;
670  NekDouble gamma_1_2 = (gamma - 1.0) / 2.0;
671 
672  // Loop on all the points of that edge
673  for (int j = 0; j < npoints; j++)
674  {
675  while ((abs(errV) > toll) || (abs(errTheta) > toll))
676  {
677 
678  VV = V * V;
679  sint = sin(theta);
680  c = sqrt(1.0 - gamma_1_2 * VV);
681  k = V / sint;
682  phi = 1.0 / k;
683  pp = phi * phi;
684  J = 1.0 / c + 1.0 / (3.0 * c * c * c) +
685  1.0 / (5.0 * c * c * c * c * c) -
686  0.5 * log((1.0 + c) / (1.0 - c));
687 
688  r = pow(c, 1.0 / gamma_1_2);
689  xi = 1.0 / (2.0 * r) * (1.0 / VV - 2.0 * pp) + J / 2.0;
690  yi = phi / (r * V) * sqrt(1.0 - VV * pp);
691  par1 = 25.0 - 5.0 * VV;
692  ss = sint * sint;
693 
694  Fx = xi - x0[j];
695  Fy = yi - x1[j];
696 
697  J11 = 39062.5 / pow(par1, 3.5) *
698  (1.0 / VV - 2.0 / VV * ss) * V +
699  1562.5 / pow(par1, 2.5) * (-2.0 /
700  (VV * V) + 4.0 / (VV * V) * ss) +
701  12.5 / pow(par1, 1.5) * V +
702  312.5 / pow(par1, 2.5) * V +
703  7812.5 / pow(par1, 3.5) * V -
704  0.25 * (-1.0 / pow(par1, 0.5) * V /
705  (1.0 - 0.2 * pow(par1, 0.5)) - (1.0 + 0.2 *
706  pow(par1, 0.5)) / pow((1.0 - 0.2 *
707  pow(par1, 0.5)), 2.0) /
708  pow(par1, 0.5) * V) /
709  (1.0 + 0.2 * pow(par1, 0.5)) *
710  (1.0 - 0.2 * pow(par1, 0.5));
711 
712  J12 = -6250.0 / pow(par1, 2.5) / VV * sint * cos(theta);
713  J21 = -6250.0 / (VV * V) * sint / pow(par1, 2.5) *
714  pow((1.0 - ss), 0.5) + 78125.0 / V * sint /
715  pow(par1, 3.5) * pow((1.0 - ss), 0.5);
716 
717  // the matrix is singular when theta = pi/2
718  if (abs(x1[j]) < toll && abs(cos(theta)) < toll)
719  {
720 
721  J22 = -39062.5 / pow(par1, 3.5) / V +
722  3125 / pow(par1, 2.5) / (VV * V) + 12.5 /
723  pow(par1, 1.5) * V + 312.5 / pow(par1, 2.5) * V +
724  7812.5 / pow(par1, 3.5) * V -
725  0.25 * (-1.0 / pow(par1, 0.5) * V /
726  (1.0 - 0.2 * pow(par1, 0.5)) -
727  (1.0 + 0.2 * pow(par1, 0.5)) /
728  pow((1.0 - 0.2* pow(par1, 0.5)), 2.0) /
729  pow(par1, 0.5) *V) /
730  (1.0 + 0.2 * pow(par1, 0.5)) *
731  (1.0 - 0.2 * pow(par1, 0.5));
732 
733  // dV = -dV/dx * Fx
734  dV = -1.0 / J22 * Fx;
735  dtheta = 0.0;
736  theta = M_PI / 2.0;
737  }
738  else
739  {
740 
741  J22 = 3125.0 / VV * cos(theta) / pow(par1, 2.5) *
742  pow((1.0 - ss), 0.5) - 3125.0 / VV * ss /
743  pow(par1, 2.5) / pow((1.0 - ss), 0.5) *
744  cos(theta);
745 
746  det = -1.0 / (J11 * J22 - J12 * J21);
747 
748  // [dV dtheta]' = -[invJ]*[Fx Fy]'
749  dV = det * ( J22 * Fx - J12 * Fy);
750  dtheta = det * (-J21 * Fx + J11 * Fy);
751  }
752 
753  V = V + dV;
754  theta = theta + dtheta;
755 
756  errV = abs(dV);
757  errTheta = abs(dtheta);
758 
759  }
760 
761  c = sqrt(1.0 - gamma_1_2 * VV);
762  rho[j] = pow(c, 1.0 / gamma_1_2) * exp(-1.0);
763  rhou[j] = rho[j] * V * cos(theta) * exp(-1.0);
764  rhov[j] = rho[j] * V * sin(theta) * exp(-1.0);
765  P = (c * c) * rho[j] / gamma;
766  E[j] = P / (gamma - 1.0) + 0.5 *
767  (rhou[j] * rhou[j] /
768  rho[j] + rhov[j] * rhov[j] / rho[j]);
769  errV = 1.0;
770  errTheta = 1.0;
771  theta = M_PI / 4.0;
772  V = kExt * sin(theta);
773 
774  }
775 
776  // Fill the physical space
777  m_fields[0]->GetBndCondExpansions()[bcRegion]->SetPhys(rho);
778  m_fields[1]->GetBndCondExpansions()[bcRegion]->SetPhys(rhou);
779  m_fields[2]->GetBndCondExpansions()[bcRegion]->SetPhys(rhov);
780  m_fields[3]->GetBndCondExpansions()[bcRegion]->SetPhys(E);
781 
782  // Forward transform to fill the coefficients space
783  for(int i = 0; i < m_fields.num_elements(); ++i)
784  {
785  m_fields[i]->GetBndCondExpansions()[bcRegion]->
786  FwdTrans_BndConstrained(
787  m_fields[i]->GetBndCondExpansions()[bcRegion]->
788  GetPhys(),
789  m_fields[i]->GetBndCondExpansions()[bcRegion]->
790  UpdateCoeffs());
791  }
792 
793  }
794 
795  // Calculation of the initial internal values as a weighted average
796  // over the distance from the boundaries
797  int nq = m_fields[0]->GetNpoints();
798 
799  Array<OneD,NekDouble> x0(nq);
800  Array<OneD,NekDouble> x1(nq);
801  Array<OneD,NekDouble> x2(nq);
802 
803  // Get the coordinates
804  // (assuming all fields have the same discretisation)
805  m_fields[0]->GetCoords(x0, x1, x2);
806 
807  for(int j = 0; j < nq; j++)
808  {
809  NekDouble Dist = 0.0;
810  NekDouble rho = 0.0;
811  NekDouble rhou = 0.0;
812  NekDouble rhov = 0.0;
813  NekDouble E = 0.0;
814  NekDouble SumDist = 0.0;
815 
816  // Calculation of all the distances
817  // Loop on all the edges of the input file
818  for (int bcRegion = 0; bcRegion < nbnd; ++bcRegion)
819  {
820  // Get quadrature points on the edge
821  int npoints = m_fields[0]->
822  GetBndCondExpansions()[bcRegion]->GetNpoints();
823 
824  Array<OneD,NekDouble> xb0(npoints, 0.0);
825  Array<OneD,NekDouble> xb1(npoints, 0.0);
826  Array<OneD,NekDouble> xb2(npoints, 0.0);
827 
828  m_fields[0]->GetBndCondExpansions()[bcRegion]->
829  GetCoords(xb0, xb1, xb2);
830 
831  for (int k = 0; k < npoints; k++)
832  {
833  Dist = sqrt((xb0[k] - x0[j]) * (xb0[k] - x0[j]) +
834  (xb1[k] - x1[j]) * (xb1[k] - x1[j]));
835 
836  SumDist += Dist;
837  rho += Dist * (m_fields[0]->
838  GetBndCondExpansions()[bcRegion]->GetPhys())[k];
839  rhou += Dist * (m_fields[1]->
840  GetBndCondExpansions()[bcRegion]->GetPhys())[k];
841  rhov += Dist * (m_fields[2]->
842  GetBndCondExpansions()[bcRegion]->GetPhys())[k];
843  E += Dist * (m_fields[3]->
844  GetBndCondExpansions()[bcRegion]->GetPhys())[k];
845  }
846  }
847 
848  rho = rho / SumDist;
849  rhou = rhou / SumDist;
850  rhov = rhov / SumDist;
851  E = E / SumDist;
852 
853  (m_fields[0]->UpdatePhys())[j] = rho;
854  (m_fields[1]->UpdatePhys())[j] = rhou;
855  (m_fields[2]->UpdatePhys())[j] = rhov;
856  (m_fields[3]->UpdatePhys())[j] = E;
857 
858  }
859 
860  for (int i = 0 ; i < m_fields.num_elements(); i++)
861  {
862  m_fields[i]->SetPhysState(true);
863  m_fields[i]->FwdTrans(m_fields[i]->GetPhys(),
864  m_fields[i]->UpdateCoeffs());
865  }
866 
867  // Dump initial conditions to file
868  std::string outname = m_sessionName + "_initialRingleb.chk";
869  WriteFld(outname);
870  }
871 
872  /**
873  * @brief Set the boundary conditions for the Ringleb flow problem.
874  */
875  void EulerCFE::SetBoundaryRinglebFlow(
876  int bcRegion,
877  NekDouble time,
878  int cnt,
879  Array<OneD, Array<OneD, NekDouble> > &physarray)
880  {
881  int nvariables = physarray.num_elements();
882  int nTraceNumPoints = GetTraceTotPoints();
883  Array<OneD, Array<OneD, NekDouble> > Fwd(nvariables);
884 
885  // For 3DHomogenoeus1D
886  int n_planes;
887  if (m_expdim == 2 && m_HomogeneousType == eHomogeneous1D)
888  {
889  int nPointsTot = m_fields[0]->GetTotPoints();
890  int nPointsTot_plane = m_fields[0]->GetPlane(0)->GetTotPoints();
891  n_planes = nPointsTot/nPointsTot_plane;
892  nTraceNumPoints = nTraceNumPoints * n_planes;
893  }
894 
895  // Get physical values of the forward trace (from exp to phys)
896  for (int i = 0; i < nvariables; ++i)
897  {
898  Fwd[i] = Array<OneD, NekDouble>(nTraceNumPoints);
899  m_fields[i]->ExtractTracePhys(physarray[i], Fwd[i]);
900  }
901 
902  int id2, id2_plane, e_max;
903 
904  e_max = m_fields[0]->GetBndCondExpansions()[bcRegion]->GetExpSize();
905 
906  for(int e = 0; e < e_max; ++e)
907  {
908  int npoints = m_fields[0]->
909  GetBndCondExpansions()[bcRegion]->GetExp(e)->GetTotPoints();
910  int id1 = m_fields[0]->
911  GetBndCondExpansions()[bcRegion]->GetPhys_Offset(e);
912 
913  // For 3DHomogenoeus1D
914  if (m_expdim == 2 && m_HomogeneousType == eHomogeneous1D)
915  {
916  int cnt_plane = cnt/n_planes;
917  int e_plane;
918  int e_max_plane = e_max/n_planes;
919  int nTracePts_plane = GetTraceTotPoints();
920 
921  int planeID = floor((e + 0.5 )/ e_max_plane );
922  e_plane = e - e_max_plane*planeID;
923 
924  id2_plane = m_fields[0]->GetTrace()->GetPhys_Offset(
925  m_fields[0]->GetTraceMap()->
926  GetBndCondCoeffsToGlobalCoeffsMap(
927  cnt_plane + e_plane));
928  id2 = id2_plane + planeID*nTracePts_plane;
929  }
930  else // For general case
931  {
932  id2 = m_fields[0]->
933  GetTrace()->GetPhys_Offset(m_fields[0]->GetTraceMap()->
934  GetBndCondTraceToGlobalTraceMap(cnt++));
935  }
936 
937  Array<OneD,NekDouble> x0(npoints, 0.0);
938  Array<OneD,NekDouble> x1(npoints, 0.0);
939  Array<OneD,NekDouble> x2(npoints, 0.0);
940 
941  m_fields[0]->GetBndCondExpansions()[bcRegion]->
942  GetExp(e)->GetCoords(x0, x1, x2);
943 
944  // Flow parameters
945  NekDouble c, k, phi, r, J, VV, pp, sint, P, ss;
946  NekDouble J11, J12, J21, J22, det;
947  NekDouble Fx, Fy;
948  NekDouble xi, yi;
949  NekDouble dV;
950  NekDouble dtheta;
951  NekDouble par1;
952  NekDouble theta = M_PI / 4.0;
953  NekDouble kExt = 0.7;
954  NekDouble V = kExt * sin(theta);
955  NekDouble toll = 1.0e-8;
956  NekDouble errV = 1.0;
957  NekDouble errTheta = 1.0;
958  NekDouble gamma = m_gamma;
959  NekDouble gamma_1_2 = (gamma - 1.0) / 2.0;
960 
961  // Loop on all the points of that edge
962  for (int j = 0; j < npoints; j++)
963  {
964 
965  while ((abs(errV) > toll) || (abs(errTheta) > toll))
966  {
967  VV = V * V;
968  sint = sin(theta);
969  c = sqrt(1.0 - gamma_1_2 * VV);
970  k = V / sint;
971  phi = 1.0 / k;
972  pp = phi * phi;
973  J = 1.0 / c + 1.0 / (3.0 * c * c * c) +
974  1.0 / (5.0 * c * c * c * c * c) -
975  0.5 * log((1.0 + c) / (1.0 - c));
976 
977  r = pow(c, 1.0 / gamma_1_2);
978  xi = 1.0 / (2.0 * r) * (1.0 / VV - 2.0 * pp) + J / 2.0;
979  yi = phi / (r * V) * sqrt(1.0 - VV * pp);
980  par1 = 25.0 - 5.0 * VV;
981  ss = sint * sint;
982 
983  Fx = xi - x0[j];
984  Fy = yi - x1[j];
985 
986  J11 = 39062.5 / pow(par1, 3.5) *
987  (1.0 / VV - 2.0 / VV * ss) * V + 1562.5 /
988  pow(par1, 2.5) * (-2.0 / (VV * V) + 4.0 /
989  (VV * V) * ss) + 12.5 / pow(par1, 1.5) * V +
990  312.5 / pow(par1, 2.5) * V + 7812.5 /
991  pow(par1, 3.5) * V - 0.25 *
992  (-1.0 / pow(par1, 0.5) * V / (1.0 - 0.2 *
993  pow(par1, 0.5)) - (1.0 + 0.2 * pow(par1, 0.5)) /
994  pow((1.0 - 0.2 * pow(par1, 0.5)), 2.0) /
995  pow(par1, 0.5) * V) / (1.0 + 0.2 * pow(par1, 0.5)) *
996  (1.0 - 0.2 * pow(par1, 0.5));
997 
998  J12 = -6250.0 / pow(par1, 2.5) / VV * sint * cos(theta);
999  J21 = -6250.0 / (VV * V) * sint / pow(par1, 2.5) *
1000  pow((1.0 - ss), 0.5) + 78125.0 / V * sint /
1001  pow(par1, 3.5) * pow((1.0 - ss), 0.5);
1002 
1003  // the matrix is singular when theta = pi/2
1004  if (abs(x1[j]) < toll && abs(cos(theta)) < toll)
1005  {
1006  J22 = -39062.5 / pow(par1, 3.5) / V + 3125 /
1007  pow(par1, 2.5) / (VV * V) + 12.5 /
1008  pow(par1, 1.5) * V + 312.5 / pow(par1, 2.5) *
1009  V + 7812.5 / pow(par1, 3.5) * V - 0.25 *
1010  (-1.0 / pow(par1, 0.5) * V / (1.0 - 0.2 *
1011  pow(par1, 0.5)) - (1.0 + 0.2 * pow(par1, 0.5)) /
1012  pow((1.0 - 0.2 * pow(par1, 0.5)), 2.0) /
1013  pow(par1, 0.5) * V) / (1.0 + 0.2 *
1014  pow(par1, 0.5)) * (1.0 - 0.2 * pow(par1, 0.5));
1015 
1016  // dV = -dV/dx * Fx
1017  dV = -1.0 / J22 * Fx;
1018  dtheta = 0.0;
1019  theta = M_PI / 2.0;
1020  }
1021  else
1022  {
1023  J22 = 3125.0 / VV * cos(theta) / pow(par1, 2.5) *
1024  pow((1.0 - ss), 0.5) - 3125.0 / VV * ss /
1025  pow(par1, 2.5) / pow((1.0 - ss), 0.5) *
1026  cos(theta);
1027 
1028  det = -1.0 / (J11 * J22 - J12 * J21);
1029 
1030  // [dV dtheta]' = -[invJ]*[Fx Fy]'
1031  dV = det * ( J22 * Fx - J12 * Fy);
1032  dtheta = det * (-J21 * Fx + J11 * Fy);
1033  }
1034 
1035  V = V + dV;
1036  theta = theta + dtheta;
1037 
1038  errV = abs(dV);
1039  errTheta = abs(dtheta);
1040  }
1041 
1042  c = sqrt(1.0 - gamma_1_2 * VV);
1043  int kk = id2 + j;
1044  NekDouble timeramp = 200.0;
1045  std::string restartstr = "RESTART";
1046  if (time<timeramp &&
1047  !(m_session->DefinesFunction("InitialConditions") &&
1048  m_session->GetFunctionType("InitialConditions", 0) ==
1049  LibUtilities::eFunctionTypeFile))
1050  {
1051  Fwd[0][kk] = pow(c, 1.0 / gamma_1_2) *
1052  exp(-1.0 + time /timeramp);
1053 
1054  Fwd[1][kk] = Fwd[0][kk] * V * cos(theta) *
1055  exp(-1 + time / timeramp);
1056 
1057  Fwd[2][kk] = Fwd[0][kk] * V * sin(theta) *
1058  exp(-1 + time / timeramp);
1059  }
1060  else
1061  {
1062  Fwd[0][kk] = pow(c, 1.0 / gamma_1_2);
1063  Fwd[1][kk] = Fwd[0][kk] * V * cos(theta);
1064  Fwd[2][kk] = Fwd[0][kk] * V * sin(theta);
1065  }
1066 
1067  P = (c * c) * Fwd[0][kk] / gamma;
1068  Fwd[3][kk] = P / (gamma - 1.0) + 0.5 *
1069  (Fwd[1][kk] * Fwd[1][kk] / Fwd[0][kk] +
1070  Fwd[2][kk] * Fwd[2][kk] / Fwd[0][kk]);
1071 
1072  errV = 1.0;
1073  errTheta = 1.0;
1074  theta = M_PI / 4.0;
1075  V = kExt * sin(theta);
1076  }
1077 
1078  for (int i = 0; i < nvariables; ++i)
1079  {
1080  Vmath::Vcopy(npoints, &Fwd[i][id2], 1,
1081  &(m_fields[i]->GetBndCondExpansions()[bcRegion]->
1082  UpdatePhys())[id1],1);
1083  }
1084  }
1085  }
1086 }
Generated on Sun Mar 15 2015 00:11:31 for Nektar++ by   doxygen 1.8.1.2