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
33// formulation (1):
34//
35///////////////////////////////////////////////////////////////////////////////
36
37#include <iostream>
38
41
42namespace Nektar
43{
44
47 "PulseWavePropagation", PulseWavePropagation::create,
48 "Pulse Wave Propagation equation.");
49/**
50 * @class PulseWavePropagation
51 *
52 * Set up the routines based on the weak formulation from
53 * "Computational Modelling of 1D blood flow with variable
54 * mechanical properties" by S. J. Sherwin et al. The weak
55 * formulation (1) reads:
56 * \f$ \sum_{e=1}^{N_{el}} \left[ \left( \frac{\partial \mathbf{U}^{\delta}
57 * }{\partial t} , \mathbf{\psi}^{\delta} \right)_{\Omega_e} - \left(
58 * \frac{\partial \mathbf{F(\mathbf{U})}^{\delta} }
59 * {\partial x}, \mathbf{\psi}^{\delta} \right)_{\Omega_e} + \left[
60 * \mathbf{\psi}^{\delta} \cdot \{ \mathbf{F}^u -
61 * \mathbf{F}(\mathbf{U}^{\delta}) \} \right]_{x_e^l}^{x_eû} \right] = 0 \f$
62 */
66 : PulseWaveSystem(pSession, pGraph)
67{
68}
69
70void PulseWavePropagation::v_InitObject([[maybe_unused]] bool DeclareField)
71{
72 // Will set up an array of vessels/fields in PulseWaveSystem::v_InitObject
73 // so set DeclareField to false so that the fields are not set up in
74 // EquationSystem unnecessarily. Note the number of fields in Equation
75 // system is related to the number of variables. The number of vessels is
76 // therefore held in PulwWaveSystem.
78
79 if (m_session->DefinesSolverInfo("PressureArea"))
80 {
82 m_session->GetSolverInfo("PressureArea"), m_vessels, m_session);
83 }
84 else
85 {
87 "Beta", m_vessels, m_session);
88 }
89
91 {
94 }
95 else
96 {
97 ASSERTL0(false, "Implicit Pulse Wave Propagation not set up.");
98 }
99
100 // Create advection object
101 std::string advName;
102 std::string riemName;
103 switch (m_upwindTypePulse)
104 {
105 case eUpwindPulse:
106 {
107 advName = "WeakDG";
108 riemName = "UpwindPulse";
109 }
110 break;
111 default:
112 {
113 ASSERTL0(false, "populate switch statement for upwind flux");
114 }
115 break;
116 }
121 riemName, m_session);
122 m_riemannSolver->SetScalar("A0", &PulseWavePropagation::GetA0, this);
123 m_riemannSolver->SetScalar("beta", &PulseWavePropagation::GetBeta, this);
124 m_riemannSolver->SetScalar("alpha", &PulseWavePropagation::GetAlpha, this);
125 m_riemannSolver->SetScalar("N", &PulseWavePropagation::GetN, this);
126 m_riemannSolver->SetParam("rho", &PulseWavePropagation::GetRho, this);
128 this);
129
130 m_advObject->SetRiemannSolver(m_riemannSolver);
131 m_advObject->InitObject(m_session, m_fields);
132}
133
134/**
135 * Computes the right hand side of (1). The RHS is everything
136 * except the term that contains the time derivative
137 * \f$\frac{\partial \mathbf{U}}{\partial t}\f$. In case of a
138 * Discontinuous Galerkin projection, m_advObject->Advect
139 * will be called
140 *
141 */
143 const Array<OneD, const Array<OneD, NekDouble>> &inarray,
144 Array<OneD, Array<OneD, NekDouble>> &outarray, const NekDouble time)
145{
146 size_t i;
147
149
150 // Dummy array for WeakDG advection
152
153 // Output array for advection
155
156 size_t cnt = 0;
157
158 // Set up Inflow and Outflow boundary conditions.
159 SetPulseWaveBoundaryConditions(inarray, outarray, time);
160
161 // Set up any interface conditions and write into boundary condition
163
164 // do advection evaluation in all domains
165 for (size_t omega = 0; omega < m_nDomains; ++omega)
166 {
168 m_currentDomain = omega;
169 size_t nq = m_vessels[omega * m_nVariables]->GetTotPoints();
170
171 timer.Start();
172 for (i = 0; i < m_nVariables; ++i)
173 {
174 physarray[i] = inarray[i] + cnt; // note this is doing a hidden copy
175 out[i] = outarray[i] + cnt;
176 }
177
178 for (i = 0; i < m_nVariables; ++i)
179 {
180 m_fields[i] = m_vessels[omega * m_nVariables + i];
181 }
182
183 m_advObject->Advect(m_nVariables, m_fields, advVel, physarray, out,
184 time);
185 for (i = 0; i < m_nVariables; ++i)
186 {
187 Vmath::Neg(nq, out[i], 1);
188 }
189 timer.Stop();
190 timer.AccumulateRegion("PulseWavePropagation:_DoOdeRHS", 1);
191 cnt += nq;
192 }
193}
194
196 const Array<OneD, const Array<OneD, NekDouble>> &inarray,
198 [[maybe_unused]] const NekDouble time)
199{
200 // Just copy over array
201 if (inarray != outarray)
202 {
203 for (size_t i = 0; i < m_nVariables; ++i)
204 {
205 Vmath::Vcopy(inarray[i].size(), inarray[i], 1, outarray[i], 1);
206 }
207 }
208}
209
210/**
211 * Does the projection between ... space and the ... space. Also checks for
212 *Q-inflow boundary conditions at the inflow of the current arterial segment and
213 *applies the Q-inflow if specified
214 */
216 const Array<OneD, const Array<OneD, NekDouble>> &inarray,
217 [[maybe_unused]] Array<OneD, Array<OneD, NekDouble>> &outarray,
218 const NekDouble time)
219
220{
221 size_t omega;
222
224
225 size_t offset = 0;
226
227 // This will be moved to the RCR boundary condition once factory is setup
228 if (time == 0)
229 {
231
232 for (omega = 0; omega < m_nDomains; ++omega)
233 {
234 vessel[0] = m_vessels[2 * omega];
235 vessel[1] = m_vessels[2 * omega + 1];
236
237 for (size_t j = 0; j < 2; ++j)
238 {
239 std::string BCType;
240
241 if (j < vessel[0]->GetBndConditions().size())
242 {
243 BCType = vessel[0]->GetBndConditions()[j]->GetUserDefined();
244 }
245
246 // if no condition given define it to be NoUserDefined
247 if (BCType.empty() || BCType == "Interface")
248 {
249 BCType = "NoUserDefined";
250 }
251
254
255 // turn on time dependent BCs
256 if (BCType == "Q-inflow")
257 {
258 vessel[0]->GetBndConditions()[j]->SetIsTimeDependent(true);
259 }
260 else if (BCType == "A-inflow")
261 {
262 vessel[0]->GetBndConditions()[j]->SetIsTimeDependent(true);
263 }
264 else if (BCType == "U-inflow")
265 {
266 vessel[1]->GetBndConditions()[j]->SetIsTimeDependent(true);
267 }
268 else if (BCType == "RCR-terminal")
269 {
270 vessel[0]->GetBndConditions()[j]->SetIsTimeDependent(true);
271 }
272 }
273 }
274 }
275
277
278 // Loop over all vessels and set boundary conditions
280 for (omega = 0; omega < m_nDomains; ++omega)
281 {
282 timer.Start();
283 for (size_t n = 0; n < 2; ++n)
284 {
285 m_Boundary[2 * omega + n]->DoBoundary(
286 inarray, m_A_0, m_beta, m_alpha, time, omega, offset, n);
287 }
288
289 offset += m_vessels[2 * omega]->GetTotPoints();
290 timer.Stop();
291 timer.AccumulateRegion("PulseWavePropagation:_SetBCs", 1);
292 }
293}
294
295/**
296 * Calculates the second term of the weak form (1): \f$
297 * \left( \frac{\partial \mathbf{F(\mathbf{U})}^{\delta}
298 * }{\partial x}, \mathbf{\psi}^{\delta} \right)_{\Omega_e}
299 * \f$
300 * The variables of the system are $\mathbf{U} = [A,u]^T$
301 * physfield[0] = A physfield[1] = u
302 * flux[0] = F[0] = A*u flux[1] = F[1] = u^2/2 + p/rho
303 */
305 const Array<OneD, Array<OneD, NekDouble>> &physfield,
307{
308 size_t nq = m_vessels[m_currentDomain * m_nVariables]->GetTotPoints();
309 NekDouble domain = m_currentDomain;
311 Array<OneD, NekDouble> dAUdx(nq);
312 NekDouble viscoelasticGradient = 0.0;
313
315
316 timer.Start();
317 for (size_t j = 0; j < nq; ++j)
318 {
319 flux[0][0][j] = physfield[0][j] * physfield[1][j];
320 }
321 timer.Stop();
322 timer.AccumulateRegion("PulseWavePropagation:GetFluxVector-flux", 3);
323
324 // d/dx of AU, for the viscoelastic tube law and extra fields
325 m_fields[0]->PhysDeriv(flux[0][0], dAUdx);
326
327 for (size_t j = 0; j < nq; ++j)
328 {
329 if ((j == 0) || (j == nq - 1))
330 {
331 viscoelasticGradient = dAUdx[j];
332 }
333 else
334 {
335 viscoelasticGradient = (dAUdx[j] + dAUdx[j + 1]) / 2;
336 }
337
338 m_pressureArea->GetPressure(m_pressure[domain][j], m_beta[domain][j],
339 physfield[0][j], m_A_0[domain][j],
340 viscoelasticGradient, m_gamma[domain][j],
341 m_alpha[domain][j]);
342
343 flux[1][0][j] = physfield[1][j] * physfield[1][j] / 2 +
344 m_pressure[domain][j] / m_rho;
345 }
346
347 m_session->MatchSolverInfo("OutputExtraFields", "True", extraFields, true);
348
349 if (extraFields)
350 {
351 /*
352 Calculates wave speed and characteristic variables.
353
354 Ideally this should be moved to PulseWaveSystem, but it's easiest to
355 implement here.
356 */
357 size_t counter = 0;
358
359 m_PWV[domain] = Array<OneD, NekDouble>(nq);
360 m_W1[domain] = Array<OneD, NekDouble>(nq);
361 m_W2[domain] = Array<OneD, NekDouble>(nq);
362
363 for (size_t j = 0; j < nq; ++j)
364 {
365 m_pressureArea->GetC(m_PWV[domain][j], m_beta[domain][j],
366 physfield[0][counter + j], m_A_0[domain][j],
367 m_alpha[domain][j]);
368 m_pressureArea->GetW1(m_W1[domain][j], physfield[1][counter + j],
369 m_beta[domain][j], physfield[0][counter + j],
370 m_A_0[domain][j], m_alpha[domain][j]);
371 m_pressureArea->GetW2(m_W2[domain][j], physfield[1][counter + j],
372 m_beta[domain][j], physfield[0][counter + j],
373 m_A_0[domain][j], m_alpha[domain][j]);
374 }
375
376 counter += nq;
377 }
378}
379
381{
383}
384
386{
388}
389
391{
393}
394
396{
398}
399
401{
402 return m_rho;
403}
404
406{
407 return m_nDomains;
408}
409
410/**
411 * Print summary routine, calls virtual routine reimplemented in
412 * UnsteadySystem
413 */
415{
417}
418
419} // namespace Nektar
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:208
tKey RegisterCreatorFunction(tKey idKey, CreatorFunction classCreator, std::string pDesc="")
Register a class with the factory.
tBaseSharedPtr CreateInstance(tKey idKey, tParam... args)
Create an instance of the class referred to by idKey.
void DefineProjection(FuncPointerT func, ObjectPointerT obj)
void DefineOdeRhs(FuncPointerT func, ObjectPointerT obj)
void AccumulateRegion(std::string, int iolevel=0)
Accumulate elapsed time for a region.
Definition: Timer.cpp:70
Array< OneD, NekDouble > & GetAlpha()
void v_InitObject(bool DeclareField=false) override
Init object for UnsteadySystem class.
SolverUtils::RiemannSolverSharedPtr m_riemannSolver
PulseWavePropagation(const LibUtilities::SessionReaderSharedPtr &pSession, const SpatialDomains::MeshGraphSharedPtr &pGraph)
SolverUtils::AdvectionSharedPtr m_advObject
void SetPulseWaveBoundaryConditions(const Array< OneD, const Array< OneD, NekDouble > > &inarray, Array< OneD, Array< OneD, NekDouble > > &outarray, const NekDouble time)
void DoOdeProjection(const Array< OneD, const Array< OneD, NekDouble > > &inarray, Array< OneD, Array< OneD, NekDouble > > &outarray, const NekDouble time)
Array< OneD, PulseWaveBoundarySharedPtr > m_Boundary
Array< OneD, NekDouble > & GetN()
Array< OneD, NekDouble > & GetA0()
void GetFluxVector(const Array< OneD, Array< OneD, NekDouble > > &physfield, Array< OneD, Array< OneD, Array< OneD, NekDouble > > > &flux)
DG Pulse Wave Propagation routines:
void v_GenerateSummary(SolverUtils::SummaryList &s) override
void DoOdeRhs(const Array< OneD, const Array< OneD, NekDouble > > &inarray, Array< OneD, Array< OneD, NekDouble > > &outarray, const NekDouble time)
Array< OneD, NekDouble > & GetBeta()
static EquationSystemSharedPtr create(const LibUtilities::SessionReaderSharedPtr &pSession, const SpatialDomains::MeshGraphSharedPtr &pGraph)
Creates an instance of this class.
static std::string className
Name of class.
Base class for unsteady solvers.
Array< OneD, Array< OneD, NekDouble > > m_A_0
Array< OneD, Array< OneD, NekDouble > > m_beta_trace
Array< OneD, Array< OneD, NekDouble > > m_W2
Array< OneD, Array< OneD, NekDouble > > m_trace_fwd_normal
PulseWavePressureAreaSharedPtr m_pressureArea
Array< OneD, Array< OneD, NekDouble > > m_W1
void v_InitObject(bool DeclareField=false) override
Array< OneD, Array< OneD, NekDouble > > m_alpha_trace
UpwindTypePulse m_upwindTypePulse
Array< OneD, Array< OneD, NekDouble > > m_A_0_trace
Array< OneD, Array< OneD, NekDouble > > m_pressure
Array< OneD, Array< OneD, NekDouble > > m_gamma
Array< OneD, Array< OneD, NekDouble > > m_alpha
Array< OneD, Array< OneD, NekDouble > > m_PWV
Array< OneD, MultiRegions::ExpListSharedPtr > m_vessels
void EnforceInterfaceConditions(const Array< OneD, const Array< OneD, NekDouble > > &fields)
Array< OneD, Array< OneD, NekDouble > > m_beta
int m_spacedim
Spatial dimension (>= expansion dim).
Array< OneD, MultiRegions::ExpListSharedPtr > m_fields
Array holding all dependent variables.
LibUtilities::SessionReaderSharedPtr m_session
The session reader.
SOLVER_UTILS_EXPORT void SetBoundaryConditions(NekDouble time)
Evaluates the boundary conditions at the given time.
LibUtilities::TimeIntegrationSchemeOperators m_ode
The time integration scheme operators to use.
bool m_explicitAdvection
Indicates if explicit or implicit treatment of advection is used.
SOLVER_UTILS_EXPORT void v_GenerateSummary(SummaryList &s) override
Print a summary of time stepping parameters.
std::shared_ptr< SessionReader > SessionReaderSharedPtr
AdvectionFactory & GetAdvectionFactory()
Gets the factory for initialising advection objects.
Definition: Advection.cpp:43
std::vector< std::pair< std::string, std::string > > SummaryList
Definition: Misc.h:46
EquationSystemFactory & GetEquationSystemFactory()
RiemannSolverFactory & GetRiemannSolverFactory()
std::shared_ptr< MeshGraph > MeshGraphSharedPtr
Definition: MeshGraph.h:174
PressureAreaFactory & GetPressureAreaFactory()
@ eUpwindPulse
simple upwinding scheme
BoundaryFactory & GetBoundaryFactory()
double NekDouble
void Neg(int n, T *x, const int incx)
Negate x = -x.
Definition: Vmath.hpp:292
void Vcopy(int n, const T *x, const int incx, T *y, const int incy)
Definition: Vmath.hpp:825