Nektar++
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
QInflow.cpp
Go to the documentation of this file.
1 ///////////////////////////////////////////////////////////////////////////////
2 //
3 // File QInflow.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: QInflow class
33 //
34 ///////////////////////////////////////////////////////////////////////////////
35 
37 
38 namespace Nektar
39 {
40 
41  std::string QInflow::className
43  "Qinflow",
45  "Inflow boundary condition");
46 
47  /**
48  *
49  */
50  QInflow::QInflow(Array<OneD, MultiRegions::ExpListSharedPtr> pVessel,
52  PulseWavePressureAreaSharedPtr pressureArea)
53  : PulseWaveBoundary(pVessel,pSession,pressureArea)
54  {
55  }
56 
57  /**
58  *
59  */
61  {
62 
63  }
64 
66  const Array<OneD,const Array<OneD, NekDouble> > &inarray,
67  Array<OneD, Array<OneD, NekDouble> > &A_0,
68  Array<OneD, Array<OneD, NekDouble> > &beta,
69  const NekDouble time,
70  int omega,
71  int offset,
72  int n)
73  {
74  NekDouble Q, A_u, u_u;
75  NekDouble A_r, u_r, A_l, u_l;
76  Array<OneD, MultiRegions::ExpListSharedPtr> vessel(2);
77 
78  vessel[0] = m_vessels[2*omega];
79  vessel[1] = m_vessels[2*omega+1];
80 
81  vessel[0]->EvaluateBoundaryConditions(time);
82 
83  // Note: The Q value is contained in A in the
84  // inputfile, the value in u has to be 1.0
85  //ASSERTL0((vessel[1]->UpdateBndCondExpansion(n))->UpdatePhys()[0] == 1.0, "For the Q-inflow BC the value in u must be 1.0");
86 
87  // Get the values of all variables needed for the Riemann problem
88  Q = (vessel[0]->UpdateBndCondExpansion(n))->GetCoeffs()[0];
89 
90  A_r = inarray[0][offset];
91  u_r = inarray[1][offset];
92 
93  // Call the Q-inflow Riemann solver
94  Q_RiemannSolver(Q,A_r,u_r,A_0[omega][0],beta[omega][0],A_u,u_u);
95 
96  // Set the boundary conditions to prescribe
97  A_l=A_r;
98  u_l=2*u_u-u_r;
99 
100  // Store the updated values in the boundary condition
101  (vessel[0]->UpdateBndCondExpansion(n))->UpdatePhys()[0] = A_l;
102  (vessel[1]->UpdateBndCondExpansion(n))->UpdatePhys()[0] = u_l;
103  }
104 
105  /**
106  * Q-inflow Riemann solver for pulse wave propagation. This Riemann solver is called
107  * by DoOdeProjection in case of a Q-inflow boundary condition. It is based on the
108  * conservation of mass and total pressure and on the characteristic information. For
109  * further details see "Pulse wave propagation in the human vascular system", section 3.4.1
110  * Returns the upwinded quantities \f$(A_u,u_u)\f$ and stores them into the boundary values
111  */
113  NekDouble &Au,NekDouble &uu)
114  {
115  NekDouble W2 = 0.0;
116  NekDouble A_calc = 0.0;
117  NekDouble fa = 0.0;
118  NekDouble dfa = 0.0;
119  NekDouble delta_A_calc = 0.0;
120  NekDouble rho = m_rho;
121 
122  int proceed = 1;
123  int iter = 0;
124  int MAX_ITER = 200;
125 
126  // Tolerances for the algorithm
127  NekDouble Tol = 1.0e-10;
128 
129  // Riemann invariant \f$W_2(Ar,ur)\f$
130  W2 = u_r - 4*sqrt(beta/(2*rho))*sqrt(sqrt(A_r));
131 
132  // Newton Iteration (Area only)
133  A_calc = A_r;
134  while ((proceed) && (iter < MAX_ITER))
135  {
136  iter =iter+1;
137 
138  fa = Q - W2*A_calc - A_calc*4*sqrt(beta/(2*rho))*sqrt(sqrt(A_calc));
139  dfa = -W2 - 5*sqrt(beta/(2*rho))*sqrt(sqrt(A_calc));
140  delta_A_calc = fa/dfa;
141  A_calc = A_calc - delta_A_calc;
142 
143  if (sqrt(delta_A_calc*delta_A_calc) < Tol)
144  proceed = 0;
145  }
146 
147  // Obtain u_u and A_u
148  uu = W2+4*sqrt(beta/(2*rho))*sqrt(sqrt(A_calc));
149  Au = A_calc;
150  }
151 
152 }