Nektar++
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
RinglebFlowBC.cpp
Go to the documentation of this file.
1 ///////////////////////////////////////////////////////////////////////////////
2 //
3 // File: RinglebFlowBC.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: Ringleb flow boundary condition
33 //
34 ///////////////////////////////////////////////////////////////////////////////
35 
36 #include "RinglebFlowBC.h"
38 
39 using namespace std;
40 
41 namespace Nektar
42 {
43 
44 std::string RinglebFlowBC::className = GetCFSBndCondFactory().
45  RegisterCreatorFunction("RinglebFlow",
46  RinglebFlowBC::create,
47  "Ringleb flow boundary condition.");
48 
49 RinglebFlowBC::RinglebFlowBC(const LibUtilities::SessionReaderSharedPtr& pSession,
51  const Array<OneD, Array<OneD, NekDouble> >& pTraceNormals,
52  const int pSpaceDim,
53  const int bcRegion,
54  const int cnt)
55  : CFSBndCond(pSession, pFields, pTraceNormals, pSpaceDim, bcRegion, cnt)
56 {
57  m_expdim = pFields[0]->GetGraph()->GetMeshDimension();
58 
59  m_homo1D = false;
60  if (m_session->DefinesSolverInfo("HOMOGENEOUS"))
61  {
62  std::string HomoStr = m_session->GetSolverInfo("HOMOGENEOUS");
63  if ((HomoStr == "HOMOGENEOUS1D") || (HomoStr == "Homogeneous1D")
64  || (HomoStr == "1D") || (HomoStr == "Homo1D"))
65  {
66  m_homo1D = true;
67  }
68  }
69 }
70 
73  Array<OneD, Array<OneD, NekDouble> > &physarray,
74  const NekDouble &time)
75 {
76  int nvariables = physarray.num_elements();
77 
78  // For 3DHomogenoeus1D
79  int n_planes = 1;
80  if (m_expdim == 2 && m_homo1D)
81  {
82  int nPointsTot = m_fields[0]->GetTotPoints();
83  int nPointsTot_plane = m_fields[0]->GetPlane(0)->GetTotPoints();
84  n_planes = nPointsTot/nPointsTot_plane;
85  }
86 
87  int id2, id2_plane, e_max;
88 
89  e_max = m_fields[0]->GetBndCondExpansions()[m_bcRegion]->GetExpSize();
90 
91  for(int e = 0; e < e_max; ++e)
92  {
93  int npoints = m_fields[0]->
94  GetBndCondExpansions()[m_bcRegion]->GetExp(e)->GetTotPoints();
95  int id1 = m_fields[0]->
96  GetBndCondExpansions()[m_bcRegion]->GetPhys_Offset(e);
97 
98  // For 3DHomogenoeus1D
99  if (m_expdim == 2 && m_homo1D)
100  {
101  int m_offset_plane = m_offset/n_planes;
102  int e_plane;
103  int e_max_plane = e_max/n_planes;
104  int nTracePts_plane = m_fields[0]->GetTrace()->GetNpoints();
105 
106  int planeID = floor((e + 0.5 )/ e_max_plane );
107  e_plane = e - e_max_plane*planeID;
108 
109  id2_plane = m_fields[0]->GetTrace()->GetPhys_Offset(
110  m_fields[0]->GetTraceMap()->
111  GetBndCondCoeffsToGlobalCoeffsMap(
112  m_offset_plane + e_plane));
113  id2 = id2_plane + planeID*nTracePts_plane;
114  }
115  else // For general case
116  {
117  id2 = m_fields[0]->
118  GetTrace()->GetPhys_Offset(m_fields[0]->GetTraceMap()->
119  GetBndCondTraceToGlobalTraceMap(m_offset+e));
120  }
121 
122  Array<OneD,NekDouble> x0(npoints, 0.0);
123  Array<OneD,NekDouble> x1(npoints, 0.0);
124  Array<OneD,NekDouble> x2(npoints, 0.0);
125 
126  m_fields[0]->GetBndCondExpansions()[m_bcRegion]->
127  GetExp(e)->GetCoords(x0, x1, x2);
128 
129  // Flow parameters
130  NekDouble c, k, phi, r, J, VV, pp, sint, P, ss;
131  NekDouble J11, J12, J21, J22, det;
132  NekDouble Fx, Fy;
133  NekDouble xi, yi;
134  NekDouble dV;
135  NekDouble dtheta;
136  NekDouble par1;
137  NekDouble theta = M_PI / 4.0;
138  NekDouble kExt = 0.7;
139  NekDouble V = kExt * sin(theta);
140  NekDouble toll = 1.0e-8;
141  NekDouble errV = 1.0;
142  NekDouble errTheta = 1.0;
143  NekDouble gamma = m_gamma;
144  NekDouble gamma_1_2 = (gamma - 1.0) / 2.0;
145 
146  // Loop on all the points of that edge
147  for (int j = 0; j < npoints; j++)
148  {
149 
150  while ((abs(errV) > toll) || (abs(errTheta) > toll))
151  {
152  VV = V * V;
153  sint = sin(theta);
154  c = sqrt(1.0 - gamma_1_2 * VV);
155  k = V / sint;
156  phi = 1.0 / k;
157  pp = phi * phi;
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));
161 
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;
166  ss = sint * sint;
167 
168  Fx = xi - x0[j];
169  Fy = yi - x1[j];
170 
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));
182 
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);
187 
188  // the matrix is singular when theta = pi/2
189  if (abs(x1[j]) < toll && abs(cos(theta)) < toll)
190  {
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));
200 
201  // dV = -dV/dx * Fx
202  dV = -1.0 / J22 * Fx;
203  dtheta = 0.0;
204  theta = M_PI / 2.0;
205  }
206  else
207  {
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) *
211  cos(theta);
212 
213  det = -1.0 / (J11 * J22 - J12 * J21);
214 
215  // [dV dtheta]' = -[invJ]*[Fx Fy]'
216  dV = det * ( J22 * Fx - J12 * Fy);
217  dtheta = det * (-J21 * Fx + J11 * Fy);
218  }
219 
220  V = V + dV;
221  theta = theta + dtheta;
222 
223  errV = abs(dV);
224  errTheta = abs(dtheta);
225  }
226 
227  c = sqrt(1.0 - gamma_1_2 * VV);
228  int kk = id2 + j;
229  NekDouble timeramp = 200.0;;
230  if (time<timeramp &&
231  !(m_session->DefinesFunction("InitialConditions") &&
232  m_session->GetFunctionType("InitialConditions", 0) ==
234  {
235  Fwd[0][kk] = pow(c, 1.0 / gamma_1_2) *
236  exp(-1.0 + time /timeramp);
237 
238  Fwd[1][kk] = Fwd[0][kk] * V * cos(theta) *
239  exp(-1 + time / timeramp);
240 
241  Fwd[2][kk] = Fwd[0][kk] * V * sin(theta) *
242  exp(-1 + time / timeramp);
243  }
244  else
245  {
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);
249  }
250 
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]);
255 
256  errV = 1.0;
257  errTheta = 1.0;
258  theta = M_PI / 4.0;
259  V = kExt * sin(theta);
260  }
261 
262  for (int i = 0; i < nvariables; ++i)
263  {
264  Vmath::Vcopy(npoints, &Fwd[i][id2], 1,
265  &(m_fields[i]->GetBndCondExpansions()[m_bcRegion]->
266  UpdatePhys())[id1],1);
267  }
268  }
269 }
270 
271 }
STL namespace.
LibUtilities::SessionReaderSharedPtr m_session
Session reader.
Definition: CFSBndCond.h:84
Array< OneD, MultiRegions::ExpListSharedPtr > m_fields
Array of fields.
Definition: CFSBndCond.h:86
boost::shared_ptr< SessionReader > SessionReaderSharedPtr
Definition: MeshPartition.h:51
CFSBndCondFactory & GetCFSBndCondFactory()
Declaration of the boundary condition factory singleton.
Definition: CFSBndCond.cpp:42
virtual void v_Apply(Array< OneD, Array< OneD, NekDouble > > &Fwd, Array< OneD, Array< OneD, NekDouble > > &physarray, const NekDouble &time)
Encapsulates the user-defined boundary conditions for compressible flow solver.
Definition: CFSBndCond.h:71
double NekDouble
NekDouble m_gamma
Parameters of the flow.
Definition: CFSBndCond.h:95
int m_offset
Offset.
Definition: CFSBndCond.h:103
void Vcopy(int n, const T *x, const int incx, T *y, const int incy)
Definition: Vmath.cpp:1061
int m_bcRegion
Id of the boundary region.
Definition: CFSBndCond.h:101