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