Nektar++
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
PulseWavePropagation.cpp
Go to the documentation of this file.
1 ///////////////////////////////////////////////////////////////////////////////
2 //
3 // File PulseWavePropagation.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: Pulse Wave Propagation solve routines based on the weak formulation (1):
33 //
34 ///////////////////////////////////////////////////////////////////////////////
35 
36 #include <iostream>
37 
39 
40 using namespace std;
41 
42 namespace Nektar
43 {
44  string PulseWavePropagation::className = GetEquationSystemFactory().RegisterCreatorFunction("PulseWavePropagation", PulseWavePropagation::create, "Pulse Wave Propagation equation.");
45  /**
46  * @class PulseWavePropagation
47  *
48  * Set up the routines based on the weak formulation from
49  * "Computational Modelling of 1D blood flow with variable
50  * mechanical properties" by S. J. Sherwin et al. The weak
51  * formulation (1) reads:
52  * \f$ \sum_{e=1}^{N_{el}} \left[ \left( \frac{\partial \mathbf{U}^{\delta} }{\partial t} ,
53  * \mathbf{\psi}^{\delta} \right)_{\Omega_e} - \left( \frac{\partial \mathbf{F(\mathbf{U})}^{\delta} }
54  * {\partial x}, \mathbf{\psi}^{\delta} \right)_{\Omega_e} + \left[ \mathbf{\psi}^{\delta}
55  * \cdot \{ \mathbf{F}^u - \mathbf{F}(\mathbf{U}^{\delta}) \} \right]_{x_e^l}^{x_eû} \right] = 0 \f$
56  */
57  PulseWavePropagation::PulseWavePropagation(const LibUtilities::SessionReaderSharedPtr& pSession)
58  : PulseWaveSystem(pSession)
59  {
60  }
61 
63  {
65 
67  m_pressureArea->DoPressure();
68 
70  {
73  }
74  else
75  {
76  ASSERTL0(false, "Implicit Pulse Wave Propagation not set up.");
77  }
78  }
79 
81  {
82  }
83 
84 
85  /**
86  * Computes the right hand side of (1). The RHS is everything
87  * except the term that contains the time derivative
88  * \f$\frac{\partial \mathbf{U}}{\partial t}\f$. In case of a
89  * Discontinuous Galerkin projection, the routine WeakDGAdvection
90  * will be called which then calls v_GetFluxVector and
91  * v_NumericalFlux implemented in the PulseWavePropagation class.
92  *
93  */
95  Array<OneD, Array<OneD, NekDouble> >&outarray,
96  const NekDouble time)
97  {
98  int i;
99 
102 
103  Array<OneD, NekDouble> tmpArray;
104 
105  int cnt = 0;
106 
107  // Set up Inflow and Outflow boundary conditions.
108  SetPulseWaveBoundaryConditions(inarray, outarray, time);
109 
110  // Set up any interface conditions and write into boundary condition
112 
113  // do advection evauation in all domains
114  for(int omega=0; omega < m_nDomains; ++omega)
115  {
116  m_currentDomain = omega;
117  int nq = m_vessels[omega*m_nVariables]->GetTotPoints();
118  int ncoeffs = m_vessels[omega*m_nVariables]->GetNcoeffs();
119 
120  for (i = 0; i < m_nVariables; ++i)
121  {
122  physarray[i] = inarray[i]+cnt;
123  modarray[i] = Array<OneD, NekDouble>(ncoeffs);
124  }
125 
126  for(i = 0; i < m_nVariables; ++i)
127  {
128  m_fields[i] = m_vessels[omega*m_nVariables+ i];
129  }
130 
131  WeakDGAdvection(physarray, modarray, true, true);
132 
133  for(i = 0; i < m_nVariables; ++i)
134  {
135  Vmath::Neg(ncoeffs,modarray[i],1);
136  }
137 
138  for(i = 0; i < m_nVariables; ++i)
139  {
140  m_vessels[omega*m_nVariables+i]->MultiplyByElmtInvMass(modarray[i],modarray[i]);
141  m_vessels[omega*m_nVariables+i]->BwdTrans(modarray[i],tmpArray = outarray[i]+cnt);
142  }
143  cnt += nq;
144  }
145  }
146 
148  Array<OneD, Array<OneD, NekDouble> >&outarray,
149  const NekDouble time)
150  {
151  // Just copy over array
152  for(int i = 0; i < m_nVariables; ++i)
153  {
154  Vmath::Vcopy(inarray[i].num_elements(),inarray[i],1,outarray[i],1);
155  }
156  }
157 
158 
159 
160  /**
161  * Does the projection between ... space and the ... space. Also checks for Q-inflow boundary
162  * conditions at the inflow of the current arterial segment and applies the Q-inflow if specified
163  */
165  const Array<OneD,const Array<OneD, NekDouble> >&inarray,
166  Array<OneD, Array<OneD, NekDouble> >&outarray,
167  const NekDouble time)
168 
169  {
170  int omega;
171 
173 
174  int offset=0;
175 
176  //This will be moved to the RCR boundary condition once factory is setup
177  if (time == 0)
178  {
180 
181  for(omega = 0; omega < m_nDomains; ++omega)
182  {
183  vessel[0] = m_vessels[2*omega];
184 
185  for(int j = 0; j < 2; ++j)
186  {
187  std::string BCType =vessel[0]->GetBndConditions()[j]->GetUserDefined();
188  if(BCType.empty()) // if not condition given define it to be NoUserDefined
189  {
190  BCType = "NoUserDefined";
191  }
192 
194 
195  // turn on time depedent BCs
196  if(BCType == "Q-inflow")
197  {
198  vessel[0]->GetBndConditions()[j]->SetIsTimeDependent(true);
199  }
200  else if(BCType == "RCR-terminal")
201  {
202  vessel[0]->GetBndConditions()[j]->SetIsTimeDependent(true);
203  }
204  }
205  }
206 
207  }
208 
209  SetBoundaryConditions(time);
210 
211  // Loop over all vessesls and set boundary conditions
212  for(omega = 0; omega < m_nDomains; ++omega)
213  {
214  for(int n = 0; n < 2; ++n)
215  {
216  m_Boundary[2*omega+n]->DoBoundary(inarray,m_A_0,m_beta,time,omega,offset,n);
217  }
218  offset += m_vessels[2*omega]->GetTotPoints();
219  }
220 
221  }
222 
223 
224  /**
225  * Calculates the second term of the weak form (1): \f$
226  * \left( \frac{\partial \mathbf{F(\mathbf{U})}^{\delta}
227  * }{\partial x}, \mathbf{\psi}^{\delta} \right)_{\Omega_e}
228  * \f$
229  * The variables of the system are $\mathbf{U} = [A,u]^T$
230  * physfield[0] = A physfield[1] = u
231  * flux[0] = F[0] = A*u flux[1] = F[1] = u^2/2 + p/rho
232  * p-A-relationship: p = p_ext + beta*(sqrt(A)-sqrt(A_0))
233  */
236  {
237  int nq = m_vessels[m_currentDomain*m_nVariables]->GetTotPoints();
238  NekDouble p = 0.0;
239  NekDouble p_t = 0.0;
240 
241  switch (i)
242  {
243  case 0: // Flux for A equation
244  {
245  for (int j = 0; j < nq; j++)
246  {
247  flux[0][j] = physfield[0][j]*physfield[1][j];
248  }
249  }
250  break;
251  case 1: // Flux for u equation
252  {
253  for (int j = 0; j < nq; j++)
254  {
255  ASSERTL0(physfield[0][j]>=0,"Negative A not allowed.");
256 
257  p = m_pext + m_beta[m_currentDomain][j]*
258  (sqrt(physfield[0][j]) - sqrt(m_A_0[m_currentDomain][j]));
259 
260  p_t = (physfield[1][j]*physfield[1][j])/2 + p/m_rho;
261  flux[0][j] = p_t;
262  }
263  }
264  break;
265  default:
266  ASSERTL0(false,"GetFluxVector: illegal vector index");
267  break;
268  }
269  }
270 
271 
272  /**
273  * Calculates the third term of the weak form (1): numerical flux
274  * at boundary \f$ \left[ \mathbf{\psi}^{\delta} \cdot \{
275  * \mathbf{F}^u - \mathbf{F}(\mathbf{U}^{\delta}) \}
276  * \right]_{x_e^l}^{x_eû} \f$
277  */
279  Array<OneD, Array<OneD, NekDouble> > &numflux)
280  {
281  int i;
282  int nTracePts = GetTraceTotPoints();
283 
286 
287  for (i = 0; i < m_nVariables; ++i)
288  {
289  Fwd[i] = Array<OneD, NekDouble>(nTracePts);
290  Bwd[i] = Array<OneD, NekDouble>(nTracePts);
291  }
292 
293  // Get the physical values at the trace
294  for (i = 0; i < m_nVariables; ++i)
295  {
296  m_vessels[m_currentDomain*m_nVariables+ i]->
297  GetFwdBwdTracePhys(physfield[i],Fwd[i],Bwd[i]);
298  }
299 
300  // Solve the upwinding Riemann problem within one arterial
301  // segment by calling the upwinding Riemann solver implemented
302  // in this file
303  NekDouble Aflux, uflux;
304  for (i = 0; i < nTracePts; ++i)
305  {
306  switch(m_upwindTypePulse)
307  {
308  case eUpwindPulse:
309  {
310  RiemannSolverUpwind(Fwd[0][i],Fwd[1][i],Bwd[0][i],Bwd[1][i],
311  Aflux, uflux, m_A_0_trace[m_currentDomain][i],
314  }
315  break;
316  default:
317  {
318  ASSERTL0(false,"populate switch statement for upwind flux");
319  }
320  break;
321  }
322  numflux[0][i] = Aflux;
323  numflux[1][i] = uflux;
324  }
325  }
326 
327 
328  /**
329  * Riemann solver for upwinding at an interface between two
330  * elements. Uses the characteristic variables for calculating
331  * the upwinded state \f$(A_u,u_u)\f$ from the left
332  * \f$(A_L,u_L)\f$ and right state \f$(A_R,u_R)\f$. Returns the
333  * upwinded flux $\mathbf{F}^u$ needed for the weak formulation
334  * (1). Details can be found in "Pulse wave propagation in the
335  * human vascular system", section 3.3
336  *
337  */
339  NekDouble AR,NekDouble uR,
340  NekDouble &Aflux,
341  NekDouble &uflux,
342  NekDouble A_0,
343  NekDouble beta,
344  NekDouble n)
345  {
347  Array<OneD, NekDouble> upwindedphysfield(2);
348  NekDouble cL = 0.0;
349  NekDouble cR = 0.0;
350  NekDouble rho = m_rho;
351  NekDouble pext = m_pext;
352  NekDouble p = 0.0;
353  NekDouble p_t = 0.0;
354 
355  // Compute the wave speeds. The use of the normal here allows
356  // for the definition of the characteristics to be inverted
357  // (and hence the left and right state) if n is in the -ve
358  // x-direction. This means we end up with the positive
359  // defintion of the flux which has to therefore be multiplied
360  // by the normal at the end of the methods This is a bit of a
361  // mind twister but is efficient from a coding perspective.
362  cL = sqrt(beta*sqrt(AL)/(2*rho))*n;
363  cR = sqrt(beta*sqrt(AR)/(2*rho))*n;
364 
365  ASSERTL1(fabs(cL+cR) > fabs(uL+uR),"Conditions are not sub-sonic");
366 
367  // If upwinding from left and right for subsonic domain
368  // then know characteristics immediately
369  W[0] = uL + 4*cL;
370  W[1] = uR - 4*cR;
371 
372  // Calculate conservative variables from characteristics
373  NekDouble w0mw1 = 0.25*(W[0]-W[1]);
374  NekDouble fac = rho/(2*beta);
375  w0mw1 *= w0mw1; // squared
376  w0mw1 *= w0mw1; // fourth power
377  fac *= fac; // squared
378  upwindedphysfield[0]= w0mw1*fac;
379  upwindedphysfield[1]= 0.5*(W[0] + W[1]);
380 
381  // Compute the fluxes multipled by the normal.
382  Aflux = upwindedphysfield[0] * upwindedphysfield[1]*n;
383  p = pext + beta*(sqrt(upwindedphysfield[0]) - sqrt(A_0));
384  p_t = 0.5*(upwindedphysfield[1]*upwindedphysfield[1]) + p/rho;
385  uflux = p_t*n;
386  }
387 
388  /**
389  * Print summary routine, calls virtual routine reimplemented in
390  * UnsteadySystem
391  */
393  {
395  }
396 }
Array< OneD, Array< OneD, NekDouble > > m_beta_trace
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:161
UpwindTypePulse m_upwindTypePulse
tBaseSharedPtr CreateInstance(tKey idKey BOOST_PP_COMMA_IF(MAX_PARAM) BOOST_PP_ENUM_BINARY_PARAMS(MAX_PARAM, tParam, x))
Create an instance of the class referred to by idKey.
Definition: NekFactory.hpp:162
void SetPulseWaveBoundaryConditions(const Array< OneD, const Array< OneD, NekDouble > > &inarray, Array< OneD, Array< OneD, NekDouble > > &outarray, const NekDouble time)
LibUtilities::TimeIntegrationSchemeOperators m_ode
The time integration scheme operators to use.
std::vector< std::pair< std::string, std::string > > SummaryList
Definition: Misc.h:47
Array< OneD, Array< OneD, NekDouble > > m_trace_fwd_normal
STL namespace.
boost::shared_ptr< SessionReader > SessionReaderSharedPtr
Definition: MeshPartition.h:51
void EnforceInterfaceConditions(const Array< OneD, const Array< OneD, NekDouble > > &fields)
virtual void v_GenerateSummary(SolverUtils::SummaryList &s)
BoundaryFactory & GetBoundaryFactory()
Array< OneD, Array< OneD, NekDouble > > m_A_0_trace
void DoOdeRhs(const Array< OneD, const Array< OneD, NekDouble > > &inarray, Array< OneD, Array< OneD, NekDouble > > &outarray, const NekDouble time)
Base class for unsteady solvers.
void DefineProjection(FuncPointerT func, ObjectPointerT obj)
virtual SOLVER_UTILS_EXPORT void v_GenerateSummary(SummaryList &s)
Print a summary of time stepping parameters.
void DoOdeProjection(const Array< OneD, const Array< OneD, NekDouble > > &inarray, Array< OneD, Array< OneD, NekDouble > > &outarray, const NekDouble time)
bool m_explicitAdvection
Indicates if explicit or implicit treatment of advection is used.
void DefineOdeRhs(FuncPointerT func, ObjectPointerT obj)
virtual void v_NumericalFlux(Array< OneD, Array< OneD, NekDouble > > &physfield, Array< OneD, Array< OneD, NekDouble > > &numflux)
DG Pulse Wave Propagation routines: Numerical Flux at interelemental boundaries.
virtual void v_GetFluxVector(const int i, Array< OneD, Array< OneD, NekDouble > > &physfield, Array< OneD, Array< OneD, NekDouble > > &flux)
Array< OneD, Array< OneD, NekDouble > > m_A_0
void Neg(int n, T *x, const int incx)
Negate x = -x.
Definition: Vmath.cpp:382
double NekDouble
EquationSystemFactory & GetEquationSystemFactory()
SOLVER_UTILS_EXPORT void SetBoundaryConditions(NekDouble time)
Evaluates the boundary conditions at the given time.
SOLVER_UTILS_EXPORT int GetTraceTotPoints()
Array< OneD, Array< OneD, NekDouble > > m_beta
Array< OneD, MultiRegions::ExpListSharedPtr > m_fields
Array holding all dependent variables.
LibUtilities::SessionReaderSharedPtr m_session
The session reader.
PressureAreaFactory & GetPressureAreaFactory()
Array< OneD, PulseWaveBoundarySharedPtr > m_Boundary
SOLVER_UTILS_EXPORT void WeakDGAdvection(const Array< OneD, Array< OneD, NekDouble > > &InField, Array< OneD, Array< OneD, NekDouble > > &OutField, bool NumericalFluxIncludesNormal=true, bool InFieldIsInPhysSpace=false, int nvariables=0)
Calculate the weak discontinuous Galerkin advection.
#define ASSERTL1(condition, msg)
Assert Level 1 – Debugging which is used whether in FULLDEBUG or DEBUG compilation mode...
Definition: ErrorUtil.hpp:191
void Vcopy(int n, const T *x, const int incx, T *y, const int incy)
Definition: Vmath.cpp:1047
Array< OneD, MultiRegions::ExpListSharedPtr > m_vessels
PulseWavePressureAreaSharedPtr m_pressureArea
void RiemannSolverUpwind(NekDouble AL, NekDouble uL, NekDouble AR, NekDouble uR, NekDouble &Aflux, NekDouble &uflux, NekDouble A_0, NekDouble beta, NekDouble n)
Upwinding Riemann solver for interelemental boundaries.
simple upwinding scheme
tKey RegisterCreatorFunction(tKey idKey, CreatorFunction classCreator, tDescription pDesc="")
Register a class with the factory.
Definition: NekFactory.hpp:215