Nektar++
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 // Permission is hereby granted, free of charge, to any person obtaining a
14 // copy of this software and associated documentation files (the "Software"),
15 // to deal in the Software without restriction, including without limitation
16 // the rights to use, copy, modify, merge, publish, distribute, sublicense,
17 // and/or sell copies of the Software, and to permit persons to whom the
18 // Software is furnished to do so, subject to the following conditions:
19 //
20 // The above copyright notice and this permission notice shall be included
21 // in all copies or substantial portions of the Software.
22 //
23 // THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
24 // OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
25 // FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
26 // THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
27 // LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
28 // FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
29 // DEALINGS IN THE SOFTWARE.
30 //
31 // Description: Ringleb flow boundary condition
32 //
33 ///////////////////////////////////////////////////////////////////////////////
34 
35 #include "RinglebFlowBC.h"
36 #include <MultiRegions/AssemblyMap/AssemblyMapDG.h>
37 
38 using namespace std;
39 
40 namespace Nektar
41 {
42 
43 std::string RinglebFlowBC::className =
44  GetCFSBndCondFactory().RegisterCreatorFunction(
45  "RinglebFlow", RinglebFlowBC::create,
46  "Ringleb flow boundary condition.");
47 
48 RinglebFlowBC::RinglebFlowBC(
49  const LibUtilities::SessionReaderSharedPtr &pSession,
50  const Array<OneD, MultiRegions::ExpListSharedPtr> &pFields,
51  const Array<OneD, Array<OneD, NekDouble>> &pTraceNormals,
52  const int pSpaceDim, const int bcRegion, const int cnt)
53  : CFSBndCond(pSession, pFields, pTraceNormals, pSpaceDim, bcRegion, cnt)
54 {
55  m_expdim = pFields[0]->GetGraph()->GetMeshDimension();
56 
57  m_homo1D = false;
58  if (m_session->DefinesSolverInfo("HOMOGENEOUS"))
59  {
60  std::string HomoStr = m_session->GetSolverInfo("HOMOGENEOUS");
61  if ((HomoStr == "HOMOGENEOUS1D") || (HomoStr == "Homogeneous1D") ||
62  (HomoStr == "1D") || (HomoStr == "Homo1D"))
63  {
64  m_homo1D = true;
65  }
66  }
67 }
68 
69 void RinglebFlowBC::v_Apply(Array<OneD, Array<OneD, NekDouble>> &Fwd,
70  Array<OneD, Array<OneD, NekDouble>> &physarray,
71  const NekDouble &time)
72 {
73  int nvariables = physarray.size();
74 
75  // For 3DHomogenoeus1D
76  int n_planes = 1;
77  if (m_expdim == 2 && m_homo1D)
78  {
79  int nPointsTot = m_fields[0]->GetTotPoints();
80  int nPointsTot_plane = m_fields[0]->GetPlane(0)->GetTotPoints();
81  n_planes = nPointsTot / nPointsTot_plane;
82  }
83 
84  int id2, id2_plane, e_max;
85 
86  e_max = m_fields[0]->GetBndCondExpansions()[m_bcRegion]->GetExpSize();
87 
88  for (int e = 0; e < e_max; ++e)
89  {
90  int npoints = m_fields[0]
91  ->GetBndCondExpansions()[m_bcRegion]
92  ->GetExp(e)
93  ->GetTotPoints();
94  int id1 =
95  m_fields[0]->GetBndCondExpansions()[m_bcRegion]->GetPhys_Offset(e);
96 
97  // For 3DHomogenoeus1D
98  if (m_expdim == 2 && m_homo1D)
99  {
100  int m_offset_plane = m_offset / n_planes;
101  int e_plane;
102  int e_max_plane = e_max / n_planes;
103  int nTracePts_plane = m_fields[0]->GetTrace()->GetNpoints();
104 
105  int planeID = floor((e + 0.5) / e_max_plane);
106  e_plane = e - e_max_plane * planeID;
107 
108  id2_plane = m_fields[0]->GetTrace()->GetPhys_Offset(
109  m_fields[0]->GetTraceMap()->GetBndCondIDToGlobalTraceID(
110  m_offset_plane + e_plane));
111  id2 = id2_plane + planeID * nTracePts_plane;
112  }
113  else // For general case
114  {
115  id2 = m_fields[0]->GetTrace()->GetPhys_Offset(
116  m_fields[0]->GetTraceMap()->GetBndCondIDToGlobalTraceID(
117  m_offset + e));
118  }
119 
120  Array<OneD, NekDouble> x0(npoints, 0.0);
121  Array<OneD, NekDouble> x1(npoints, 0.0);
122  Array<OneD, NekDouble> x2(npoints, 0.0);
123 
124  m_fields[0]->GetBndCondExpansions()[m_bcRegion]->GetExp(e)->GetCoords(
125  x0, x1, x2);
126 
127  // Flow parameters
128  NekDouble c, k, phi, r, J, VV, pp, sint, P, ss;
129  NekDouble J11, J12, J21, J22, det;
130  NekDouble Fx, Fy;
131  NekDouble xi, yi;
132  NekDouble dV;
133  NekDouble dtheta;
134  NekDouble par1;
135  NekDouble theta = M_PI / 4.0;
136  NekDouble kExt = 0.7;
137  NekDouble V = kExt * sin(theta);
138  NekDouble toll = 1.0e-8;
139  NekDouble errV = 1.0;
140  NekDouble errTheta = 1.0;
141  NekDouble gamma = m_gamma;
142  NekDouble gamma_1_2 = (gamma - 1.0) / 2.0;
143 
144  // Loop on all the points of that edge
145  for (int j = 0; j < npoints; j++)
146  {
147 
148  while ((abs(errV) > toll) || (abs(errTheta) > toll))
149  {
150  VV = V * V;
151  sint = sin(theta);
152  c = sqrt(1.0 - gamma_1_2 * VV);
153  k = V / sint;
154  phi = 1.0 / k;
155  pp = phi * phi;
156  J = 1.0 / c + 1.0 / (3.0 * c * c * c) +
157  1.0 / (5.0 * c * c * c * c * c) -
158  0.5 * log((1.0 + c) / (1.0 - c));
159 
160  r = pow(c, 1.0 / gamma_1_2);
161  xi = 1.0 / (2.0 * r) * (1.0 / VV - 2.0 * pp) + J / 2.0;
162  yi = phi / (r * V) * sqrt(1.0 - VV * pp);
163  par1 = 25.0 - 5.0 * VV;
164  ss = sint * sint;
165 
166  Fx = xi - x0[j];
167  Fy = yi - x1[j];
168 
169  J11 =
170  39062.5 / pow(par1, 3.5) * (1.0 / VV - 2.0 / VV * ss) * V +
171  1562.5 / pow(par1, 2.5) *
172  (-2.0 / (VV * V) + 4.0 / (VV * V) * ss) +
173  12.5 / pow(par1, 1.5) * V + 312.5 / pow(par1, 2.5) * V +
174  7812.5 / pow(par1, 3.5) * V -
175  0.25 *
176  (-1.0 / pow(par1, 0.5) * V /
177  (1.0 - 0.2 * pow(par1, 0.5)) -
178  (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) /
181  (1.0 + 0.2 * pow(par1, 0.5)) *
182  (1.0 - 0.2 * pow(par1, 0.5));
183 
184  J12 = -6250.0 / pow(par1, 2.5) / VV * sint * cos(theta);
185  J21 =
186  -6250.0 / (VV * V) * sint / pow(par1, 2.5) *
187  pow((1.0 - ss), 0.5) +
188  78125.0 / V * sint / pow(par1, 3.5) * pow((1.0 - ss), 0.5);
189 
190  // the matrix is singular when theta = pi/2
191  if (abs(x1[j]) < toll && abs(cos(theta)) < toll)
192  {
193  J22 = -39062.5 / pow(par1, 3.5) / V +
194  3125 / pow(par1, 2.5) / (VV * V) +
195  12.5 / pow(par1, 1.5) * V +
196  312.5 / pow(par1, 2.5) * V +
197  7812.5 / pow(par1, 3.5) * V -
198  0.25 *
199  (-1.0 / pow(par1, 0.5) * V /
200  (1.0 - 0.2 * pow(par1, 0.5)) -
201  (1.0 + 0.2 * pow(par1, 0.5)) /
202  pow((1.0 - 0.2 * pow(par1, 0.5)), 2.0) /
203  pow(par1, 0.5) * V) /
204  (1.0 + 0.2 * pow(par1, 0.5)) *
205  (1.0 - 0.2 * pow(par1, 0.5));
206 
207  // dV = -dV/dx * Fx
208  dV = -1.0 / J22 * Fx;
209  dtheta = 0.0;
210  theta = M_PI / 2.0;
211  }
212  else
213  {
214  J22 = 3125.0 / VV * cos(theta) / pow(par1, 2.5) *
215  pow((1.0 - ss), 0.5) -
216  3125.0 / VV * ss / pow(par1, 2.5) /
217  pow((1.0 - ss), 0.5) * cos(theta);
218 
219  det = -1.0 / (J11 * J22 - J12 * J21);
220 
221  // [dV dtheta]' = -[invJ]*[Fx Fy]'
222  dV = det * (J22 * Fx - J12 * Fy);
223  dtheta = det * (-J21 * Fx + J11 * Fy);
224  }
225 
226  V = V + dV;
227  theta = theta + dtheta;
228 
229  errV = abs(dV);
230  errTheta = abs(dtheta);
231  }
232 
233  c = sqrt(1.0 - gamma_1_2 * VV);
234  int kk = id2 + j;
235  NekDouble timeramp = 200.0;
236  ;
237  if (time < timeramp &&
238  !(m_session->DefinesFunction("InitialConditions") &&
239  m_session->GetFunctionType("InitialConditions", 0) ==
240  LibUtilities::eFunctionTypeFile))
241  {
242  Fwd[0][kk] =
243  pow(c, 1.0 / gamma_1_2) * exp(-1.0 + time / timeramp);
244 
245  Fwd[1][kk] =
246  Fwd[0][kk] * V * cos(theta) * exp(-1 + time / timeramp);
247 
248  Fwd[2][kk] =
249  Fwd[0][kk] * V * sin(theta) * exp(-1 + time / timeramp);
250  }
251  else
252  {
253  Fwd[0][kk] = pow(c, 1.0 / gamma_1_2);
254  Fwd[1][kk] = Fwd[0][kk] * V * cos(theta);
255  Fwd[2][kk] = Fwd[0][kk] * V * sin(theta);
256  }
257 
258  P = (c * c) * Fwd[0][kk] / gamma;
259  Fwd[3][kk] = P / (gamma - 1.0) +
260  0.5 * (Fwd[1][kk] * Fwd[1][kk] / Fwd[0][kk] +
261  Fwd[2][kk] * Fwd[2][kk] / Fwd[0][kk]);
262 
263  errV = 1.0;
264  errTheta = 1.0;
265  theta = M_PI / 4.0;
266  V = kExt * sin(theta);
267  }
268 
269  for (int i = 0; i < nvariables; ++i)
270  {
271  Vmath::Vcopy(npoints, &Fwd[i][id2], 1,
272  &(m_fields[i]
273  ->GetBndCondExpansions()[m_bcRegion]
274  ->UpdatePhys())[id1],
275  1);
276  }
277  }
278 }
279 
280 } // namespace Nektar
Nektar::CFSBndCond
Encapsulates the user-defined boundary conditions for compressible flow solver.
Definition: CFSBndCond.h:70
Nektar::CFSBndCond::m_session
LibUtilities::SessionReaderSharedPtr m_session
Session reader.
Definition: CFSBndCond.h:89
Nektar::CFSBndCond::m_gamma
NekDouble m_gamma
Parameters of the flow.
Definition: CFSBndCond.h:102
Nektar::CFSBndCond::m_bcRegion
int m_bcRegion
Id of the boundary region.
Definition: CFSBndCond.h:109
Nektar::CFSBndCond::m_offset
int m_offset
Offset.
Definition: CFSBndCond.h:111
Nektar::CFSBndCond::m_fields
Array< OneD, MultiRegions::ExpListSharedPtr > m_fields
Array of fields.
Definition: CFSBndCond.h:91
Nektar::LibUtilities::NekFactory::RegisterCreatorFunction
tKey RegisterCreatorFunction(tKey idKey, CreatorFunction classCreator, std::string pDesc="")
Register a class with the factory.
Definition: NekFactory.hpp:198
Nektar::RinglebFlowBC::v_Apply
virtual void v_Apply(Array< OneD, Array< OneD, NekDouble >> &Fwd, Array< OneD, Array< OneD, NekDouble >> &physarray, const NekDouble &time) override
Definition: RinglebFlowBC.cpp:69
Nektar::LibUtilities::SessionReaderSharedPtr
std::shared_ptr< SessionReader > SessionReaderSharedPtr
Definition: SessionReader.h:115
Nektar
The above copyright notice and this permission notice shall be included.
Definition: CoupledSolver.h:2
Nektar::GetCFSBndCondFactory
CFSBndCondFactory & GetCFSBndCondFactory()
Declaration of the boundary condition factory singleton.
Definition: CFSBndCond.cpp:41
Nektar::NekDouble
double NekDouble
Definition: NektarUnivTypeDefs.hpp:43
Vmath::Vcopy
void Vcopy(int n, const T *x, const int incx, T *y, const int incy)
Definition: Vmath.cpp:1255
tinysimd::log
scalarT< T > log(scalarT< T > in)
Definition: scalar.hpp:303
tinysimd::abs
scalarT< T > abs(scalarT< T > in)
Definition: scalar.hpp:298
tinysimd::sqrt
scalarT< T > sqrt(scalarT< T > in)
Definition: scalar.hpp:294
Generated on Sat Mar 25 2023 12:35:52 for Nektar++ by doxygen 1.9.1