Nektar++
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 namespace Nektar
41 {
42  string PulseWavePropagation::className = GetEquationSystemFactory().RegisterCreatorFunction("PulseWavePropagation", PulseWavePropagation::create, "Pulse Wave Propagation equation.");
43  /**
44  * @class PulseWavePropagation
45  *
46  * Set up the routines based on the weak formulation from
47  * "Computational Modelling of 1D blood flow with variable
48  * mechanical properties" by S. J. Sherwin et al. The weak
49  * formulation (1) reads:
50  * \f$ \sum_{e=1}^{N_{el}} \left[ \left( \frac{\partial \mathbf{U}^{\delta} }{\partial t} ,
51  * \mathbf{\psi}^{\delta} \right)_{\Omega_e} - \left( \frac{\partial \mathbf{F(\mathbf{U})}^{\delta} }
52  * {\partial x}, \mathbf{\psi}^{\delta} \right)_{\Omega_e} + \left[ \mathbf{\psi}^{\delta}
53  * \cdot \{ \mathbf{F}^u - \mathbf{F}(\mathbf{U}^{\delta}) \} \right]_{x_e^l}^{x_eû} \right] = 0 \f$
54  */
56  : PulseWaveSystem(pSession)
57  {
58  }
59 
61  {
63 
65  m_pressureArea->DoPressure();
66 
68  {
71  }
72  else
73  {
74  ASSERTL0(false, "Implicit Pulse Wave Propagation not set up.");
75  }
76  }
77 
79  {
80  }
81 
82 
83  /**
84  * Computes the right hand side of (1). The RHS is everything
85  * except the term that contains the time derivative
86  * \f$\frac{\partial \mathbf{U}}{\partial t}\f$. In case of a
87  * Discontinuous Galerkin projection, the routine WeakDGAdvection
88  * will be called which then calls v_GetFluxVector and
89  * v_NumericalFlux implemented in the PulseWavePropagation class.
90  *
91  */
93  Array<OneD, Array<OneD, NekDouble> >&outarray,
94  const NekDouble time)
95  {
96  int i;
97 
100 
101  Array<OneD, NekDouble> tmpArray;
102 
103  int cnt = 0;
104 
105  // Set up Inflow and Outflow boundary conditions.
106  SetPulseWaveBoundaryConditions(inarray, outarray, time);
107 
108  // Set up any interface conditions and write into boundary condition
110 
111  // do advection evauation in all domains
112  for(int omega=0; omega < m_nDomains; ++omega)
113  {
114  m_currentDomain = omega;
115  int nq = m_vessels[omega*m_nVariables]->GetTotPoints();
116  int ncoeffs = m_vessels[omega*m_nVariables]->GetNcoeffs();
117 
118  for (i = 0; i < m_nVariables; ++i)
119  {
120  physarray[i] = inarray[i]+cnt;
121  modarray[i] = Array<OneD, NekDouble>(ncoeffs);
122  }
123 
124  for(i = 0; i < m_nVariables; ++i)
125  {
126  m_fields[i] = m_vessels[omega*m_nVariables+ i];
127  }
128 
129  WeakDGAdvection(physarray, modarray, true, true);
130 
131  for(i = 0; i < m_nVariables; ++i)
132  {
133  Vmath::Neg(ncoeffs,modarray[i],1);
134  }
135 
136  for(i = 0; i < m_nVariables; ++i)
137  {
138  m_vessels[omega*m_nVariables+i]->MultiplyByElmtInvMass(modarray[i],modarray[i]);
139  m_vessels[omega*m_nVariables+i]->BwdTrans(modarray[i],tmpArray = outarray[i]+cnt);
140  }
141  cnt += nq;
142  }
143  }
144 
146  Array<OneD, Array<OneD, NekDouble> >&outarray,
147  const NekDouble time)
148  {
149  // Just copy over array
150  for(int i = 0; i < m_nVariables; ++i)
151  {
152  Vmath::Vcopy(inarray[i].num_elements(),inarray[i],1,outarray[i],1);
153  }
154  }
155 
156 
157 
158  /**
159  * Does the projection between ... space and the ... space. Also checks for Q-inflow boundary
160  * conditions at the inflow of the current arterial segment and applies the Q-inflow if specified
161  */
163  const Array<OneD,const Array<OneD, NekDouble> >&inarray,
164  Array<OneD, Array<OneD, NekDouble> >&outarray,
165  const NekDouble time)
166 
167  {
168  int omega;
169 
171 
172  int offset=0;
173 
174  //This will be moved to the RCR boundary condition once factory is setup
175  if (time == 0)
176  {
178 
179  for(omega = 0; omega < m_nDomains; ++omega)
180  {
181  vessel[0] = m_vessels[2*omega];
182 
183  for(int j = 0; j < 2; ++j)
184  {
185  std::string BCType =vessel[0]->GetBndConditions()[j]->GetUserDefined();
186  if(BCType.empty()) // if not condition given define it to be NoUserDefined
187  {
188  BCType = "NoUserDefined";
189  }
190 
192 
193  // turn on time depedent BCs
194  if(BCType == "Q-inflow")
195  {
196  vessel[0]->GetBndConditions()[j]->SetIsTimeDependent(true);
197  }
198  else if(BCType == "RCR-terminal")
199  {
200  vessel[0]->GetBndConditions()[j]->SetIsTimeDependent(true);
201  }
202  }
203  }
204 
205  }
206 
207  SetBoundaryConditions(time);
208 
209  // Loop over all vessesls and set boundary conditions
210  for(omega = 0; omega < m_nDomains; ++omega)
211  {
212  for(int n = 0; n < 2; ++n)
213  {
214  m_Boundary[2*omega+n]->DoBoundary(inarray,m_A_0,m_beta,time,omega,offset,n);
215  }
216  offset += m_vessels[2*omega]->GetTotPoints();
217  }
218 
219  }
220 
221 
222  /**
223  * Calculates the second term of the weak form (1): \f$
224  * \left( \frac{\partial \mathbf{F(\mathbf{U})}^{\delta}
225  * }{\partial x}, \mathbf{\psi}^{\delta} \right)_{\Omega_e}
226  * \f$
227  * The variables of the system are $\mathbf{U} = [A,u]^T$
228  * physfield[0] = A physfield[1] = u
229  * flux[0] = F[0] = A*u flux[1] = F[1] = u^2/2 + p/rho
230  * p-A-relationship: p = p_ext + beta*(sqrt(A)-sqrt(A_0))
231  */
234  {
235  int nq = m_vessels[m_currentDomain*m_nVariables]->GetTotPoints();
236  NekDouble p = 0.0;
237  NekDouble p_t = 0.0;
238 
239  switch (i)
240  {
241  case 0: // Flux for A equation
242  {
243  for (int j = 0; j < nq; j++)
244  {
245  flux[0][j] = physfield[0][j]*physfield[1][j];
246  }
247  }
248  break;
249  case 1: // Flux for u equation
250  {
251  for (int j = 0; j < nq; j++)
252  {
253  ASSERTL0(physfield[0][j]>=0,"Negative A not allowed.");
254 
255  p = m_pext + m_beta[m_currentDomain][j]*
256  (sqrt(physfield[0][j]) - sqrt(m_A_0[m_currentDomain][j]));
257 
258  p_t = (physfield[1][j]*physfield[1][j])/2 + p/m_rho;
259  flux[0][j] = p_t;
260  }
261  }
262  break;
263  default:
264  ASSERTL0(false,"GetFluxVector: illegal vector index");
265  break;
266  }
267  }
268 
269 
270  /**
271  * Calculates the third term of the weak form (1): numerical flux
272  * at boundary \f$ \left[ \mathbf{\psi}^{\delta} \cdot \{
273  * \mathbf{F}^u - \mathbf{F}(\mathbf{U}^{\delta}) \}
274  * \right]_{x_e^l}^{x_eû} \f$
275  */
277  Array<OneD, Array<OneD, NekDouble> > &numflux)
278  {
279  int i;
280  int nTracePts = GetTraceTotPoints();
281 
284 
285  for (i = 0; i < m_nVariables; ++i)
286  {
287  Fwd[i] = Array<OneD, NekDouble>(nTracePts);
288  Bwd[i] = Array<OneD, NekDouble>(nTracePts);
289  }
290 
291  // Get the physical values at the trace
292  for (i = 0; i < m_nVariables; ++i)
293  {
294  m_vessels[m_currentDomain*m_nVariables+ i]->
295  GetFwdBwdTracePhys(physfield[i],Fwd[i],Bwd[i]);
296  }
297 
298  // Solve the upwinding Riemann problem within one arterial
299  // segment by calling the upwinding Riemann solver implemented
300  // in this file
301  NekDouble Aflux, uflux;
302  for (i = 0; i < nTracePts; ++i)
303  {
304  switch(m_upwindTypePulse)
305  {
306  case eUpwindPulse:
307  {
308  RiemannSolverUpwind(Fwd[0][i],Fwd[1][i],Bwd[0][i],Bwd[1][i],
309  Aflux, uflux, m_A_0_trace[m_currentDomain][i],
312  }
313  break;
314  default:
315  {
316  ASSERTL0(false,"populate switch statement for upwind flux");
317  }
318  break;
319  }
320  numflux[0][i] = Aflux;
321  numflux[1][i] = uflux;
322  }
323  }
324 
325 
326  /**
327  * Riemann solver for upwinding at an interface between two
328  * elements. Uses the characteristic variables for calculating
329  * the upwinded state \f$(A_u,u_u)\f$ from the left
330  * \f$(A_L,u_L)\f$ and right state \f$(A_R,u_R)\f$. Returns the
331  * upwinded flux $\mathbf{F}^u$ needed for the weak formulation
332  * (1). Details can be found in "Pulse wave propagation in the
333  * human vascular system", section 3.3
334  *
335  */
337  NekDouble AR,NekDouble uR,
338  NekDouble &Aflux,
339  NekDouble &uflux,
340  NekDouble A_0,
341  NekDouble beta,
342  NekDouble n)
343  {
345  Array<OneD, NekDouble> upwindedphysfield(2);
346  NekDouble cL = 0.0;
347  NekDouble cR = 0.0;
348  NekDouble rho = m_rho;
349  NekDouble pext = m_pext;
350  NekDouble p = 0.0;
351  NekDouble p_t = 0.0;
352 
353  // Compute the wave speeds. The use of the normal here allows
354  // for the definition of the characteristics to be inverted
355  // (and hence the left and right state) if n is in the -ve
356  // x-direction. This means we end up with the positive
357  // defintion of the flux which has to therefore be multiplied
358  // by the normal at the end of the methods This is a bit of a
359  // mind twister but is efficient from a coding perspective.
360  cL = sqrt(beta*sqrt(AL)/(2*rho))*n;
361  cR = sqrt(beta*sqrt(AR)/(2*rho))*n;
362 
363  ASSERTL1(fabs(cL+cR) > fabs(uL+uR),"Conditions are not sub-sonic");
364 
365  // If upwinding from left and right for subsonic domain
366  // then know characteristics immediately
367  W[0] = uL + 4*cL;
368  W[1] = uR - 4*cR;
369 
370  // Calculate conservative variables from characteristics
371  NekDouble w0mw1 = 0.25*(W[0]-W[1]);
372  NekDouble fac = rho/(2*beta);
373  w0mw1 *= w0mw1; // squared
374  w0mw1 *= w0mw1; // fourth power
375  fac *= fac; // squared
376  upwindedphysfield[0]= w0mw1*fac;
377  upwindedphysfield[1]= 0.5*(W[0] + W[1]);
378 
379  // Compute the fluxes multipled by the normal.
380  Aflux = upwindedphysfield[0] * upwindedphysfield[1]*n;
381  p = pext + beta*(sqrt(upwindedphysfield[0]) - sqrt(A_0));
382  p_t = 0.5*(upwindedphysfield[1]*upwindedphysfield[1]) + p/rho;
383  uflux = p_t*n;
384  }
385 
386  /**
387  * Print summary routine, calls virtual routine reimplemented in
388  * UnsteadySystem
389  */
391  {
393  }
394 }
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)
static std::string className
Name of class.
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
boost::shared_ptr< SessionReader > SessionReaderSharedPtr
Definition: MeshPartition.h:50
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)
PulseWavePropagation(const LibUtilities::SessionReaderSharedPtr &pSession)
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
static EquationSystemSharedPtr create(const LibUtilities::SessionReaderSharedPtr &pSession)
Creates an instance of this class.
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:1038
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