From 87459379cc58f5ef4f134506d23c1639b6e9501e Mon Sep 17 00:00:00 2001 From: "Documenter.jl" Date: Mon, 13 Jan 2025 12:29:55 +0000 Subject: [PATCH] build based on 0b07dff --- dev/.documenter-siteinfo.json | 2 +- dev/index.html | 8 ++++---- dev/math/index.html | 2 +- dev/search_index.js | 2 +- dev/tuto/index.html | 2 +- 5 files changed, 8 insertions(+), 8 deletions(-) diff --git a/dev/.documenter-siteinfo.json b/dev/.documenter-siteinfo.json index 2395046..9243052 100644 --- a/dev/.documenter-siteinfo.json +++ b/dev/.documenter-siteinfo.json @@ -1 +1 @@ -{"documenter":{"julia_version":"1.11.2","generation_timestamp":"2025-01-13T12:26:28","documenter_version":"1.8.0"}} \ No newline at end of file +{"documenter":{"julia_version":"1.11.2","generation_timestamp":"2025-01-13T12:29:50","documenter_version":"1.8.0"}} \ No newline at end of file diff --git a/dev/index.html b/dev/index.html index 3fc8918..a741c53 100644 --- a/dev/index.html +++ b/dev/index.html @@ -12,18 +12,18 @@ R &= \sqrt{2A} \end{aligned}\right.\]

Both models were designed to be used with Trixi.jl, a flexible and high-performance framework for solving systems of conservation laws using the Discontinuous Galerkin (DG) method.

Features

Installation

To install BloodFlowTrixi.jl, use the following commands in Julia:

julia> ]
 pkg> add Trixi
-pkg> add https://github.com/your-repo/BloodFlowTrixi.jl

Future Plans

short term

long term

License

This package is licensed under the MIT license.

Acknowledgments

This package was developed as part of my PhD research in applied mathematics, focusing on mathematical modeling and numerical simulation of blood flow in arteries. Special thanks to the developers of Trixi.jl, whose framework was invaluable in implementing and testing these models.

BloodFlowTrixi.BloodFlowTrixi — Module
Package BloodFlowTrixi v0.0.1

This package implements 1D and 2D blood flow models for arterial circulation using Trixi.jl, enabling efficient numerical simulation and analysis.

Docs under https://yolhan83.github.io/BloodFlowTrixi.jl

source
BloodFlowTrixi.BloodFlowEquations1D — Type
BloodFlowEquations1D(;h,rho=1.0,xi=0.25,nu=0.04)

Blood Flow equations in one space dimension. This model describes the dynamics of blood flow along a compliant artery using one-dimensional equations derived from the Navier-Stokes equations. The equations account for conservation of mass and momentum, incorporating the effect of arterial compliance and frictional losses.

The governing equations are given by

\[\left\{\begin{aligned} +pkg> add https://github.com/your-repo/BloodFlowTrixi.jl

Future Plans

short term

  • Add second order 1D model.
  • Design prim variables for 1D and 2D models.
  • Add proper tests for 1D and 2D models.
  • Add 3D representations of the solutions for 1D and 2D models.
  • Design easy to use interfaces for users to define their own initial and boundary conditions and source terms.

long term

  • Add 3D fluid-structure interaction models for complex arterial geometries.
  • Design support for artery networks and simulate vascular networks using the 2D and 1D model.
  • Autodiff support for 1D and 2D models for parameter optimization.

License

This package is licensed under the MIT license.

Acknowledgments

This package was developed as part of my PhD research in applied mathematics, focusing on mathematical modeling and numerical simulation of blood flow in arteries. Special thanks to the developers of Trixi.jl, whose framework was invaluable in implementing and testing these models.

BloodFlowTrixi.BloodFlowTrixi — Module
Package BloodFlowTrixi v0.0.1

This package implements 1D and 2D blood flow models for arterial circulation using Trixi.jl, enabling efficient numerical simulation and analysis.

Docs under https://yolhan83.github.io/BloodFlowTrixi.jl

source
BloodFlowTrixi.BloodFlowEquations1D — Type
BloodFlowEquations1D(;h,rho=1.0,xi=0.25,nu=0.04)

Blood Flow equations in one space dimension. This model describes the dynamics of blood flow along a compliant artery using one-dimensional equations derived from the Navier-Stokes equations. The equations account for conservation of mass and momentum, incorporating the effect of arterial compliance and frictional losses.

The governing equations are given by

\[\left\{\begin{aligned} \frac{\partial a}{\partial t} + \frac{\partial}{\partial x}(Q) &= 0 \\ \frac{\partial Q}{\partial t} + \frac{\partial}{\partial x}\left(\frac{Q^2}{A} + A P(a)\right) &= P(a) \frac{\partial A}{\partial x} - 2 \pi R k \frac Q {A}\\ P(a) &= P_{ext} + \frac{Eh\sqrt{\pi}}{1-\xi^2}\frac{\sqrt{A} - \sqrt{A_0}}{A_0} \\ R &= \sqrt{\frac{A}{\pi}} -\end{aligned}\right.\]

source
BloodFlowTrixi.BloodFlowEquations2D — Type
BloodFlowEquations2D(;h,rho=1.0,xi=0.25)

Defines the two-dimensional blood flow equations derived from the Navier-Stokes equations in curvilinear coordinates under the thin-artery assumption. This model describes the dynamics of blood flow along a compliant artery in two spatial dimensions (s, θ).

Parameters

  • h::T: Wall thickness of the artery.
  • rho::T: Fluid density (default 1.0).
  • xi::T: Poisson's ratio (default 0.25).
  • nu::T: Viscosity coefficient.

The governing equations account for conservation of mass and momentum, incorporating the effects of arterial compliance, curvature, and frictional losses.

\[\left\{\begin{aligned} +\end{aligned}\right.\]

source
BloodFlowTrixi.BloodFlowEquations2D — Type
BloodFlowEquations2D(;h,rho=1.0,xi=0.25)

Defines the two-dimensional blood flow equations derived from the Navier-Stokes equations in curvilinear coordinates under the thin-artery assumption. This model describes the dynamics of blood flow along a compliant artery in two spatial dimensions (s, θ).

Parameters

  • h::T: Wall thickness of the artery.
  • rho::T: Fluid density (default 1.0).
  • xi::T: Poisson's ratio (default 0.25).
  • nu::T: Viscosity coefficient.

The governing equations account for conservation of mass and momentum, incorporating the effects of arterial compliance, curvature, and frictional losses.

\[\left\{\begin{aligned} \frac{\partial a}{\partial t} + \frac{\partial}{\partial \theta}\left( \frac{Q_{R\theta}}{A} \right) + \frac{\partial}{\partial s}(Q_s) &= 0 \\ \frac{\partial Q_{R\theta}}{\partial t} + \frac{\partial}{\partial \theta}\left(\frac{Q_{R\theta}^2}{2A^2} + A P(a)\right) + \frac{\partial}{\partial s}\left( \frac{Q_{R\theta}Q_s}{A} \right) &= P(a) \frac{\partial A}{\partial \theta} - 2 R k \frac{Q_{R\theta}}{A} + \frac{2R}{3} \mathcal{C}\sin \theta \frac{Q_s^2}{A} \\ \frac{\partial Q_{s}}{\partial t} + \frac{\partial}{\partial \theta}\left(\frac{Q_{R\theta} Q_s}{A^2} \right) + \frac{\partial}{\partial s}\left( \frac{Q_s^2}{A} - \frac{Q_{R\theta}^2}{2A^2} + A P(a) \right) &= P(a) \frac{\partial A}{\partial s} - R k \frac{Q_s}{A} - \frac{2R}{3} \mathcal{C}\sin \theta \frac{Q_s Q_{R\theta}}{A^2} \\ P(a) &= P_{ext} + \frac{Eh}{\sqrt{2}\left(1-\xi^2\right)}\frac{\sqrt{A} - \sqrt{A_0}}{A_0} \\ R &= \sqrt{2A} -\end{aligned}\right.\]

source
Trixi.DissipationLocalLaxFriedrichs — Method
(dissipation::Trixi.DissipationLocalLaxFriedrichs)(u_ll, u_rr, orientation_or_normal_direction, eq::BloodFlowEquations1D)

Calculates the dissipation term using the Local Lax-Friedrichs method.

Parameters

  • u_ll: Left state vector.
  • u_rr: Right state vector.
  • orientation_or_normal_direction: Orientation or normal direction.
  • eq: Instance of BloodFlowEquations1D.

Returns

Dissipation vector.

source
BloodFlowTrixi.boundary_condition_outflow — Method
boundary_condition_outflow(u_inner, orientation_or_normal, direction, x, t, surface_flux_function, eq::BloodFlowEquations1D)

Implements the outflow boundary condition, assuming that there is no reflection at the boundary.

Parameters

  • u_inner: State vector inside the domain near the boundary.
  • orientation_or_normal: Normal orientation of the boundary.
  • direction: Integer indicating the direction of the boundary.
  • x: Position vector.
  • t: Time.
  • surface_flux_function: Function to compute flux at the boundary.
  • eq: Instance of BloodFlowEquations1D.

Returns

Computed boundary flux.

source
BloodFlowTrixi.boundary_condition_pressure_in — Method
boundary_condition_pressure_in(u_inner, orientation_or_normal, direction, x, t, surface_flux_function, eq::BloodFlowEquations1D)

Implements a pressure inflow boundary condition where the inflow pressure varies with time.

Parameters

  • u_inner: State vector inside the domain near the boundary.
  • orientation_or_normal: Normal orientation of the boundary.
  • direction: Integer indicating the boundary direction.
  • x: Position vector.
  • t: Time scalar.
  • surface_flux_function: Function to compute flux at the boundary.
  • eq: Instance of BloodFlowEquations1D.

Returns

Computed boundary flux with inflow pressure specified by:

\[P_{in} = \begin{cases} +\end{aligned}\right.\]

source
Trixi.DissipationLocalLaxFriedrichs — Method
(dissipation::Trixi.DissipationLocalLaxFriedrichs)(u_ll, u_rr, orientation_or_normal_direction, eq::BloodFlowEquations1D)

Calculates the dissipation term using the Local Lax-Friedrichs method.

Parameters

  • u_ll: Left state vector.
  • u_rr: Right state vector.
  • orientation_or_normal_direction: Orientation or normal direction.
  • eq: Instance of BloodFlowEquations1D.

Returns

Dissipation vector.

source
BloodFlowTrixi.boundary_condition_outflow — Method
boundary_condition_outflow(u_inner, orientation_or_normal, direction, x, t, surface_flux_function, eq::BloodFlowEquations1D)

Implements the outflow boundary condition, assuming that there is no reflection at the boundary.

Parameters

  • u_inner: State vector inside the domain near the boundary.
  • orientation_or_normal: Normal orientation of the boundary.
  • direction: Integer indicating the direction of the boundary.
  • x: Position vector.
  • t: Time.
  • surface_flux_function: Function to compute flux at the boundary.
  • eq: Instance of BloodFlowEquations1D.

Returns

Computed boundary flux.

source
BloodFlowTrixi.boundary_condition_pressure_in — Method
boundary_condition_pressure_in(u_inner, orientation_or_normal, direction, x, t, surface_flux_function, eq::BloodFlowEquations1D)

Implements a pressure inflow boundary condition where the inflow pressure varies with time.

Parameters

  • u_inner: State vector inside the domain near the boundary.
  • orientation_or_normal: Normal orientation of the boundary.
  • direction: Integer indicating the boundary direction.
  • x: Position vector.
  • t: Time scalar.
  • surface_flux_function: Function to compute flux at the boundary.
  • eq: Instance of BloodFlowEquations1D.

Returns

Computed boundary flux with inflow pressure specified by:

\[P_{in} = \begin{cases} 2 \times 10^4 \sin^2(\pi t / 0.125) & \text{if } t < 0.125 \\ 0 & \text{otherwise} -\end{cases}\]

The corresponding inflow area A_{in} is computed using the inverse pressure relation, and the boundary state is constructed accordingly.

source
BloodFlowTrixi.boundary_condition_slip_wall — Method
boundary_condition_slip_wall(u_inner, orientation_or_normal, direction, x, t, surface_flux_function, eq::BloodFlowEquations1D)

Implements a slip wall boundary condition where the normal component of velocity is reflected.

Parameters

  • u_inner: State vector inside the domain near the boundary.
  • orientation_or_normal: Normal orientation of the boundary.
  • direction: Integer indicating the direction of the boundary.
  • x: Position vector.
  • t: Time.
  • surface_flux_function: Function to compute flux at the boundary.
  • eq: Instance of BloodFlowEquations1D.

Returns

Computed boundary flux at the slip wall.

source
BloodFlowTrixi.flux_nonconservative — Method
flux_nonconservative(u_ll, u_rr, orientation::Integer, eq::BloodFlowEquations1D)

Computes the non-conservative flux for the model, used for handling discontinuities in pressure.

Parameters

  • u_ll: Left state vector.
  • u_rr: Right state vector.
  • orientation::Integer: Orientation index.
  • eq: Instance of BloodFlowEquations1D.

Returns

Non-conservative flux vector.

source
BloodFlowTrixi.friction — Method
friction(u, x, eq::BloodFlowEquations1D)

Calculates the friction term for the blood flow equations, which represents viscous resistance to flow along the artery wall.

Parameters

  • u: State vector containing cross-sectional area and flow rate.
  • x: Position along the artery.
  • eq::BloodFlowEquations1D: Instance of the blood flow model.

Returns

Friction coefficient as a scalar.

source
BloodFlowTrixi.initial_condition_simple — Method
initial_condition_simple(x, t, eq::BloodFlowEquations1D; R0=2.0)

Generates a simple initial condition with a specified initial radius R0.

Parameters

  • x: Position vector.
  • t: Time scalar.
  • eq::BloodFlowEquations1D: Instance of the blood flow model.
  • R0: Initial radius (default: 2.0).

Returns

State vector with zero initial area perturbation, zero flow rate, constant elasticity modulus, and reference area computed as A_0 = \pi R_0^2.

This initial condition is suitable for basic tests without complex dynamics.

source
BloodFlowTrixi.inv_pressure — Method
inv_pressure(p, u, eq::BloodFlowEquations1D)

Computes the inverse relation of pressure to cross-sectional area.

Parameters

  • p: Pressure.
  • u: State vector.
  • eq: Instance of BloodFlowEquations1D.

Returns

Cross-sectional area corresponding to the given pressure.

source
BloodFlowTrixi.pressure — Method
pressure(u, eq::BloodFlowEquations1D)

Computes the pressure given the state vector based on the compliance of the artery.

Parameters

  • u: State vector.
  • eq: Instance of BloodFlowEquations1D.

Returns

Pressure as a scalar.

source
BloodFlowTrixi.pressure_der — Method
pressure_der(u, eq::BloodFlowEquations1D)

Computes the derivative of pressure with respect to cross-sectional area.

Parameters

  • u: State vector.
  • eq: Instance of BloodFlowEquations1D.

Returns

Derivative of pressure.

source
BloodFlowTrixi.radius — Method
radius(u, eq::BloodFlowEquations1D)

Computes the radius of the artery based on the cross-sectional area.

Parameters

  • u: State vector.
  • eq: Instance of BloodFlowEquations1D.

Returns

Radius as a scalar.

source
BloodFlowTrixi.source_term_simple — Method
source_term_simple(u, x, t, eq::BloodFlowEquations1D)

Computes a simple source term for the blood flow model, focusing on frictional effects.

Parameters

  • u: State vector containing area perturbation, flow rate, elasticity modulus, and reference area.
  • x: Position vector.
  • t: Time scalar.
  • eq::BloodFlowEquations1D: Instance of the blood flow model.

Returns

Source terms vector where:

  • s_1 = 0 (no source for area perturbation).
  • s_2 represents the friction term given by s_2 = \frac{2 \pi k Q}{R A}.

Friction coefficient k is computed using the friction function, and the radius R is obtained using the radius function.

source
Trixi.cons2prim — Method
Trixi.cons2prim(u, eq::BloodFlowEquations1D)

Converts the conserved variables to primitive variables.

Parameters

  • u: State vector.
  • eq: Instance of BloodFlowEquations1D.

Returns

Primitive variable vector.

source
Trixi.flux — Method
Trixi.flux(u, orientation::Integer, eq::BloodFlowEquations1D)

Computes the flux vector for the conservation laws of the blood flow model.

Parameters

  • u: State vector.
  • orientation::Integer: Orientation index for flux computation.
  • eq: Instance of BloodFlowEquations1D.

Returns

Flux vector as an SVector.

source
Trixi.initial_condition_convergence_test — Method
initial_condition_convergence_test(x, t, eq::BloodFlowEquations1D)

Generates a smooth initial condition for convergence tests of the blood flow equations.

Parameters

  • x: Position vector.
  • t: Time scalar.
  • eq::BloodFlowEquations1D: Instance of the blood flow model.

Returns

Initial condition state vector with zero initial area perturbation, sinusoidal flow rate, a constant elasticity modulus, and reference area.

Details

The returned initial condition has:

  • Zero perturbation in area (a = 0).
  • A sinusoidal flow rate given by Q = sin(\pi x t).
  • A constant elasticity modulus E.
  • A reference cross-sectional area A_0 = \pi R_0^2 for R_0 = 1.

This initial condition can be used to verify the accuracy and stability of numerical solvers.

source
Trixi.max_abs_speed_naive — Method
Trixi.max_abs_speed_naive(u_ll, u_rr, orientation::Integer, eq::BloodFlowEquations1D)

Calculates the maximum absolute speed for wave propagation in the blood flow model using a naive approach.

Parameters

  • u_ll: Left state vector.
  • u_rr: Right state vector.
  • orientation::Integer: Orientation index.
  • eq: Instance of BloodFlowEquations1D.

Returns

Maximum absolute speed.

source
Trixi.prim2cons — Method
Trixi.prim2cons(u, eq::BloodFlowEquations1D)

Converts the primitive variables to conserved variables.

Parameters

  • u: Primitive variable vector.
  • eq: Instance of BloodFlowEquations1D.

Returns

Conserved variable vector.

source
Trixi.source_terms_convergence_test — Method
source_terms_convergence_test(u, x, t, eq::BloodFlowEquations1D)

Computes the source terms for convergence tests of the blood flow equations.

Parameters

  • u: State vector containing area perturbation, flow rate, elasticity modulus, and reference area.
  • x: Position vector.
  • t: Time scalar.
  • eq::BloodFlowEquations1D: Instance of the blood flow model.

Returns

Source terms vector.

Details

The source terms are derived based on the smooth initial condition and friction effects:

  • s_1 represents the source term for area perturbation and is given by s_1 = \pi t \cos(\pi x t).
  • s_2 represents the source term for the flow rate and includes contributions from spatial and temporal variations as well as friction effects.

The radius R is computed using the radius function, and the friction coefficient k is obtained using the friction function.

This function is useful for evaluating the correctness of source term handling in numerical solvers.

source
+\end{cases}\]

The corresponding inflow area A_{in} is computed using the inverse pressure relation, and the boundary state is constructed accordingly.

source
BloodFlowTrixi.boundary_condition_slip_wall — Method
boundary_condition_slip_wall(u_inner, orientation_or_normal, direction, x, t, surface_flux_function, eq::BloodFlowEquations1D)

Implements a slip wall boundary condition where the normal component of velocity is reflected.

Parameters

  • u_inner: State vector inside the domain near the boundary.
  • orientation_or_normal: Normal orientation of the boundary.
  • direction: Integer indicating the direction of the boundary.
  • x: Position vector.
  • t: Time.
  • surface_flux_function: Function to compute flux at the boundary.
  • eq: Instance of BloodFlowEquations1D.

Returns

Computed boundary flux at the slip wall.

source
BloodFlowTrixi.flux_nonconservative — Method
flux_nonconservative(u_ll, u_rr, orientation::Integer, eq::BloodFlowEquations1D)

Computes the non-conservative flux for the model, used for handling discontinuities in pressure.

Parameters

  • u_ll: Left state vector.
  • u_rr: Right state vector.
  • orientation::Integer: Orientation index.
  • eq: Instance of BloodFlowEquations1D.

Returns

Non-conservative flux vector.

source
BloodFlowTrixi.friction — Method
friction(u, x, eq::BloodFlowEquations1D)

Calculates the friction term for the blood flow equations, which represents viscous resistance to flow along the artery wall.

Parameters

  • u: State vector containing cross-sectional area and flow rate.
  • x: Position along the artery.
  • eq::BloodFlowEquations1D: Instance of the blood flow model.

Returns

Friction coefficient as a scalar.

source
BloodFlowTrixi.initial_condition_simple — Method
initial_condition_simple(x, t, eq::BloodFlowEquations1D; R0=2.0)

Generates a simple initial condition with a specified initial radius R0.

Parameters

  • x: Position vector.
  • t: Time scalar.
  • eq::BloodFlowEquations1D: Instance of the blood flow model.
  • R0: Initial radius (default: 2.0).

Returns

State vector with zero initial area perturbation, zero flow rate, constant elasticity modulus, and reference area computed as A_0 = \pi R_0^2.

This initial condition is suitable for basic tests without complex dynamics.

source
BloodFlowTrixi.inv_pressure — Method
inv_pressure(p, u, eq::BloodFlowEquations1D)

Computes the inverse relation of pressure to cross-sectional area.

Parameters

  • p: Pressure.
  • u: State vector.
  • eq: Instance of BloodFlowEquations1D.

Returns

Cross-sectional area corresponding to the given pressure.

source
BloodFlowTrixi.pressure — Method
pressure(u, eq::BloodFlowEquations1D)

Computes the pressure given the state vector based on the compliance of the artery.

Parameters

  • u: State vector.
  • eq: Instance of BloodFlowEquations1D.

Returns

Pressure as a scalar.

source
BloodFlowTrixi.pressure_der — Method
pressure_der(u, eq::BloodFlowEquations1D)

Computes the derivative of pressure with respect to cross-sectional area.

Parameters

  • u: State vector.
  • eq: Instance of BloodFlowEquations1D.

Returns

Derivative of pressure.

source
BloodFlowTrixi.radius — Method
radius(u, eq::BloodFlowEquations1D)

Computes the radius of the artery based on the cross-sectional area.

Parameters

  • u: State vector.
  • eq: Instance of BloodFlowEquations1D.

Returns

Radius as a scalar.

source
BloodFlowTrixi.source_term_simple — Method
source_term_simple(u, x, t, eq::BloodFlowEquations1D)

Computes a simple source term for the blood flow model, focusing on frictional effects.

Parameters

  • u: State vector containing area perturbation, flow rate, elasticity modulus, and reference area.
  • x: Position vector.
  • t: Time scalar.
  • eq::BloodFlowEquations1D: Instance of the blood flow model.

Returns

Source terms vector where:

  • s_1 = 0 (no source for area perturbation).
  • s_2 represents the friction term given by s_2 = \frac{2 \pi k Q}{R A}.

Friction coefficient k is computed using the friction function, and the radius R is obtained using the radius function.

source
Trixi.cons2prim — Method
Trixi.cons2prim(u, eq::BloodFlowEquations1D)

Converts the conserved variables to primitive variables.

Parameters

  • u: State vector.
  • eq: Instance of BloodFlowEquations1D.

Returns

Primitive variable vector.

source
Trixi.flux — Method
Trixi.flux(u, orientation::Integer, eq::BloodFlowEquations1D)

Computes the flux vector for the conservation laws of the blood flow model.

Parameters

  • u: State vector.
  • orientation::Integer: Orientation index for flux computation.
  • eq: Instance of BloodFlowEquations1D.

Returns

Flux vector as an SVector.

source
Trixi.initial_condition_convergence_test — Method
initial_condition_convergence_test(x, t, eq::BloodFlowEquations1D)

Generates a smooth initial condition for convergence tests of the blood flow equations.

Parameters

  • x: Position vector.
  • t: Time scalar.
  • eq::BloodFlowEquations1D: Instance of the blood flow model.

Returns

Initial condition state vector with zero initial area perturbation, sinusoidal flow rate, a constant elasticity modulus, and reference area.

Details

The returned initial condition has:

  • Zero perturbation in area (a = 0).
  • A sinusoidal flow rate given by Q = sin(\pi x t).
  • A constant elasticity modulus E.
  • A reference cross-sectional area A_0 = \pi R_0^2 for R_0 = 1.

This initial condition can be used to verify the accuracy and stability of numerical solvers.

source
Trixi.max_abs_speed_naive — Method
Trixi.max_abs_speed_naive(u_ll, u_rr, orientation::Integer, eq::BloodFlowEquations1D)

Calculates the maximum absolute speed for wave propagation in the blood flow model using a naive approach.

Parameters

  • u_ll: Left state vector.
  • u_rr: Right state vector.
  • orientation::Integer: Orientation index.
  • eq: Instance of BloodFlowEquations1D.

Returns

Maximum absolute speed.

source
Trixi.prim2cons — Method
Trixi.prim2cons(u, eq::BloodFlowEquations1D)

Converts the primitive variables to conserved variables.

Parameters

  • u: Primitive variable vector.
  • eq: Instance of BloodFlowEquations1D.

Returns

Conserved variable vector.

source
Trixi.source_terms_convergence_test — Method
source_terms_convergence_test(u, x, t, eq::BloodFlowEquations1D)

Computes the source terms for convergence tests of the blood flow equations.

Parameters

  • u: State vector containing area perturbation, flow rate, elasticity modulus, and reference area.
  • x: Position vector.
  • t: Time scalar.
  • eq::BloodFlowEquations1D: Instance of the blood flow model.

Returns

Source terms vector.

Details

The source terms are derived based on the smooth initial condition and friction effects:

  • s_1 represents the source term for area perturbation and is given by s_1 = \pi t \cos(\pi x t).
  • s_2 represents the source term for the flow rate and includes contributions from spatial and temporal variations as well as friction effects.

The radius R is computed using the radius function, and the friction coefficient k is obtained using the friction function.

This function is useful for evaluating the correctness of source term handling in numerical solvers.

source
diff --git a/dev/math/index.html b/dev/math/index.html index acaf356..3258b33 100644 --- a/dev/math/index.html +++ b/dev/math/index.html @@ -1,2 +1,2 @@ -Mathematics · BloodFlowTrixi.jl
GitHub

Modèles mathématiques 1D et 2D pour l’écoulement sanguin

Modèle 1D

Le modèle 1D repose sur une intégration par section des équations de Navier-Stokes sous l’hypothèse d’un écoulement incompressible dans des artères supposées fines. Ce modèle est particulièrement adapté pour des études globales du réseau artériel, où la géométrie est approximativement linéaire ou faiblement courbée.

Hypothèses et simplifications

  • L’écoulement est considéré comme incompressible.
  • L’artère est modélisée comme un tube cylindrique de section variable en fonction de la pression.
  • Un profil de vitesse parabolique est utilisé, permettant une moyennisation sur la section transversale de l’artère.

Équations principales

Les équations dérivées sont un système d’équations hyperboliques aux dérivées partielles décrivant la conservation de la masse et de la quantité de mouvement :

  1. Conservation de la masse :

\[∂_t A + ∂_x Q = 0\]

  1. Conservation de la quantité de mouvement :

\[∂_t Q + ∂_x \left( \alpha \frac{Q^2}{A} + \frac{1}{\rho} A P(A, x) \right) = \frac{1}{\rho} P(A, x) ∂_x A - K \frac{Q}{A}\]

Énergie et relation d’entropie du modèle 1D

L’énergie associée au système est donnée par :

\[E(t, x) = \frac{A u_x^2}{2} + \frac{1}{\rho} A P(A, x) - \frac{\beta(x)}{3 \rho A_0(x)} A^{3/2}\]

La relation d’entropie vérifiée par cette énergie est :

\[∂_t E + ∂_x \left( \left( E + \frac{\beta(x)}{3 \rho A_0(x)} A^{3/2} \right) u_x \right) = ∂_x \left( 3 \nu A ∂_x \left( \frac{Q}{A} \right) \right) u_x + \frac{2 \pi R k}{1 - R k / 4 \nu} u_x^2 ≤ 0\]

Sous des conditions aux limites nulles :

\[∂_t \left( \int_0^L E \, dx \right) = - 3 \nu \int_0^L A (∂_x u_x)^2 \, dx - \frac{2 \pi R k}{1 - R k / 4 \nu} \int_0^L u_x^2 \, dx < 0\]

Modèle 2D

Le modèle 2D est dérivé à partir d’une intégration radiale des équations de Navier-Stokes, permettant de mieux représenter les effets locaux dans des configurations géométriques complexes, comme les bifurcations artérielles et les anévrismes sévères.

Hypothèses et simplifications

  • L’écoulement est supposé incompressible.
  • La géométrie de l’artère est décrite à l’aide d’un système de coordonnées curvilignes (( s, heta )).
  • Le profil de vitesse est obtenu sans recourir à un ansatz spécifique.

Équations principales

  1. Conservation de la masse :

\[∂_t A + ∂_θ \left( \frac{Q_{Rθ}}{A} \right) + ∂_s(Q_s) = 0\]

  1. Conservation de la quantité de mouvement (composante radiale et axiale) :

\[∂_t (Q_{Rθ}) + ∂_θ \left( \frac{Q_{Rθ}^2}{2 A^2} + A P \right) + ∂_s \left( \frac{Q_{Rθ} Q_s}{A} \right) = \frac{2 R}{3} C \sin θ \frac{Q_s^2}{A} + \frac{2 R k Q_{Rθ}}{A} + ∂_θ (A P)\]

\[∂_t (Q_s) + ∂_θ \left( \frac{Q_s Q_{Rθ}}{A^2} \right) + ∂_s \left( \frac{Q_s^2}{A} - \frac{Q_{Rθ}^2}{2 A^2} + A P \right) = - \frac{2 R}{3} C \sin θ \frac{Q_{Rθ} Q_s}{A^2} + \frac{k R Q_s}{A} + ∂_s (A P)\]

Énergie et relation d’entropie du modèle 2D

L’énergie associée au système est donnée par :

\[E(t, θ, s) = A \left( \frac{9}{8} u_θ^2 + \frac{u_s^2}{2} + p \right) - p̃\]

La relation d’entropie correspondante est :

\[∂_t E + ∂_θ \left( \frac{3}{2} \frac{u_θ}{R} \left( E + p̃ - \frac{9}{16} A u_θ^2 \right) \right) + ∂_s \left( u_s \left( E + p̃ - \frac{9}{16} A u_θ^2 \right) \right) = \frac{9}{4} R k u_θ^2 + k R u_s^2 ≤ 0\]

Cette relation garantit que l’énergie décroît localement dans le temps, ce qui assure la stabilité du modèle.

Comparaison des modèles 1D et 2D

  • Modèle 1D :
    • Rapide et efficace pour des simulations globales sur de grands réseaux artériels.
    • Bien adapté pour des géométries simples ou faiblement courbées.
    • Coût de calcul très faible.
  • Modèle 2D :
    • Plus précis pour des géométries complexes (bifurcations, anévrismes).
    • Permet de mieux capturer les effets locaux et les interactions fluide-structure.
    • Coût de calcul modéré par rapport aux modèles tridimensionnels (NS-FSI 3D).

L’utilisation combinée de ces deux modèles permet une alternative efficace aux simulations 3D, tout en offrant un bon compromis entre précision et coût de calcul.

+Mathematics · BloodFlowTrixi.jl
GitHub

Modèles mathématiques 1D et 2D pour l’écoulement sanguin

Modèle 1D

Le modèle 1D repose sur une intégration par section des équations de Navier-Stokes sous l’hypothèse d’un écoulement incompressible dans des artères supposées fines. Ce modèle est particulièrement adapté pour des études globales du réseau artériel, où la géométrie est approximativement linéaire ou faiblement courbée.

Hypothèses et simplifications

  • L’écoulement est considéré comme incompressible.
  • L’artère est modélisée comme un tube cylindrique de section variable en fonction de la pression.
  • Un profil de vitesse parabolique est utilisé, permettant une moyennisation sur la section transversale de l’artère.

Équations principales

Les équations dérivées sont un système d’équations hyperboliques aux dérivées partielles décrivant la conservation de la masse et de la quantité de mouvement :

  1. Conservation de la masse :

\[∂_t A + ∂_x Q = 0\]

  1. Conservation de la quantité de mouvement :

\[∂_t Q + ∂_x \left( \alpha \frac{Q^2}{A} + \frac{1}{\rho} A P(A, x) \right) - \partial_x \left( A \partial_x\left(\frac Q A\right) \right) = \frac{1}{\rho} P(A, x) ∂_x A - K \frac{Q}{A}\]

Énergie et relation d’entropie du modèle 1D

L’énergie associée au système est donnée par :

\[E(t, x) = \frac{A u_x^2}{2} + \frac{1}{\rho} A P(A, x) - \frac{\beta(x)}{3 \rho A_0(x)} A^{3/2}\]

La relation d’entropie vérifiée par cette énergie est :

\[∂_t E + ∂_x \left( \left( E + \frac{\beta(x)}{3 \rho A_0(x)} A^{3/2} \right) u_x \right) = ∂_x \left( 3 \nu A ∂_x \left( \frac{Q}{A} \right) \right) u_x + \frac{2 \pi R k}{1 - R k / 4 \nu} u_x^2 ≤ 0\]

Sous des conditions aux limites nulles :

\[∂_t \left( \int_0^L E \, dx \right) = - 3 \nu \int_0^L A (∂_x u_x)^2 \, dx - \frac{2 \pi R k}{1 - R k / 4 \nu} \int_0^L u_x^2 \, dx < 0\]

Modèle 2D

Le modèle 2D est dérivé à partir d’une intégration radiale des équations de Navier-Stokes, permettant de mieux représenter les effets locaux dans des configurations géométriques complexes, comme les bifurcations artérielles et les anévrismes sévères.

Hypothèses et simplifications

  • L’écoulement est supposé incompressible.
  • La géométrie de l’artère est décrite à l’aide d’un système de coordonnées curvilignes (( s, heta )).
  • Le profil de vitesse est obtenu sans recourir à un ansatz spécifique.

Équations principales

  1. Conservation de la masse :

\[∂_t A + ∂_θ \left( \frac{Q_{Rθ}}{A} \right) + ∂_s(Q_s) = 0\]

  1. Conservation de la quantité de mouvement (composante radiale et axiale) :

\[∂_t (Q_{Rθ}) + ∂_θ \left( \frac{Q_{Rθ}^2}{2 A^2} + A P \right) + ∂_s \left( \frac{Q_{Rθ} Q_s}{A} \right) = \frac{2 R}{3} C \sin θ \frac{Q_s^2}{A} + \frac{2 R k Q_{Rθ}}{A} + ∂_θ (A P)\]

\[∂_t (Q_s) + ∂_θ \left( \frac{Q_s Q_{Rθ}}{A^2} \right) + ∂_s \left( \frac{Q_s^2}{A} - \frac{Q_{Rθ}^2}{2 A^2} + A P \right) = - \frac{2 R}{3} C \sin θ \frac{Q_{Rθ} Q_s}{A^2} + \frac{k R Q_s}{A} + ∂_s (A P)\]

Énergie et relation d’entropie du modèle 2D

L’énergie associée au système est donnée par :

\[E(t, θ, s) = A \left( \frac{9}{8} u_θ^2 + \frac{u_s^2}{2} + p \right) - p̃\]

La relation d’entropie correspondante est :

\[∂_t E + ∂_θ \left( \frac{3}{2} \frac{u_θ}{R} \left( E + p̃ - \frac{9}{16} A u_θ^2 \right) \right) + ∂_s \left( u_s \left( E + p̃ - \frac{9}{16} A u_θ^2 \right) \right) = \frac{9}{4} R k u_θ^2 + k R u_s^2 ≤ 0\]

Cette relation garantit que l’énergie décroît localement dans le temps, ce qui assure la stabilité du modèle.

Comparaison des modèles 1D et 2D

  • Modèle 1D :
    • Rapide et efficace pour des simulations globales sur de grands réseaux artériels.
    • Bien adapté pour des géométries simples ou faiblement courbées.
    • Coût de calcul très faible.
  • Modèle 2D :
    • Plus précis pour des géométries complexes (bifurcations, anévrismes).
    • Permet de mieux capturer les effets locaux et les interactions fluide-structure.
    • Coût de calcul modéré par rapport aux modèles tridimensionnels (NS-FSI 3D).

L’utilisation combinée de ces deux modèles permet une alternative efficace aux simulations 3D, tout en offrant un bon compromis entre précision et coût de calcul.

diff --git a/dev/search_index.js b/dev/search_index.js index 8008bcc..552b845 100644 --- a/dev/search_index.js +++ b/dev/search_index.js @@ -1,3 +1,3 @@ var documenterSearchIndex = {"docs": -[{"location":"tuto/#Tutorial-for-the-1D-model","page":"Tutorial","title":"Tutorial for the 1D model","text":"","category":"section"},{"location":"tuto/","page":"Tutorial","title":"Tutorial","text":"In this section, we describe how to use BloodFlowTrixi.jl with Trixi.jl. This tutorial will guide you through setting up and running a 1D blood flow simulation, including mesh creation, boundary conditions, numerical fluxes, and visualization of results.","category":"page"},{"location":"tuto/#Packages","page":"Tutorial","title":"Packages","text":"","category":"section"},{"location":"tuto/","page":"Tutorial","title":"Tutorial","text":"Before starting, ensure that the required packages are loaded:","category":"page"},{"location":"tuto/","page":"Tutorial","title":"Tutorial","text":"using Trixi\nusing BloodFlowTrixi\nusing OrdinaryDiffEq\nusing Plots","category":"page"},{"location":"tuto/","page":"Tutorial","title":"Tutorial","text":"First, we need to choose the equation that describes the blood flow dynamics:","category":"page"},{"location":"tuto/","page":"Tutorial","title":"Tutorial","text":"eq = BloodFlowEquations1D(; h=0.1)","category":"page"},{"location":"tuto/","page":"Tutorial","title":"Tutorial","text":"Here, h represents a parameter related to the initial condition or model scaling.","category":"page"},{"location":"tuto/#Mesh-and-boundary-conditions","page":"Tutorial","title":"Mesh and boundary conditions","text":"","category":"section"},{"location":"tuto/","page":"Tutorial","title":"Tutorial","text":"We begin by defining a one-dimensional Tree mesh, which discretizes the spatial domain:","category":"page"},{"location":"tuto/","page":"Tutorial","title":"Tutorial","text":"mesh = TreeMesh(0.0, 40.0, initial_refinement_level=6, n_cells_max=10^4, periodicity=false)","category":"page"},{"location":"tuto/","page":"Tutorial","title":"Tutorial","text":"This generates a non-periodic mesh for the interval 0 40, with 2^initial-refinement-level+1-1 cells. The parameter initial_refinement_level controls the initial number of cells, while n_cells_max specifies the maximum number of cells allowed during mesh refinement.","category":"page"},{"location":"tuto/","page":"Tutorial","title":"Tutorial","text":"In Trixi.jl, the Tree mesh has two labeled boundaries: ***xneg*** (left boundary) and ***xpos*** (right boundary). These labels are used to apply boundary conditions:","category":"page"},{"location":"tuto/","page":"Tutorial","title":"Tutorial","text":"bc = (\n x_neg = boundary_condition_pressure_in,\n x_pos = Trixi.BoundaryConditionDoNothing()\n)","category":"page"},{"location":"tuto/","page":"Tutorial","title":"Tutorial","text":"boundary_condition_pressure_in applies a pressure inflow condition at the left boundary.\nTrixi.BoundaryConditionDoNothing() specifies a \"do nothing\" boundary condition at the right boundary, meaning no flux is imposed.","category":"page"},{"location":"tuto/#Boundary-condition-implementation","page":"Tutorial","title":"Boundary condition implementation","text":"","category":"section"},{"location":"tuto/","page":"Tutorial","title":"Tutorial","text":"The inflow boundary condition is defined as:","category":"page"},{"location":"tuto/","page":"Tutorial","title":"Tutorial","text":"boundary_condition_pressure_in(u_inner, orientation_or_normal, direction, x, t, surface_flux_function, eq::BloodFlowEquations1D)","category":"page"},{"location":"tuto/","page":"Tutorial","title":"Tutorial","text":"This function applies a time-dependent pressure inflow condition.","category":"page"},{"location":"tuto/#Parameters","page":"Tutorial","title":"Parameters","text":"","category":"section"},{"location":"tuto/","page":"Tutorial","title":"Tutorial","text":"u_inner: State vector inside the domain near the boundary.\norientation_or_normal: Normal orientation of the boundary.\ndirection: Integer indicating the boundary direction.\nx: Position vector.\nt: Time scalar.\nsurface_flux_function: Function to compute flux at the boundary.\neq: Instance of BloodFlowEquations1D.","category":"page"},{"location":"tuto/#Returns","page":"Tutorial","title":"Returns","text":"","category":"section"},{"location":"tuto/","page":"Tutorial","title":"Tutorial","text":"The boundary flux is computed based on the inflow pressure: $ P{\\text{in}} = \\begin{cases} 2 \\times 10^4 \\sin^2\\left(\\frac{\\pi t}{0.125}\\right) & \\text{if } t < 0.125 \\\n0 & \\text{otherwise} \\end{cases} $ This time-dependent inflow pressure mimics a pulsatile flow, typical in arterial blood flow. The inflow area A{\\text{in}}$ is determined using the inverse pressure relation, ensuring consistency with the physical model.","category":"page"},{"location":"tuto/#Numerical-flux","page":"Tutorial","title":"Numerical flux","text":"","category":"section"},{"location":"tuto/","page":"Tutorial","title":"Tutorial","text":"To compute fluxes at cell interfaces, we use a combination of conservative and non-conservative fluxes:","category":"page"},{"location":"tuto/","page":"Tutorial","title":"Tutorial","text":"volume_flux = (flux_lax_friedrichs, flux_nonconservative)\nsurface_flux = (flux_lax_friedrichs, flux_nonconservative)","category":"page"},{"location":"tuto/","page":"Tutorial","title":"Tutorial","text":"flux_lax_friedrichs is a standard numerical flux for hyperbolic conservation laws.\nflux_nonconservative handles the non-conservative terms in the model, particularly those related to pressure discontinuities.","category":"page"},{"location":"tuto/","page":"Tutorial","title":"Tutorial","text":"The non-conservative flux function is defined as:","category":"page"},{"location":"tuto/","page":"Tutorial","title":"Tutorial","text":"flux_nonconservative(u_ll, u_rr, orientation::Integer, eq::BloodFlowEquations1D)","category":"page"},{"location":"tuto/#Parameters-2","page":"Tutorial","title":"Parameters","text":"","category":"section"},{"location":"tuto/","page":"Tutorial","title":"Tutorial","text":"u_ll: Left state vector.\nu_rr: Right state vector.\norientation::Integer: Orientation index.\neq: Instance of BloodFlowEquations1D.","category":"page"},{"location":"tuto/#Returns-2","page":"Tutorial","title":"Returns","text":"","category":"section"},{"location":"tuto/","page":"Tutorial","title":"Tutorial","text":"The function returns the non-conservative flux vector, which is essential for capturing sharp pressure changes in the simulation.","category":"page"},{"location":"tuto/#Basis-functions-and-Shock-Capturing-DG-scheme","page":"Tutorial","title":"Basis functions and Shock Capturing DG scheme","text":"","category":"section"},{"location":"tuto/","page":"Tutorial","title":"Tutorial","text":"To approximate the solution, we use polynomial basis functions:","category":"page"},{"location":"tuto/","page":"Tutorial","title":"Tutorial","text":"basis = LobattoLegendreBasis(2)","category":"page"},{"location":"tuto/","page":"Tutorial","title":"Tutorial","text":"This defines a Lobatto-Legendre basis of polynomial degree 2, which is commonly used in high-order methods like Discontinuous Galerkin (DG) schemes.","category":"page"},{"location":"tuto/","page":"Tutorial","title":"Tutorial","text":"We then define an indicator for shock capturing, focusing on the first variable (area perturbation $a$):","category":"page"},{"location":"tuto/","page":"Tutorial","title":"Tutorial","text":"id = IndicatorHennemannGassner(eq, basis; variable=first)","category":"page"},{"location":"tuto/","page":"Tutorial","title":"Tutorial","text":"This indicator helps detect shocks or discontinuities in the solution and applies appropriate stabilization.","category":"page"},{"location":"tuto/","page":"Tutorial","title":"Tutorial","text":"The solver is defined as:","category":"page"},{"location":"tuto/","page":"Tutorial","title":"Tutorial","text":"vol = VolumeIntegralShockCapturingHG(id, volume_flux_dg=volume_flux, volume_flux_fv=surface_flux)\nsolver = DGSEM(basis, surface_flux, vol)","category":"page"},{"location":"tuto/","page":"Tutorial","title":"Tutorial","text":"Here, DGSEM represents the Discontinuous Galerkin Spectral Element Method, a high-order accurate scheme suitable for hyperbolic problems.","category":"page"},{"location":"tuto/#Semi-discretization","page":"Tutorial","title":"Semi-discretization","text":"","category":"section"},{"location":"tuto/","page":"Tutorial","title":"Tutorial","text":"We are now ready to semi-discretize the problem:","category":"page"},{"location":"tuto/","page":"Tutorial","title":"Tutorial","text":"semi = SemidiscretizationHyperbolic(\n mesh,\n eq,\n initial_condition_simple,\n source_terms = source_term_simple,\n solver,\n boundary_conditions=bc\n)","category":"page"},{"location":"tuto/","page":"Tutorial","title":"Tutorial","text":"This step sets up the semi-discretized form of the PDE, which will be advanced in time using an ODE solver.","category":"page"},{"location":"tuto/#Source-term","page":"Tutorial","title":"Source term","text":"","category":"section"},{"location":"tuto/","page":"Tutorial","title":"Tutorial","text":"The source term accounts for additional forces acting on the blood flow, such as friction:","category":"page"},{"location":"tuto/","page":"Tutorial","title":"Tutorial","text":"source_term_simple(u, x, t, eq::BloodFlowEquations1D)","category":"page"},{"location":"tuto/#Parameters-3","page":"Tutorial","title":"Parameters","text":"","category":"section"},{"location":"tuto/","page":"Tutorial","title":"Tutorial","text":"u: State vector containing area perturbation, flow rate, elasticity modulus, and reference area.\nx: Position vector.\nt: Time scalar.\neq::BloodFlowEquations1D: Instance of the blood flow model.","category":"page"},{"location":"tuto/#Returns-3","page":"Tutorial","title":"Returns","text":"","category":"section"},{"location":"tuto/","page":"Tutorial","title":"Tutorial","text":"The source term vector is given by:","category":"page"},{"location":"tuto/","page":"Tutorial","title":"Tutorial","text":"s_1 = 0\n(no source for area perturbation).\ns_2 = frac2 pi k QR A\n, representing frictional effects.","category":"page"},{"location":"tuto/","page":"Tutorial","title":"Tutorial","text":"The friction coefficient k is computed using a model-specific friction function, and the radius R is obtained from the state vector using the radius function.","category":"page"},{"location":"tuto/#Initial-condition","page":"Tutorial","title":"Initial condition","text":"","category":"section"},{"location":"tuto/","page":"Tutorial","title":"Tutorial","text":"The initial condition specifies the starting state of the simulation:","category":"page"},{"location":"tuto/","page":"Tutorial","title":"Tutorial","text":"initial_condition_simple(x, t, eq::BloodFlowEquations1D; R0=2.0)","category":"page"},{"location":"tuto/","page":"Tutorial","title":"Tutorial","text":"This function generates a simple initial condition with a uniform radius R0.","category":"page"},{"location":"tuto/#Parameters-4","page":"Tutorial","title":"Parameters","text":"","category":"section"},{"location":"tuto/","page":"Tutorial","title":"Tutorial","text":"x: Position vector.\nt: Time scalar.\neq::BloodFlowEquations1D: Instance of the blood flow model.\nR0: Initial radius (default: 2.0).","category":"page"},{"location":"tuto/#Returns-4","page":"Tutorial","title":"Returns","text":"","category":"section"},{"location":"tuto/","page":"Tutorial","title":"Tutorial","text":"The function returns a state vector with:","category":"page"},{"location":"tuto/","page":"Tutorial","title":"Tutorial","text":"Zero initial area perturbation.\nZero initial flow rate.\nConstant elasticity modulus.\nReference area A_0 = pi R_0^2.","category":"page"},{"location":"tuto/","page":"Tutorial","title":"Tutorial","text":"This simple initial condition is suitable for testing the model without introducing complex dynamics.","category":"page"},{"location":"tuto/#Run-the-simulation","page":"Tutorial","title":"Run the simulation","text":"","category":"section"},{"location":"tuto/","page":"Tutorial","title":"Tutorial","text":"First, we discretize the problem in time:","category":"page"},{"location":"tuto/","page":"Tutorial","title":"Tutorial","text":"Trixi.default_analysis_integrals(::BloodFlowEquations1D) = ()\ntspan = (0.0, 0.5)\node = semidiscretize(semi, tspan)","category":"page"},{"location":"tuto/","page":"Tutorial","title":"Tutorial","text":"Here, tspan defines the time interval for the simulation.","category":"page"},{"location":"tuto/","page":"Tutorial","title":"Tutorial","text":"Next, we add some callbacks to monitor the simulation:","category":"page"},{"location":"tuto/","page":"Tutorial","title":"Tutorial","text":"summary_callback = SummaryCallback()\nanalysis_callback = AnalysisCallback(semi, interval=200)\nstepsize_callback = StepsizeCallback(; cfl=0.5)\ncallbacks = CallbackSet(summary_callback, analysis_callback, stepsize_callback)","category":"page"},{"location":"tuto/","page":"Tutorial","title":"Tutorial","text":"SummaryCallback provides a summary of the simulation progress.\nAnalysisCallback computes analysis metrics at specified intervals.\nStepsizeCallback adjusts the time step based on the CFL condition.","category":"page"},{"location":"tuto/","page":"Tutorial","title":"Tutorial","text":"Finally, we solve the problem:","category":"page"},{"location":"tuto/","page":"Tutorial","title":"Tutorial","text":"dt = stepsize_callback(ode)\nsol = solve(ode, SSPRK33(), dt=dt, dtmax=1e-4, dtmin=1e-11,\n save_everystep=false, saveat=0.002, callback=callbacks)","category":"page"},{"location":"tuto/","page":"Tutorial","title":"Tutorial","text":"Here, SSPRK33() is a third-order Strong Stability Preserving Runge-Kutta method, suitable for hyperbolic PDEs.","category":"page"},{"location":"tuto/#Plot-the-results","page":"Tutorial","title":"Plot the results","text":"","category":"section"},{"location":"tuto/","page":"Tutorial","title":"Tutorial","text":"The results can be visualized using the following code:","category":"page"},{"location":"tuto/","page":"Tutorial","title":"Tutorial","text":"@gif for i in eachindex(sol)\n a1 = sol[i][1:4:end]\n Q1 = sol[i][2:4:end]\n A01 = sol[i][4:4:end]\n A1 = A01 .+ a1\n plot(Q1 ./ A1, lw=4, color=:red, ylim=(-10, 50), label=\"velocity\", legend=:bottomleft)\nend","category":"page"},{"location":"tuto/","page":"Tutorial","title":"Tutorial","text":"This code generates an animated GIF showing the evolution of the velocity profile over time. The velocity is computed as QA, where Q is the flow rate, and A is the cross-sectional area.","category":"page"},{"location":"tuto/#Plain-code","page":"Tutorial","title":"Plain code","text":"","category":"section"},{"location":"tuto/","page":"Tutorial","title":"Tutorial","text":"using Trixi\nusing BloodFlowTrixi\nusing OrdinaryDiffEq,Plots\neq = BloodFlowEquations1D(;h=0.1)\nmesh = TreeMesh(0.0,40.0,initial_refinement_level=6,n_cells_max=10^4,periodicity=false)\nbc = (\n x_neg = boundary_condition_pressure_in,\n x_pos = Trixi.BoundaryConditionDoNothing()\n )\nvolume_flux = (flux_lax_friedrichs,flux_nonconservative)\nsurface_flux = (flux_lax_friedrichs,flux_nonconservative)\nbasis = LobattoLegendreBasis(2)\nid = IndicatorHennemannGassner(eq,basis;variable=first)\nvol = VolumeIntegralShockCapturingHG(id,volume_flux_dg=volume_flux,volume_flux_fv=surface_flux)\nsolver = DGSEM(basis,surface_flux,vol)\nsemi = SemidiscretizationHyperbolic(mesh,\neq,\ninitial_condition_simple,\nsource_terms = source_term_simple,\nsolver,\nboundary_conditions=bc)\nTrixi.default_analysis_integrals(::BloodFlowEquations1D) = ()\ntspan = (0.0, 0.5)\node = semidiscretize(semi, tspan)\nsummary_callback = SummaryCallback()\nanalysis_callback = AnalysisCallback(semi, interval = 200)\nstepsize_callback = StepsizeCallback(; cfl=0.5)\ncallbacks = CallbackSet(summary_callback,analysis_callback,stepsize_callback)\ndt = stepsize_callback(ode)\nsol = solve(ode, SSPRK33(), dt = dt, dtmax = 1e-4,dtmin = 1e-11,\n save_everystep = false,saveat = 0.002, callback = callbacks)\n\n@gif for i in eachindex(sol)\n a1 = sol[i][1:4:end]\n Q1 = sol[i][2:4:end]\n A01 = sol[i][4:4:end]\n A1 = A01.+a1\n plot(Q1./A1,lw=4,color=:red,ylim=(-10,50),label=\"velocity\",legend=:bottomleft)\nend","category":"page"},{"location":"tuto/","page":"Tutorial","title":"Tutorial","text":"(Image: Alt Text)","category":"page"},{"location":"tuto/#Tutorial-for-the-2D-model","page":"Tutorial","title":"Tutorial for the 2D model","text":"","category":"section"},{"location":"tuto/","page":"Tutorial","title":"Tutorial","text":"In this section, we describe how to use BloodFlowTrixi.jl with Trixi.jl. This tutorial will guide you through setting up and running a 2D blood flow simulation, including mesh creation, boundary conditions, numerical fluxes, and visualization of results.","category":"page"},{"location":"tuto/#Packages-2","page":"Tutorial","title":"Packages","text":"","category":"section"},{"location":"tuto/","page":"Tutorial","title":"Tutorial","text":"Before starting, ensure that the required packages are loaded:","category":"page"},{"location":"tuto/","page":"Tutorial","title":"Tutorial","text":"using Trixi\nusing BloodFlowTrixi\nusing OrdinaryDiffEq\nusing Plots","category":"page"},{"location":"tuto/","page":"Tutorial","title":"Tutorial","text":"First, we need to choose the equation that describes the blood flow dynamics:","category":"page"},{"location":"tuto/","page":"Tutorial","title":"Tutorial","text":"eq = BloodFlowEquations2D(; h=0.1)","category":"page"},{"location":"tuto/","page":"Tutorial","title":"Tutorial","text":"Here, h represents a parameter related to the initial condition or model scaling.","category":"page"},{"location":"tuto/#Mesh-and-boundary-conditions-2","page":"Tutorial","title":"Mesh and boundary conditions","text":"","category":"section"},{"location":"tuto/","page":"Tutorial","title":"Tutorial","text":"We begin by defining a two-dimensional P4est mesh, which discretizes the spatial domain:","category":"page"},{"location":"tuto/","page":"Tutorial","title":"Tutorial","text":"mesh = P4estMesh(\n (2,4),\n polydeg= 2,\n coordinates_min =(0.0,0.0),\n coordinates_max = (2*pi,40.0),\n initial_refinement_level = 4,\n periodicity = (true, false)\n)","category":"page"},{"location":"tuto/","page":"Tutorial","title":"Tutorial","text":"This generates a non-periodic mesh for the domain 02pi times 0 40, with 2times 4times 4^textinitial-refinement-level cells. ","category":"page"},{"location":"tuto/","page":"Tutorial","title":"Tutorial","text":"In Trixi.jl, the P4est mesh has four labeled boundaries: ***xneg*** (left boundary), ***xpos*** (right boundary), ***yneg*** (bottom boundary), and ***ypos*** (top boundary). These labels are used to apply boundary conditions:","category":"page"},{"location":"tuto/","page":"Tutorial","title":"Tutorial","text":"bc = Dict(\n :y_neg =>boundary_condition_pressure_in,\n :y_pos => Trixi.BoundaryConditionDoNothing()\n )\n","category":"page"},{"location":"tuto/","page":"Tutorial","title":"Tutorial","text":"boundary_condition_pressure_in applies a pressure inflow condition at the bottom boundary.\nTrixi.BoundaryConditionDoNothing() specifies a \"do nothing\" boundary condition at the right boundary, meaning no flux is imposed.","category":"page"},{"location":"tuto/#Boundary-condition-implementation-2","page":"Tutorial","title":"Boundary condition implementation","text":"","category":"section"},{"location":"tuto/","page":"Tutorial","title":"Tutorial","text":"The inflow boundary condition is defined as:","category":"page"},{"location":"tuto/","page":"Tutorial","title":"Tutorial","text":"boundary_condition_pressure_in(u_inner, normal, x, t, surface_flux_function, eq::BloodFlowEquations2D)","category":"page"},{"location":"tuto/","page":"Tutorial","title":"Tutorial","text":"This function applies a time-dependent pressure inflow condition.","category":"page"},{"location":"tuto/#Parameters-5","page":"Tutorial","title":"Parameters","text":"","category":"section"},{"location":"tuto/","page":"Tutorial","title":"Tutorial","text":"u_inner: State vector inside the domain near the boundary.\nnormal: Normal of the boundary.\nx: Position vector.\nt: Time scalar.\nsurface_flux_function: Function to compute flux at the boundary.\neq: Instance of BloodFlowEquations2D.","category":"page"},{"location":"tuto/#Returns-5","page":"Tutorial","title":"Returns","text":"","category":"section"},{"location":"tuto/","page":"Tutorial","title":"Tutorial","text":"The boundary flux is computed based on the inflow pressure: $ P{\\text{in}} = \\begin{cases} 2 \\times 10^4 \\sin^2\\left(\\frac{\\pi t}{0.125}\\right) & \\text{if } t < 0.125 \\\n0 & \\text{otherwise} \\end{cases} $ This time-dependent inflow pressure mimics a pulsatile flow, typical in arterial blood flow. The inflow area A{\\text{in}}$ is determined using the inverse pressure relation, ensuring consistency with the physical model.","category":"page"},{"location":"tuto/#Numerical-flux-2","page":"Tutorial","title":"Numerical flux","text":"","category":"section"},{"location":"tuto/","page":"Tutorial","title":"Tutorial","text":"To compute fluxes at cell interfaces, we use a combination of conservative and non-conservative fluxes:","category":"page"},{"location":"tuto/","page":"Tutorial","title":"Tutorial","text":"volume_flux = (flux_lax_friedrichs, flux_nonconservative)\nsurface_flux = (flux_lax_friedrichs, flux_nonconservative)","category":"page"},{"location":"tuto/","page":"Tutorial","title":"Tutorial","text":"flux_lax_friedrichs is a standard numerical flux for hyperbolic conservation laws.\nflux_nonconservative handles the non-conservative terms in the model, particularly those related to pressure discontinuities.","category":"page"},{"location":"tuto/","page":"Tutorial","title":"Tutorial","text":"The non-conservative flux function is defined as:","category":"page"},{"location":"tuto/","page":"Tutorial","title":"Tutorial","text":"flux_nonconservative(u_ll, u_rr, normal::Integer, eq::BloodFlowEquations2D)","category":"page"},{"location":"tuto/#Parameters-6","page":"Tutorial","title":"Parameters","text":"","category":"section"},{"location":"tuto/","page":"Tutorial","title":"Tutorial","text":"u_ll: Left state vector.\nu_rr: Right state vector.\nnormal: normal vector.\neq: Instance of BloodFlowEquations2D.","category":"page"},{"location":"tuto/#Returns-6","page":"Tutorial","title":"Returns","text":"","category":"section"},{"location":"tuto/","page":"Tutorial","title":"Tutorial","text":"The function returns the non-conservative flux vector, which is essential for capturing sharp pressure changes in the simulation.","category":"page"},{"location":"tuto/#Basis-functions-and-Shock-Capturing-DG-scheme-2","page":"Tutorial","title":"Basis functions and Shock Capturing DG scheme","text":"","category":"section"},{"location":"tuto/","page":"Tutorial","title":"Tutorial","text":"To approximate the solution, we use polynomial basis functions:","category":"page"},{"location":"tuto/","page":"Tutorial","title":"Tutorial","text":"basis = LobattoLegendreBasis(2)","category":"page"},{"location":"tuto/","page":"Tutorial","title":"Tutorial","text":"This defines a Lobatto-Legendre basis of polynomial degree 2, which is commonly used in high-order methods like Discontinuous Galerkin (DG) schemes.","category":"page"},{"location":"tuto/","page":"Tutorial","title":"Tutorial","text":"We then define an indicator for shock capturing, focusing on the first variable (area perturbation $a$):","category":"page"},{"location":"tuto/","page":"Tutorial","title":"Tutorial","text":"id = IndicatorHennemannGassner(eq, basis; variable=first)","category":"page"},{"location":"tuto/","page":"Tutorial","title":"Tutorial","text":"This indicator helps detect shocks or discontinuities in the solution and applies appropriate stabilization.","category":"page"},{"location":"tuto/","page":"Tutorial","title":"Tutorial","text":"The solver is defined as:","category":"page"},{"location":"tuto/","page":"Tutorial","title":"Tutorial","text":"vol = VolumeIntegralShockCapturingHG(id, volume_flux_dg=volume_flux, volume_flux_fv=surface_flux)\nsolver = DGSEM(basis, surface_flux, vol)","category":"page"},{"location":"tuto/","page":"Tutorial","title":"Tutorial","text":"Here, DGSEM represents the Discontinuous Galerkin Spectral Element Method, a high-order accurate scheme suitable for hyperbolic problems.","category":"page"},{"location":"tuto/#Semi-discretization-2","page":"Tutorial","title":"Semi-discretization","text":"","category":"section"},{"location":"tuto/","page":"Tutorial","title":"Tutorial","text":"We are now ready to semi-discretize the problem:","category":"page"},{"location":"tuto/","page":"Tutorial","title":"Tutorial","text":"semi = SemidiscretizationHyperbolic(\n mesh,\n eq,\n initial_condition_simple,\n source_terms = source_term_simple,\n solver,\n boundary_conditions=bc\n)","category":"page"},{"location":"tuto/","page":"Tutorial","title":"Tutorial","text":"This step sets up the semi-discretized form of the PDE, which will be advanced in time using an ODE solver.","category":"page"},{"location":"tuto/#Source-term-2","page":"Tutorial","title":"Source term","text":"","category":"section"},{"location":"tuto/","page":"Tutorial","title":"Tutorial","text":"The source term accounts for additional forces acting on the blood flow, such as friction:","category":"page"},{"location":"tuto/","page":"Tutorial","title":"Tutorial","text":"source_term_simple(u, x, t, eq::BloodFlowEquations2D)","category":"page"},{"location":"tuto/#Parameters-7","page":"Tutorial","title":"Parameters","text":"","category":"section"},{"location":"tuto/","page":"Tutorial","title":"Tutorial","text":"u: State vector containing area perturbation, flow rate, elasticity modulus, and reference area.\nx: Position vector.\nt: Time scalar.\neq::BloodFlowEquations2D: Instance of the blood flow model.","category":"page"},{"location":"tuto/#Returns-7","page":"Tutorial","title":"Returns","text":"","category":"section"},{"location":"tuto/","page":"Tutorial","title":"Tutorial","text":"The source term vector is given by:","category":"page"},{"location":"tuto/","page":"Tutorial","title":"Tutorial","text":"s_1 = 0\n(no source for area perturbation).\ns_2 = frac2R3mathcalC sin theta fracQ_s^2A + frac3Rk Q_RθA\n.\ns_3 = -frac2R3mathcalC sin theta fracQ_s Q_RthetaA + fracRkQ_sA\n.","category":"page"},{"location":"tuto/","page":"Tutorial","title":"Tutorial","text":"The friction coefficient k is computed using a model-specific friction function, and the radius R is obtained from the state vector using the radius function. Also, the curvature mathcalC is computed using the curvature function and is equal to 1 here.","category":"page"},{"location":"tuto/#Initial-condition-2","page":"Tutorial","title":"Initial condition","text":"","category":"section"},{"location":"tuto/","page":"Tutorial","title":"Tutorial","text":"The initial condition specifies the starting state of the simulation:","category":"page"},{"location":"tuto/","page":"Tutorial","title":"Tutorial","text":"initial_condition_simple(x, t, eq::BloodFlowEquations2D; R0=2.0)","category":"page"},{"location":"tuto/","page":"Tutorial","title":"Tutorial","text":"This function generates a simple initial condition with a uniform radius R0.","category":"page"},{"location":"tuto/#Parameters-8","page":"Tutorial","title":"Parameters","text":"","category":"section"},{"location":"tuto/","page":"Tutorial","title":"Tutorial","text":"x: Position vector.\nt: Time scalar.\neq::BloodFlowEquations2D: Instance of the blood flow model.\nR0: Initial radius (default: 2.0).","category":"page"},{"location":"tuto/#Returns-8","page":"Tutorial","title":"Returns","text":"","category":"section"},{"location":"tuto/","page":"Tutorial","title":"Tutorial","text":"The function returns a state vector with:","category":"page"},{"location":"tuto/","page":"Tutorial","title":"Tutorial","text":"Zero initial area perturbation.\nZero initial flow rate (in theta and s directions).\nConstant elasticity modulus.\nReference area A_0 = fracR_0^22.","category":"page"},{"location":"tuto/","page":"Tutorial","title":"Tutorial","text":"This simple initial condition is suitable for testing the model without introducing complex dynamics.","category":"page"},{"location":"tuto/#Run-the-simulation-2","page":"Tutorial","title":"Run the simulation","text":"","category":"section"},{"location":"tuto/","page":"Tutorial","title":"Tutorial","text":"First, we discretize the problem in time:","category":"page"},{"location":"tuto/","page":"Tutorial","title":"Tutorial","text":"Trixi.default_analysis_integrals(::BloodFlowEquations2D) = ()\ntspan = (0.0, 0.3)\node = semidiscretize(semi, tspan)","category":"page"},{"location":"tuto/","page":"Tutorial","title":"Tutorial","text":"Here, tspan defines the time interval for the simulation.","category":"page"},{"location":"tuto/","page":"Tutorial","title":"Tutorial","text":"Next, we add some callbacks to monitor the simulation:","category":"page"},{"location":"tuto/","page":"Tutorial","title":"Tutorial","text":"summary_callback = SummaryCallback()\nanalysis_callback = AnalysisCallback(semi, interval=200)\nstepsize_callback = StepsizeCallback(; cfl=0.5)\ncallbacks = CallbackSet(summary_callback, analysis_callback, stepsize_callback)","category":"page"},{"location":"tuto/","page":"Tutorial","title":"Tutorial","text":"SummaryCallback provides a summary of the simulation progress.\nAnalysisCallback computes analysis metrics at specified intervals.\nStepsizeCallback adjusts the time step based on the CFL condition.","category":"page"},{"location":"tuto/","page":"Tutorial","title":"Tutorial","text":"Finally, we solve the problem:","category":"page"},{"location":"tuto/","page":"Tutorial","title":"Tutorial","text":"dt = stepsize_callback(ode)\nsol = solve(ode, SSPRK33(), dt=dt, dtmax=1e-4, dtmin=1e-11,\n save_everystep=false, saveat=0.003, callback=callbacks)","category":"page"},{"location":"tuto/","page":"Tutorial","title":"Tutorial","text":"Here, SSPRK33() is a third-order Strong Stability Preserving Runge-Kutta method, suitable for hyperbolic PDEs.","category":"page"},{"location":"tuto/#Plot-the-results-2","page":"Tutorial","title":"Plot the results","text":"","category":"section"},{"location":"tuto/","page":"Tutorial","title":"Tutorial","text":"The results can be visualized using the following code:","category":"page"},{"location":"tuto/","page":"Tutorial","title":"Tutorial","text":"@gif for i in eachindex(sol)\npd = PlotData2D(sol[i],semi,solution_variables=cons2prim)\nplt1 = Plots.plot(pd[\"A\"],aspect_ratio=0.2)\nplt2 = Plots.plot(pd[\"wtheta\"],aspect_ratio=0.2)\nplt3 = Plots.plot(pd[\"ws\"],aspect_ratio=0.2)\nplt4 = Plots.plot(pd[\"P\"],aspect_ratio=0.2)\nplot(plt1,plt2,plt3,plt4,layout=(2,2))\nend","category":"page"},{"location":"tuto/","page":"Tutorial","title":"Tutorial","text":"This code generates an animated GIF showing the evolution of the velocity profile over time. The velocity is computed as QA, where Q is the flow rate, and A is the cross-sectional area.","category":"page"},{"location":"tuto/#Plain-code-2","page":"Tutorial","title":"Plain code","text":"","category":"section"},{"location":"tuto/","page":"Tutorial","title":"Tutorial","text":"using Trixi\nusing BloodFlowTrixi\nusing OrdinaryDiffEq,Plots\neq = BloodFlowEquations2D(;h=0.1)\nmesh = P4estMesh(\n (2,4),\n polydeg= 2,\n coordinates_min =(0.0,0.0),\n coordinates_max = (2*pi,40.0),\n initial_refinement_level = 4,\n periodicity = (true, false)\n)\nbc = Dict(\n :y_neg =>boundary_condition_pressure_in,\n :y_pos => Trixi.BoundaryConditionDoNothing()\n )\nvolume_flux = (flux_lax_friedrichs,flux_nonconservative)\nsurface_flux = (flux_lax_friedrichs,flux_nonconservative)\nbasis = LobattoLegendreBasis(2)\nid = IndicatorHennemannGassner(eq,basis;variable=first)\nvol =VolumeIntegralShockCapturingHG(id,volume_flux_dg = surface_flux,volume_flux_fv = volume_flux) \nsolver = DGSEM(basis,surface_flux,vol)\nsemi = SemidiscretizationHyperbolic(mesh,\neq,\ninitial_condition_simple,\nsource_terms = source_term_simple,\nsolver,\nboundary_conditions = bc)\nTrixi.default_analysis_integrals(::BloodFlowEquations2D) = ()\ntspan = (0.0, 0.3)\node = semidiscretize(semi, tspan)\nsummary_callback = SummaryCallback()\nanalysis_callback = AliveCallback(analysis_interval=1000)\nstepsize_callback = StepsizeCallback(; cfl=0.5)\ncallbacks = CallbackSet(summary_callback,analysis_callback,stepsize_callback)\ndt = stepsize_callback(ode)\nsol = solve(ode, SSPRK33(),dt=dt, dtmax = 1e-4,dtmin = 1e-12,save_everystep = false,saveat = 0.003, callback = callbacks)\n@gif for i in eachindex(sol)\npd = PlotData2D(sol[i],semi,solution_variables=cons2prim)\nplt1 = Plots.plot(pd[\"A\"],aspect_ratio=0.2)\nplt2 = Plots.plot(pd[\"wtheta\"],aspect_ratio=0.2)\nplt3 = Plots.plot(pd[\"ws\"],aspect_ratio=0.2)\nplt4 = Plots.plot(pd[\"P\"],aspect_ratio=0.2)\nplot(plt1,plt2,plt3,plt4,layout=(2,2))\nend","category":"page"},{"location":"tuto/","page":"Tutorial","title":"Tutorial","text":"(Image: Alt Text)","category":"page"},{"location":"math/#Modèles-mathématiques-1D-et-2D-pour-l’écoulement-sanguin","page":"Mathematics","title":"Modèles mathématiques 1D et 2D pour l’écoulement sanguin","text":"","category":"section"},{"location":"math/#Modèle-1D","page":"Mathematics","title":"Modèle 1D","text":"","category":"section"},{"location":"math/","page":"Mathematics","title":"Mathematics","text":"Le modèle 1D repose sur une intégration par section des équations de Navier-Stokes sous l’hypothèse d’un écoulement incompressible dans des artères supposées fines. Ce modèle est particulièrement adapté pour des études globales du réseau artériel, où la géométrie est approximativement linéaire ou faiblement courbée.","category":"page"},{"location":"math/#Hypothèses-et-simplifications","page":"Mathematics","title":"Hypothèses et simplifications","text":"","category":"section"},{"location":"math/","page":"Mathematics","title":"Mathematics","text":"L’écoulement est considéré comme incompressible.\nL’artère est modélisée comme un tube cylindrique de section variable en fonction de la pression.\nUn profil de vitesse parabolique est utilisé, permettant une moyennisation sur la section transversale de l’artère.","category":"page"},{"location":"math/#Équations-principales","page":"Mathematics","title":"Équations principales","text":"","category":"section"},{"location":"math/","page":"Mathematics","title":"Mathematics","text":"Les équations dérivées sont un système d’équations hyperboliques aux dérivées partielles décrivant la conservation de la masse et de la quantité de mouvement :","category":"page"},{"location":"math/","page":"Mathematics","title":"Mathematics","text":"Conservation de la masse :","category":"page"},{"location":"math/","page":"Mathematics","title":"Mathematics","text":"_t A + _x Q = 0","category":"page"},{"location":"math/","page":"Mathematics","title":"Mathematics","text":"Conservation de la quantité de mouvement :","category":"page"},{"location":"math/","page":"Mathematics","title":"Mathematics","text":"_t Q + _x left( alpha fracQ^2A + frac1rho A P(A x) right) = frac1rho P(A x) _x A - K fracQA","category":"page"},{"location":"math/#Énergie-et-relation-d’entropie-du-modèle-1D","page":"Mathematics","title":"Énergie et relation d’entropie du modèle 1D","text":"","category":"section"},{"location":"math/","page":"Mathematics","title":"Mathematics","text":"L’énergie associée au système est donnée par :","category":"page"},{"location":"math/","page":"Mathematics","title":"Mathematics","text":"E(t x) = fracA u_x^22 + frac1rho A P(A x) - fracbeta(x)3 rho A_0(x) A^32","category":"page"},{"location":"math/","page":"Mathematics","title":"Mathematics","text":"La relation d’entropie vérifiée par cette énergie est :","category":"page"},{"location":"math/","page":"Mathematics","title":"Mathematics","text":"_t E + _x left( left( E + fracbeta(x)3 rho A_0(x) A^32 right) u_x right) = _x left( 3 nu A _x left( fracQA right) right) u_x + frac2 pi R k1 - R k 4 nu u_x^2 0","category":"page"},{"location":"math/","page":"Mathematics","title":"Mathematics","text":"Sous des conditions aux limites nulles :","category":"page"},{"location":"math/","page":"Mathematics","title":"Mathematics","text":"_t left( int_0^L E dx right) = - 3 nu int_0^L A (_x u_x)^2 dx - frac2 pi R k1 - R k 4 nu int_0^L u_x^2 dx 0","category":"page"},{"location":"math/","page":"Mathematics","title":"Mathematics","text":"","category":"page"},{"location":"math/#Modèle-2D","page":"Mathematics","title":"Modèle 2D","text":"","category":"section"},{"location":"math/","page":"Mathematics","title":"Mathematics","text":"Le modèle 2D est dérivé à partir d’une intégration radiale des équations de Navier-Stokes, permettant de mieux représenter les effets locaux dans des configurations géométriques complexes, comme les bifurcations artérielles et les anévrismes sévères.","category":"page"},{"location":"math/#Hypothèses-et-simplifications-2","page":"Mathematics","title":"Hypothèses et simplifications","text":"","category":"section"},{"location":"math/","page":"Mathematics","title":"Mathematics","text":"L’écoulement est supposé incompressible.\nLa géométrie de l’artère est décrite à l’aide d’un système de coordonnées curvilignes (( s, \theta )).\nLe profil de vitesse est obtenu sans recourir à un ansatz spécifique.","category":"page"},{"location":"math/#Équations-principales-2","page":"Mathematics","title":"Équations principales","text":"","category":"section"},{"location":"math/","page":"Mathematics","title":"Mathematics","text":"Conservation de la masse :","category":"page"},{"location":"math/","page":"Mathematics","title":"Mathematics","text":"_t A + _θ left( fracQ_RθA right) + _s(Q_s) = 0","category":"page"},{"location":"math/","page":"Mathematics","title":"Mathematics","text":"Conservation de la quantité de mouvement (composante radiale et axiale) :","category":"page"},{"location":"math/","page":"Mathematics","title":"Mathematics","text":"_t (Q_Rθ) + _θ left( fracQ_Rθ^22 A^2 + A P right) + _s left( fracQ_Rθ Q_sA right) = frac2 R3 C sin θ fracQ_s^2A + frac2 R k Q_RθA + _θ (A P)","category":"page"},{"location":"math/","page":"Mathematics","title":"Mathematics","text":"_t (Q_s) + _θ left( fracQ_s Q_RθA^2 right) + _s left( fracQ_s^2A - fracQ_Rθ^22 A^2 + A P right) = - frac2 R3 C sin θ fracQ_Rθ Q_sA^2 + frack R Q_sA + _s (A P)","category":"page"},{"location":"math/#Énergie-et-relation-d’entropie-du-modèle-2D","page":"Mathematics","title":"Énergie et relation d’entropie du modèle 2D","text":"","category":"section"},{"location":"math/","page":"Mathematics","title":"Mathematics","text":"L’énergie associée au système est donnée par :","category":"page"},{"location":"math/","page":"Mathematics","title":"Mathematics","text":"E(t θ s) = A left( frac98 u_θ^2 + fracu_s^22 + p right) - p","category":"page"},{"location":"math/","page":"Mathematics","title":"Mathematics","text":"La relation d’entropie correspondante est :","category":"page"},{"location":"math/","page":"Mathematics","title":"Mathematics","text":"_t E + _θ left( frac32 fracu_θR left( E + p - frac916 A u_θ^2 right) right) + _s left( u_s left( E + p - frac916 A u_θ^2 right) right) = frac94 R k u_θ^2 + k R u_s^2 0","category":"page"},{"location":"math/","page":"Mathematics","title":"Mathematics","text":"Cette relation garantit que l’énergie décroît localement dans le temps, ce qui assure la stabilité du modèle.","category":"page"},{"location":"math/","page":"Mathematics","title":"Mathematics","text":"","category":"page"},{"location":"math/#Comparaison-des-modèles-1D-et-2D","page":"Mathematics","title":"Comparaison des modèles 1D et 2D","text":"","category":"section"},{"location":"math/","page":"Mathematics","title":"Mathematics","text":"Modèle 1D :\nRapide et efficace pour des simulations globales sur de grands réseaux artériels.\nBien adapté pour des géométries simples ou faiblement courbées.\nCoût de calcul très faible.\nModèle 2D :\nPlus précis pour des géométries complexes (bifurcations, anévrismes).\nPermet de mieux capturer les effets locaux et les interactions fluide-structure.\nCoût de calcul modéré par rapport aux modèles tridimensionnels (NS-FSI 3D).","category":"page"},{"location":"math/","page":"Mathematics","title":"Mathematics","text":"L’utilisation combinée de ces deux modèles permet une alternative efficace aux simulations 3D, tout en offrant un bon compromis entre précision et coût de calcul.","category":"page"},{"location":"","page":"Home","title":"Home","text":"CurrentModule = BloodFlowTrixi","category":"page"},{"location":"#BloodFlowTrixi.jl","page":"Home","title":"BloodFlowTrixi.jl","text":"","category":"section"},{"location":"","page":"Home","title":"Home","text":"BloodFlowTrixi.jl is a Julia package that implements one-dimensional (1D) and two-dimensional (2D) blood flow models for arterial circulation. These models are derived from the Navier-Stokes equations and were developed as part of my PhD research in applied mathematics, focusing on cardiovascular pathologies such as aneurysms and stenoses.","category":"page"},{"location":"#Description","page":"Home","title":"Description","text":"","category":"section"},{"location":"","page":"Home","title":"Home","text":"This package provides:","category":"page"},{"location":"","page":"Home","title":"Home","text":"1D Blood Flow Model: This model describes blood flow along compliant arteries in a single spatial dimension. It was derived under the assumption of axisymmetric flow and accounts for arterial compliance, inertia, and frictional losses.\nMore details about this model can be found in my corresponding publication: [Article 1D]","category":"page"},{"location":"","page":"Home","title":"Home","text":"leftbeginaligned\n fracpartial apartial t + fracpartialpartial x(Q) = 0 \n fracpartial Qpartial t + fracpartialpartial xleft(fracQ^2A + A P(a)right) = P(a) fracpartial Apartial x - 2 pi R k frac Q A\n P(a) = P_ext + fracEhsqrtpi1-xi^2fracsqrtA - sqrtA_0A_0 \n R = sqrtfracApi\nendalignedright","category":"page"},{"location":"","page":"Home","title":"Home","text":"2D Blood Flow Model: The 2D model extends the Navier-Stokes equations under the thin-artery assumption, allowing for simulations in complex arterial geometries using curvilinear coordinates. It captures both longitudinal and angular dynamics, making it more accurate than classical 1D models while being less computationally expensive than full 3D models.\nThis model is described in detail in: [Article 2D]","category":"page"},{"location":"","page":"Home","title":"Home","text":"leftbeginaligned\n fracpartial apartial t + fracpartialpartial thetaleft( fracQ_RthetaA right) + fracpartialpartial s(Q_s) = 0 \n fracpartial Q_Rthetapartial t + fracpartialpartial thetaleft(fracQ_Rtheta^22A^2 + A P(a)right) + fracpartialpartial sleft( fracQ_RthetaQ_sA right) = P(a) fracpartial Apartial theta - 2 R k fracQ_RthetaA + frac2R3 mathcalCsin theta fracQ_s^2A \n fracpartial Q_spartial t + fracpartialpartial thetaleft(fracQ_Rtheta Q_sA^2 right) + fracpartialpartial sleft( fracQ_s^2A - fracQ_Rtheta^22A^2 + A P(a) right) = P(a) fracpartial Apartial s - R k fracQ_sA - frac2R3 mathcalCsin theta fracQ_s Q_RthetaA^2 \n P(a) = P_ext + fracEhsqrt2left(1-xi^2right)fracsqrtA - sqrtA_0A_0 \n R = sqrt2A\nendalignedright","category":"page"},{"location":"","page":"Home","title":"Home","text":"Both models were designed to be used with Trixi.jl, a flexible and high-performance framework for solving systems of conservation laws using the Discontinuous Galerkin (DG) method.","category":"page"},{"location":"#Features","page":"Home","title":"Features","text":"","category":"section"},{"location":"","page":"Home","title":"Home","text":"1D and 2D models for arterial blood flow.\nDerived from the Navier-Stokes equations with appropriate assumptions for compliant arteries.\nTo be used with Trixi.jl for DG-based numerical simulations.\nSupport for curvilinear geometries and compliant wall dynamics.","category":"page"},{"location":"#Installation","page":"Home","title":"Installation","text":"","category":"section"},{"location":"","page":"Home","title":"Home","text":"To install BloodFlowTrixi.jl, use the following commands in Julia:","category":"page"},{"location":"","page":"Home","title":"Home","text":"julia> ]\npkg> add Trixi\npkg> add https://github.com/your-repo/BloodFlowTrixi.jl","category":"page"},{"location":"#Future-Plans","page":"Home","title":"Future Plans","text":"","category":"section"},{"location":"","page":"Home","title":"Home","text":"short term","category":"page"},{"location":"","page":"Home","title":"Home","text":"Add second order 1D model.\nDesign prim variables for 1D and 2D models.\nAdd proper tests for 1D and 2D models.\nAdd 3D representations of the solutions for 1D and 2D models.\nDesign easy to use interfaces for users to define their own initial and boundary conditions and source terms.","category":"page"},{"location":"","page":"Home","title":"Home","text":"long term","category":"page"},{"location":"","page":"Home","title":"Home","text":"Add 3D fluid-structure interaction models for complex arterial geometries.\nDesign support for artery networks and simulate vascular networks using the 2D and 1D model.\nAutodiff support for 1D and 2D models for parameter optimization.","category":"page"},{"location":"#License","page":"Home","title":"License","text":"","category":"section"},{"location":"","page":"Home","title":"Home","text":"This package is licensed under the MIT license.","category":"page"},{"location":"#Acknowledgments","page":"Home","title":"Acknowledgments","text":"","category":"section"},{"location":"","page":"Home","title":"Home","text":"This package was developed as part of my PhD research in applied mathematics, focusing on mathematical modeling and numerical simulation of blood flow in arteries. Special thanks to the developers of Trixi.jl, whose framework was invaluable in implementing and testing these models.","category":"page"},{"location":"","page":"Home","title":"Home","text":"","category":"page"},{"location":"","page":"Home","title":"Home","text":"Modules = [BloodFlowTrixi]","category":"page"},{"location":"#BloodFlowTrixi.BloodFlowTrixi","page":"Home","title":"BloodFlowTrixi.BloodFlowTrixi","text":"Package BloodFlowTrixi v0.0.1\n\nThis package implements 1D and 2D blood flow models for arterial circulation using Trixi.jl, enabling efficient numerical simulation and analysis.\n\nDocs under https://yolhan83.github.io/BloodFlowTrixi.jl\n\n\n\n\n\n","category":"module"},{"location":"#BloodFlowTrixi.BloodFlowEquations1D","page":"Home","title":"BloodFlowTrixi.BloodFlowEquations1D","text":"BloodFlowEquations1D(;h,rho=1.0,xi=0.25,nu=0.04)\n\nBlood Flow equations in one space dimension. This model describes the dynamics of blood flow along a compliant artery using one-dimensional equations derived from the Navier-Stokes equations. The equations account for conservation of mass and momentum, incorporating the effect of arterial compliance and frictional losses.\n\nThe governing equations are given by\n\nleftbeginaligned\n fracpartial apartial t + fracpartialpartial x(Q) = 0 \n fracpartial Qpartial t + fracpartialpartial xleft(fracQ^2A + A P(a)right) = P(a) fracpartial Apartial x - 2 pi R k frac Q A\n P(a) = P_ext + fracEhsqrtpi1-xi^2fracsqrtA - sqrtA_0A_0 \n R = sqrtfracApi\nendalignedright\n\n\n\n\n\n","category":"type"},{"location":"#BloodFlowTrixi.BloodFlowEquations2D","page":"Home","title":"BloodFlowTrixi.BloodFlowEquations2D","text":"BloodFlowEquations2D(;h,rho=1.0,xi=0.25)\n\nDefines the two-dimensional blood flow equations derived from the Navier-Stokes equations in curvilinear coordinates under the thin-artery assumption. This model describes the dynamics of blood flow along a compliant artery in two spatial dimensions (s, θ).\n\nParameters\n\nh::T: Wall thickness of the artery.\nrho::T: Fluid density (default 1.0).\nxi::T: Poisson's ratio (default 0.25).\nnu::T: Viscosity coefficient.\n\nThe governing equations account for conservation of mass and momentum, incorporating the effects of arterial compliance, curvature, and frictional losses.\n\nleftbeginaligned\n fracpartial apartial t + fracpartialpartial thetaleft( fracQ_RthetaA right) + fracpartialpartial s(Q_s) = 0 \n fracpartial Q_Rthetapartial t + fracpartialpartial thetaleft(fracQ_Rtheta^22A^2 + A P(a)right) + fracpartialpartial sleft( fracQ_RthetaQ_sA right) = P(a) fracpartial Apartial theta - 2 R k fracQ_RthetaA + frac2R3 mathcalCsin theta fracQ_s^2A \n fracpartial Q_spartial t + fracpartialpartial thetaleft(fracQ_Rtheta Q_sA^2 right) + fracpartialpartial sleft( fracQ_s^2A - fracQ_Rtheta^22A^2 + A P(a) right) = P(a) fracpartial Apartial s - R k fracQ_sA - frac2R3 mathcalCsin theta fracQ_s Q_RthetaA^2 \n P(a) = P_ext + fracEhsqrt2left(1-xi^2right)fracsqrtA - sqrtA_0A_0 \n R = sqrt2A\nendalignedright\n\n\n\n\n\n","category":"type"},{"location":"#Trixi.DissipationLocalLaxFriedrichs-Tuple{Any, Any, Any, BloodFlowEquations1D}","page":"Home","title":"Trixi.DissipationLocalLaxFriedrichs","text":"(dissipation::Trixi.DissipationLocalLaxFriedrichs)(u_ll, u_rr, orientation_or_normal_direction, eq::BloodFlowEquations1D)\n\nCalculates the dissipation term using the Local Lax-Friedrichs method.\n\nParameters\n\nu_ll: Left state vector.\nu_rr: Right state vector.\norientation_or_normal_direction: Orientation or normal direction.\neq: Instance of BloodFlowEquations1D.\n\nReturns\n\nDissipation vector.\n\n\n\n\n\n","category":"method"},{"location":"#BloodFlowTrixi.boundary_condition_outflow-Tuple{Any, Any, Any, Any, Any, Any, BloodFlowEquations1D}","page":"Home","title":"BloodFlowTrixi.boundary_condition_outflow","text":"boundary_condition_outflow(u_inner, orientation_or_normal, direction, x, t, surface_flux_function, eq::BloodFlowEquations1D)\n\nImplements the outflow boundary condition, assuming that there is no reflection at the boundary.\n\nParameters\n\nu_inner: State vector inside the domain near the boundary.\norientation_or_normal: Normal orientation of the boundary.\ndirection: Integer indicating the direction of the boundary.\nx: Position vector.\nt: Time.\nsurface_flux_function: Function to compute flux at the boundary.\neq: Instance of BloodFlowEquations1D.\n\nReturns\n\nComputed boundary flux.\n\n\n\n\n\n","category":"method"},{"location":"#BloodFlowTrixi.boundary_condition_pressure_in-Tuple{Any, Any, Any, Any, Any, Any, BloodFlowEquations1D}","page":"Home","title":"BloodFlowTrixi.boundary_condition_pressure_in","text":"boundary_condition_pressure_in(u_inner, orientation_or_normal, direction, x, t, surface_flux_function, eq::BloodFlowEquations1D)\n\nImplements a pressure inflow boundary condition where the inflow pressure varies with time.\n\nParameters\n\nu_inner: State vector inside the domain near the boundary.\norientation_or_normal: Normal orientation of the boundary.\ndirection: Integer indicating the boundary direction.\nx: Position vector.\nt: Time scalar.\nsurface_flux_function: Function to compute flux at the boundary.\neq: Instance of BloodFlowEquations1D.\n\nReturns\n\nComputed boundary flux with inflow pressure specified by:\n\nP_in = begincases\n2 times 10^4 sin^2(pi t 0125) textif t 0125 \n0 textotherwise\nendcases\n\nThe corresponding inflow area A_{in} is computed using the inverse pressure relation, and the boundary state is constructed accordingly.\n\n\n\n\n\n","category":"method"},{"location":"#BloodFlowTrixi.boundary_condition_slip_wall-Tuple{Any, Any, Any, Any, Any, Any, BloodFlowEquations1D}","page":"Home","title":"BloodFlowTrixi.boundary_condition_slip_wall","text":"boundary_condition_slip_wall(u_inner, orientation_or_normal, direction, x, t, surface_flux_function, eq::BloodFlowEquations1D)\n\nImplements a slip wall boundary condition where the normal component of velocity is reflected.\n\nParameters\n\nu_inner: State vector inside the domain near the boundary.\norientation_or_normal: Normal orientation of the boundary.\ndirection: Integer indicating the direction of the boundary.\nx: Position vector.\nt: Time.\nsurface_flux_function: Function to compute flux at the boundary.\neq: Instance of BloodFlowEquations1D.\n\nReturns\n\nComputed boundary flux at the slip wall.\n\n\n\n\n\n","category":"method"},{"location":"#BloodFlowTrixi.flux_nonconservative-Tuple{Any, Any, Integer, BloodFlowEquations1D}","page":"Home","title":"BloodFlowTrixi.flux_nonconservative","text":"flux_nonconservative(u_ll, u_rr, orientation::Integer, eq::BloodFlowEquations1D)\n\nComputes the non-conservative flux for the model, used for handling discontinuities in pressure.\n\nParameters\n\nu_ll: Left state vector.\nu_rr: Right state vector.\norientation::Integer: Orientation index.\neq: Instance of BloodFlowEquations1D.\n\nReturns\n\nNon-conservative flux vector.\n\n\n\n\n\n","category":"method"},{"location":"#BloodFlowTrixi.friction-Tuple{Any, Any, BloodFlowEquations1D}","page":"Home","title":"BloodFlowTrixi.friction","text":"friction(u, x, eq::BloodFlowEquations1D)\n\nCalculates the friction term for the blood flow equations, which represents viscous resistance to flow along the artery wall.\n\nParameters\n\nu: State vector containing cross-sectional area and flow rate.\nx: Position along the artery.\neq::BloodFlowEquations1D: Instance of the blood flow model.\n\nReturns\n\nFriction coefficient as a scalar.\n\n\n\n\n\n","category":"method"},{"location":"#BloodFlowTrixi.initial_condition_simple-Tuple{Any, Any, BloodFlowEquations1D}","page":"Home","title":"BloodFlowTrixi.initial_condition_simple","text":"initial_condition_simple(x, t, eq::BloodFlowEquations1D; R0=2.0)\n\nGenerates a simple initial condition with a specified initial radius R0.\n\nParameters\n\nx: Position vector.\nt: Time scalar.\neq::BloodFlowEquations1D: Instance of the blood flow model.\nR0: Initial radius (default: 2.0).\n\nReturns\n\nState vector with zero initial area perturbation, zero flow rate, constant elasticity modulus, and reference area computed as A_0 = \\pi R_0^2.\n\nThis initial condition is suitable for basic tests without complex dynamics.\n\n\n\n\n\n","category":"method"},{"location":"#BloodFlowTrixi.inv_pressure-Tuple{Any, Any, BloodFlowEquations1D}","page":"Home","title":"BloodFlowTrixi.inv_pressure","text":"inv_pressure(p, u, eq::BloodFlowEquations1D)\n\nComputes the inverse relation of pressure to cross-sectional area.\n\nParameters\n\np: Pressure.\nu: State vector.\neq: Instance of BloodFlowEquations1D.\n\nReturns\n\nCross-sectional area corresponding to the given pressure.\n\n\n\n\n\n","category":"method"},{"location":"#BloodFlowTrixi.pressure-Tuple{Any, BloodFlowEquations1D}","page":"Home","title":"BloodFlowTrixi.pressure","text":"pressure(u, eq::BloodFlowEquations1D)\n\nComputes the pressure given the state vector based on the compliance of the artery.\n\nParameters\n\nu: State vector.\neq: Instance of BloodFlowEquations1D.\n\nReturns\n\nPressure as a scalar.\n\n\n\n\n\n","category":"method"},{"location":"#BloodFlowTrixi.pressure_der-Tuple{Any, BloodFlowEquations1D}","page":"Home","title":"BloodFlowTrixi.pressure_der","text":"pressure_der(u, eq::BloodFlowEquations1D)\n\nComputes the derivative of pressure with respect to cross-sectional area.\n\nParameters\n\nu: State vector.\neq: Instance of BloodFlowEquations1D.\n\nReturns\n\nDerivative of pressure.\n\n\n\n\n\n","category":"method"},{"location":"#BloodFlowTrixi.radius-Tuple{Any, BloodFlowEquations1D}","page":"Home","title":"BloodFlowTrixi.radius","text":"radius(u, eq::BloodFlowEquations1D)\n\nComputes the radius of the artery based on the cross-sectional area.\n\nParameters\n\nu: State vector.\neq: Instance of BloodFlowEquations1D.\n\nReturns\n\nRadius as a scalar.\n\n\n\n\n\n","category":"method"},{"location":"#BloodFlowTrixi.source_term_simple-Tuple{Any, Any, Any, BloodFlowEquations1D}","page":"Home","title":"BloodFlowTrixi.source_term_simple","text":"source_term_simple(u, x, t, eq::BloodFlowEquations1D)\n\nComputes a simple source term for the blood flow model, focusing on frictional effects.\n\nParameters\n\nu: State vector containing area perturbation, flow rate, elasticity modulus, and reference area.\nx: Position vector.\nt: Time scalar.\neq::BloodFlowEquations1D: Instance of the blood flow model.\n\nReturns\n\nSource terms vector where:\n\ns_1 = 0 (no source for area perturbation).\ns_2 represents the friction term given by s_2 = \\frac{2 \\pi k Q}{R A}.\n\nFriction coefficient k is computed using the friction function, and the radius R is obtained using the radius function.\n\n\n\n\n\n","category":"method"},{"location":"#Trixi.cons2prim-Tuple{Any, BloodFlowEquations1D}","page":"Home","title":"Trixi.cons2prim","text":"Trixi.cons2prim(u, eq::BloodFlowEquations1D)\n\nConverts the conserved variables to primitive variables.\n\nParameters\n\nu: State vector.\neq: Instance of BloodFlowEquations1D.\n\nReturns\n\nPrimitive variable vector.\n\n\n\n\n\n","category":"method"},{"location":"#Trixi.flux-Tuple{Any, Integer, BloodFlowEquations1D}","page":"Home","title":"Trixi.flux","text":"Trixi.flux(u, orientation::Integer, eq::BloodFlowEquations1D)\n\nComputes the flux vector for the conservation laws of the blood flow model.\n\nParameters\n\nu: State vector.\norientation::Integer: Orientation index for flux computation.\neq: Instance of BloodFlowEquations1D.\n\nReturns\n\nFlux vector as an SVector.\n\n\n\n\n\n","category":"method"},{"location":"#Trixi.initial_condition_convergence_test-Tuple{Any, Any, BloodFlowEquations1D}","page":"Home","title":"Trixi.initial_condition_convergence_test","text":"initial_condition_convergence_test(x, t, eq::BloodFlowEquations1D)\n\nGenerates a smooth initial condition for convergence tests of the blood flow equations.\n\nParameters\n\nx: Position vector.\nt: Time scalar.\neq::BloodFlowEquations1D: Instance of the blood flow model.\n\nReturns\n\nInitial condition state vector with zero initial area perturbation, sinusoidal flow rate, a constant elasticity modulus, and reference area.\n\nDetails\n\nThe returned initial condition has:\n\nZero perturbation in area (a = 0).\nA sinusoidal flow rate given by Q = sin(\\pi x t).\nA constant elasticity modulus E.\nA reference cross-sectional area A_0 = \\pi R_0^2 for R_0 = 1.\n\nThis initial condition can be used to verify the accuracy and stability of numerical solvers.\n\n\n\n\n\n","category":"method"},{"location":"#Trixi.max_abs_speed_naive-Tuple{Any, Any, Integer, BloodFlowEquations1D}","page":"Home","title":"Trixi.max_abs_speed_naive","text":"Trixi.max_abs_speed_naive(u_ll, u_rr, orientation::Integer, eq::BloodFlowEquations1D)\n\nCalculates the maximum absolute speed for wave propagation in the blood flow model using a naive approach.\n\nParameters\n\nu_ll: Left state vector.\nu_rr: Right state vector.\norientation::Integer: Orientation index.\neq: Instance of BloodFlowEquations1D.\n\nReturns\n\nMaximum absolute speed.\n\n\n\n\n\n","category":"method"},{"location":"#Trixi.prim2cons-Tuple{Any, BloodFlowEquations1D}","page":"Home","title":"Trixi.prim2cons","text":"Trixi.prim2cons(u, eq::BloodFlowEquations1D)\n\nConverts the primitive variables to conserved variables.\n\nParameters\n\nu: Primitive variable vector.\neq: Instance of BloodFlowEquations1D.\n\nReturns\n\nConserved variable vector.\n\n\n\n\n\n","category":"method"},{"location":"#Trixi.source_terms_convergence_test-Tuple{Any, Any, Any, BloodFlowEquations1D}","page":"Home","title":"Trixi.source_terms_convergence_test","text":"source_terms_convergence_test(u, x, t, eq::BloodFlowEquations1D)\n\nComputes the source terms for convergence tests of the blood flow equations.\n\nParameters\n\nu: State vector containing area perturbation, flow rate, elasticity modulus, and reference area.\nx: Position vector.\nt: Time scalar.\neq::BloodFlowEquations1D: Instance of the blood flow model.\n\nReturns\n\nSource terms vector.\n\nDetails\n\nThe source terms are derived based on the smooth initial condition and friction effects:\n\ns_1 represents the source term for area perturbation and is given by s_1 = \\pi t \\cos(\\pi x t).\ns_2 represents the source term for the flow rate and includes contributions from spatial and temporal variations as well as friction effects.\n\nThe radius R is computed using the radius function, and the friction coefficient k is obtained using the friction function.\n\nThis function is useful for evaluating the correctness of source term handling in numerical solvers.\n\n\n\n\n\n","category":"method"}] +[{"location":"tuto/#Tutorial-for-the-1D-model","page":"Tutorial","title":"Tutorial for the 1D model","text":"","category":"section"},{"location":"tuto/","page":"Tutorial","title":"Tutorial","text":"In this section, we describe how to use BloodFlowTrixi.jl with Trixi.jl. This tutorial will guide you through setting up and running a 1D blood flow simulation, including mesh creation, boundary conditions, numerical fluxes, and visualization of results.","category":"page"},{"location":"tuto/#Packages","page":"Tutorial","title":"Packages","text":"","category":"section"},{"location":"tuto/","page":"Tutorial","title":"Tutorial","text":"Before starting, ensure that the required packages are loaded:","category":"page"},{"location":"tuto/","page":"Tutorial","title":"Tutorial","text":"using Trixi\nusing BloodFlowTrixi\nusing OrdinaryDiffEq\nusing Plots","category":"page"},{"location":"tuto/","page":"Tutorial","title":"Tutorial","text":"First, we need to choose the equation that describes the blood flow dynamics:","category":"page"},{"location":"tuto/","page":"Tutorial","title":"Tutorial","text":"eq = BloodFlowEquations1D(; h=0.1)","category":"page"},{"location":"tuto/","page":"Tutorial","title":"Tutorial","text":"Here, h represents a parameter related to the initial condition or model scaling.","category":"page"},{"location":"tuto/#Mesh-and-boundary-conditions","page":"Tutorial","title":"Mesh and boundary conditions","text":"","category":"section"},{"location":"tuto/","page":"Tutorial","title":"Tutorial","text":"We begin by defining a one-dimensional Tree mesh, which discretizes the spatial domain:","category":"page"},{"location":"tuto/","page":"Tutorial","title":"Tutorial","text":"mesh = TreeMesh(0.0, 40.0, initial_refinement_level=6, n_cells_max=10^4, periodicity=false)","category":"page"},{"location":"tuto/","page":"Tutorial","title":"Tutorial","text":"This generates a non-periodic mesh for the interval 0 40, with 2^initial-refinement-level+1-1 cells. The parameter initial_refinement_level controls the initial number of cells, while n_cells_max specifies the maximum number of cells allowed during mesh refinement.","category":"page"},{"location":"tuto/","page":"Tutorial","title":"Tutorial","text":"In Trixi.jl, the Tree mesh has two labeled boundaries: ***xneg*** (left boundary) and ***xpos*** (right boundary). These labels are used to apply boundary conditions:","category":"page"},{"location":"tuto/","page":"Tutorial","title":"Tutorial","text":"bc = (\n x_neg = boundary_condition_pressure_in,\n x_pos = Trixi.BoundaryConditionDoNothing()\n)","category":"page"},{"location":"tuto/","page":"Tutorial","title":"Tutorial","text":"boundary_condition_pressure_in applies a pressure inflow condition at the left boundary.\nTrixi.BoundaryConditionDoNothing() specifies a \"do nothing\" boundary condition at the right boundary, meaning no flux is imposed.","category":"page"},{"location":"tuto/#Boundary-condition-implementation","page":"Tutorial","title":"Boundary condition implementation","text":"","category":"section"},{"location":"tuto/","page":"Tutorial","title":"Tutorial","text":"The inflow boundary condition is defined as:","category":"page"},{"location":"tuto/","page":"Tutorial","title":"Tutorial","text":"boundary_condition_pressure_in(u_inner, orientation_or_normal, direction, x, t, surface_flux_function, eq::BloodFlowEquations1D)","category":"page"},{"location":"tuto/","page":"Tutorial","title":"Tutorial","text":"This function applies a time-dependent pressure inflow condition.","category":"page"},{"location":"tuto/#Parameters","page":"Tutorial","title":"Parameters","text":"","category":"section"},{"location":"tuto/","page":"Tutorial","title":"Tutorial","text":"u_inner: State vector inside the domain near the boundary.\norientation_or_normal: Normal orientation of the boundary.\ndirection: Integer indicating the boundary direction.\nx: Position vector.\nt: Time scalar.\nsurface_flux_function: Function to compute flux at the boundary.\neq: Instance of BloodFlowEquations1D.","category":"page"},{"location":"tuto/#Returns","page":"Tutorial","title":"Returns","text":"","category":"section"},{"location":"tuto/","page":"Tutorial","title":"Tutorial","text":"The boundary flux is computed based on the inflow pressure: $ P{\\text{in}} = \\begin{cases} 2 \\times 10^4 \\sin^2\\left(\\frac{\\pi t}{0.125}\\right) & \\text{if } t < 0.125 \\\n0 & \\text{otherwise} \\end{cases} $ This time-dependent inflow pressure mimics a pulsatile flow, typical in arterial blood flow. The inflow area A{\\text{in}}$ is determined using the inverse pressure relation, ensuring consistency with the physical model.","category":"page"},{"location":"tuto/#Numerical-flux","page":"Tutorial","title":"Numerical flux","text":"","category":"section"},{"location":"tuto/","page":"Tutorial","title":"Tutorial","text":"To compute fluxes at cell interfaces, we use a combination of conservative and non-conservative fluxes:","category":"page"},{"location":"tuto/","page":"Tutorial","title":"Tutorial","text":"volume_flux = (flux_lax_friedrichs, flux_nonconservative)\nsurface_flux = (flux_lax_friedrichs, flux_nonconservative)","category":"page"},{"location":"tuto/","page":"Tutorial","title":"Tutorial","text":"flux_lax_friedrichs is a standard numerical flux for hyperbolic conservation laws.\nflux_nonconservative handles the non-conservative terms in the model, particularly those related to pressure discontinuities.","category":"page"},{"location":"tuto/","page":"Tutorial","title":"Tutorial","text":"The non-conservative flux function is defined as:","category":"page"},{"location":"tuto/","page":"Tutorial","title":"Tutorial","text":"flux_nonconservative(u_ll, u_rr, orientation::Integer, eq::BloodFlowEquations1D)","category":"page"},{"location":"tuto/#Parameters-2","page":"Tutorial","title":"Parameters","text":"","category":"section"},{"location":"tuto/","page":"Tutorial","title":"Tutorial","text":"u_ll: Left state vector.\nu_rr: Right state vector.\norientation::Integer: Orientation index.\neq: Instance of BloodFlowEquations1D.","category":"page"},{"location":"tuto/#Returns-2","page":"Tutorial","title":"Returns","text":"","category":"section"},{"location":"tuto/","page":"Tutorial","title":"Tutorial","text":"The function returns the non-conservative flux vector, which is essential for capturing sharp pressure changes in the simulation.","category":"page"},{"location":"tuto/#Basis-functions-and-Shock-Capturing-DG-scheme","page":"Tutorial","title":"Basis functions and Shock Capturing DG scheme","text":"","category":"section"},{"location":"tuto/","page":"Tutorial","title":"Tutorial","text":"To approximate the solution, we use polynomial basis functions:","category":"page"},{"location":"tuto/","page":"Tutorial","title":"Tutorial","text":"basis = LobattoLegendreBasis(2)","category":"page"},{"location":"tuto/","page":"Tutorial","title":"Tutorial","text":"This defines a Lobatto-Legendre basis of polynomial degree 2, which is commonly used in high-order methods like Discontinuous Galerkin (DG) schemes.","category":"page"},{"location":"tuto/","page":"Tutorial","title":"Tutorial","text":"We then define an indicator for shock capturing, focusing on the first variable (area perturbation $a$):","category":"page"},{"location":"tuto/","page":"Tutorial","title":"Tutorial","text":"id = IndicatorHennemannGassner(eq, basis; variable=first)","category":"page"},{"location":"tuto/","page":"Tutorial","title":"Tutorial","text":"This indicator helps detect shocks or discontinuities in the solution and applies appropriate stabilization.","category":"page"},{"location":"tuto/","page":"Tutorial","title":"Tutorial","text":"The solver is defined as:","category":"page"},{"location":"tuto/","page":"Tutorial","title":"Tutorial","text":"vol = VolumeIntegralShockCapturingHG(id, volume_flux_dg=volume_flux, volume_flux_fv=surface_flux)\nsolver = DGSEM(basis, surface_flux, vol)","category":"page"},{"location":"tuto/","page":"Tutorial","title":"Tutorial","text":"Here, DGSEM represents the Discontinuous Galerkin Spectral Element Method, a high-order accurate scheme suitable for hyperbolic problems.","category":"page"},{"location":"tuto/#Semi-discretization","page":"Tutorial","title":"Semi-discretization","text":"","category":"section"},{"location":"tuto/","page":"Tutorial","title":"Tutorial","text":"We are now ready to semi-discretize the problem:","category":"page"},{"location":"tuto/","page":"Tutorial","title":"Tutorial","text":"semi = SemidiscretizationHyperbolic(\n mesh,\n eq,\n initial_condition_simple,\n source_terms = source_term_simple,\n solver,\n boundary_conditions=bc\n)","category":"page"},{"location":"tuto/","page":"Tutorial","title":"Tutorial","text":"This step sets up the semi-discretized form of the PDE, which will be advanced in time using an ODE solver.","category":"page"},{"location":"tuto/#Source-term","page":"Tutorial","title":"Source term","text":"","category":"section"},{"location":"tuto/","page":"Tutorial","title":"Tutorial","text":"The source term accounts for additional forces acting on the blood flow, such as friction:","category":"page"},{"location":"tuto/","page":"Tutorial","title":"Tutorial","text":"source_term_simple(u, x, t, eq::BloodFlowEquations1D)","category":"page"},{"location":"tuto/#Parameters-3","page":"Tutorial","title":"Parameters","text":"","category":"section"},{"location":"tuto/","page":"Tutorial","title":"Tutorial","text":"u: State vector containing area perturbation, flow rate, elasticity modulus, and reference area.\nx: Position vector.\nt: Time scalar.\neq::BloodFlowEquations1D: Instance of the blood flow model.","category":"page"},{"location":"tuto/#Returns-3","page":"Tutorial","title":"Returns","text":"","category":"section"},{"location":"tuto/","page":"Tutorial","title":"Tutorial","text":"The source term vector is given by:","category":"page"},{"location":"tuto/","page":"Tutorial","title":"Tutorial","text":"s_1 = 0\n(no source for area perturbation).\ns_2 = frac2 pi k QR A\n, representing frictional effects.","category":"page"},{"location":"tuto/","page":"Tutorial","title":"Tutorial","text":"The friction coefficient k is computed using a model-specific friction function, and the radius R is obtained from the state vector using the radius function.","category":"page"},{"location":"tuto/#Initial-condition","page":"Tutorial","title":"Initial condition","text":"","category":"section"},{"location":"tuto/","page":"Tutorial","title":"Tutorial","text":"The initial condition specifies the starting state of the simulation:","category":"page"},{"location":"tuto/","page":"Tutorial","title":"Tutorial","text":"initial_condition_simple(x, t, eq::BloodFlowEquations1D; R0=2.0)","category":"page"},{"location":"tuto/","page":"Tutorial","title":"Tutorial","text":"This function generates a simple initial condition with a uniform radius R0.","category":"page"},{"location":"tuto/#Parameters-4","page":"Tutorial","title":"Parameters","text":"","category":"section"},{"location":"tuto/","page":"Tutorial","title":"Tutorial","text":"x: Position vector.\nt: Time scalar.\neq::BloodFlowEquations1D: Instance of the blood flow model.\nR0: Initial radius (default: 2.0).","category":"page"},{"location":"tuto/#Returns-4","page":"Tutorial","title":"Returns","text":"","category":"section"},{"location":"tuto/","page":"Tutorial","title":"Tutorial","text":"The function returns a state vector with:","category":"page"},{"location":"tuto/","page":"Tutorial","title":"Tutorial","text":"Zero initial area perturbation.\nZero initial flow rate.\nConstant elasticity modulus.\nReference area A_0 = pi R_0^2.","category":"page"},{"location":"tuto/","page":"Tutorial","title":"Tutorial","text":"This simple initial condition is suitable for testing the model without introducing complex dynamics.","category":"page"},{"location":"tuto/#Run-the-simulation","page":"Tutorial","title":"Run the simulation","text":"","category":"section"},{"location":"tuto/","page":"Tutorial","title":"Tutorial","text":"First, we discretize the problem in time:","category":"page"},{"location":"tuto/","page":"Tutorial","title":"Tutorial","text":"Trixi.default_analysis_integrals(::BloodFlowEquations1D) = ()\ntspan = (0.0, 0.5)\node = semidiscretize(semi, tspan)","category":"page"},{"location":"tuto/","page":"Tutorial","title":"Tutorial","text":"Here, tspan defines the time interval for the simulation.","category":"page"},{"location":"tuto/","page":"Tutorial","title":"Tutorial","text":"Next, we add some callbacks to monitor the simulation:","category":"page"},{"location":"tuto/","page":"Tutorial","title":"Tutorial","text":"summary_callback = SummaryCallback()\nanalysis_callback = AnalysisCallback(semi, interval=200)\nstepsize_callback = StepsizeCallback(; cfl=0.5)\ncallbacks = CallbackSet(summary_callback, analysis_callback, stepsize_callback)","category":"page"},{"location":"tuto/","page":"Tutorial","title":"Tutorial","text":"SummaryCallback provides a summary of the simulation progress.\nAnalysisCallback computes analysis metrics at specified intervals.\nStepsizeCallback adjusts the time step based on the CFL condition.","category":"page"},{"location":"tuto/","page":"Tutorial","title":"Tutorial","text":"Finally, we solve the problem:","category":"page"},{"location":"tuto/","page":"Tutorial","title":"Tutorial","text":"dt = stepsize_callback(ode)\nsol = solve(ode, SSPRK33(), dt=dt, dtmax=1e-4, dtmin=1e-11,\n save_everystep=false, saveat=0.002, callback=callbacks)","category":"page"},{"location":"tuto/","page":"Tutorial","title":"Tutorial","text":"Here, SSPRK33() is a third-order Strong Stability Preserving Runge-Kutta method, suitable for hyperbolic PDEs.","category":"page"},{"location":"tuto/#Plot-the-results","page":"Tutorial","title":"Plot the results","text":"","category":"section"},{"location":"tuto/","page":"Tutorial","title":"Tutorial","text":"The results can be visualized using the following code:","category":"page"},{"location":"tuto/","page":"Tutorial","title":"Tutorial","text":"@gif for i in eachindex(sol)\n a1 = sol[i][1:4:end]\n Q1 = sol[i][2:4:end]\n A01 = sol[i][4:4:end]\n A1 = A01 .+ a1\n plot(Q1 ./ A1, lw=4, color=:red, ylim=(-10, 50), label=\"velocity\", legend=:bottomleft)\nend","category":"page"},{"location":"tuto/","page":"Tutorial","title":"Tutorial","text":"This code generates an animated GIF showing the evolution of the velocity profile over time. The velocity is computed as QA, where Q is the flow rate, and A is the cross-sectional area.","category":"page"},{"location":"tuto/#Plain-code","page":"Tutorial","title":"Plain code","text":"","category":"section"},{"location":"tuto/","page":"Tutorial","title":"Tutorial","text":"using Trixi\nusing BloodFlowTrixi\nusing OrdinaryDiffEq,Plots\neq = BloodFlowEquations1D(;h=0.1)\nmesh = TreeMesh(0.0,40.0,initial_refinement_level=6,n_cells_max=10^4,periodicity=false)\nbc = (\n x_neg = boundary_condition_pressure_in,\n x_pos = Trixi.BoundaryConditionDoNothing()\n )\nvolume_flux = (flux_lax_friedrichs,flux_nonconservative)\nsurface_flux = (flux_lax_friedrichs,flux_nonconservative)\nbasis = LobattoLegendreBasis(2)\nid = IndicatorHennemannGassner(eq,basis;variable=first)\nvol = VolumeIntegralShockCapturingHG(id,volume_flux_dg=volume_flux,volume_flux_fv=surface_flux)\nsolver = DGSEM(basis,surface_flux,vol)\nsemi = SemidiscretizationHyperbolic(mesh,\neq,\ninitial_condition_simple,\nsource_terms = source_term_simple,\nsolver,\nboundary_conditions=bc)\nTrixi.default_analysis_integrals(::BloodFlowEquations1D) = ()\ntspan = (0.0, 0.5)\node = semidiscretize(semi, tspan)\nsummary_callback = SummaryCallback()\nanalysis_callback = AnalysisCallback(semi, interval = 200)\nstepsize_callback = StepsizeCallback(; cfl=0.5)\ncallbacks = CallbackSet(summary_callback,analysis_callback,stepsize_callback)\ndt = stepsize_callback(ode)\nsol = solve(ode, SSPRK33(), dt = dt, dtmax = 1e-4,dtmin = 1e-11,\n save_everystep = false,saveat = 0.002, callback = callbacks)\n\n@gif for i in eachindex(sol)\n a1 = sol[i][1:4:end]\n Q1 = sol[i][2:4:end]\n A01 = sol[i][4:4:end]\n A1 = A01.+a1\n plot(Q1./A1,lw=4,color=:red,ylim=(-10,50),label=\"velocity\",legend=:bottomleft)\nend","category":"page"},{"location":"tuto/","page":"Tutorial","title":"Tutorial","text":"(Image: Alt Text)","category":"page"},{"location":"tuto/#Tutorial-for-the-2D-model","page":"Tutorial","title":"Tutorial for the 2D model","text":"","category":"section"},{"location":"tuto/","page":"Tutorial","title":"Tutorial","text":"In this section, we describe how to use BloodFlowTrixi.jl with Trixi.jl. This tutorial will guide you through setting up and running a 2D blood flow simulation, including mesh creation, boundary conditions, numerical fluxes, and visualization of results.","category":"page"},{"location":"tuto/#Packages-2","page":"Tutorial","title":"Packages","text":"","category":"section"},{"location":"tuto/","page":"Tutorial","title":"Tutorial","text":"Before starting, ensure that the required packages are loaded:","category":"page"},{"location":"tuto/","page":"Tutorial","title":"Tutorial","text":"using Trixi\nusing BloodFlowTrixi\nusing OrdinaryDiffEq\nusing Plots","category":"page"},{"location":"tuto/","page":"Tutorial","title":"Tutorial","text":"First, we need to choose the equation that describes the blood flow dynamics:","category":"page"},{"location":"tuto/","page":"Tutorial","title":"Tutorial","text":"eq = BloodFlowEquations2D(; h=0.1)","category":"page"},{"location":"tuto/","page":"Tutorial","title":"Tutorial","text":"Here, h represents a parameter related to the initial condition or model scaling.","category":"page"},{"location":"tuto/#Mesh-and-boundary-conditions-2","page":"Tutorial","title":"Mesh and boundary conditions","text":"","category":"section"},{"location":"tuto/","page":"Tutorial","title":"Tutorial","text":"We begin by defining a two-dimensional P4est mesh, which discretizes the spatial domain:","category":"page"},{"location":"tuto/","page":"Tutorial","title":"Tutorial","text":"mesh = P4estMesh(\n (2,4),\n polydeg= 2,\n coordinates_min =(0.0,0.0),\n coordinates_max = (2*pi,40.0),\n initial_refinement_level = 4,\n periodicity = (true, false)\n)","category":"page"},{"location":"tuto/","page":"Tutorial","title":"Tutorial","text":"This generates a non-periodic mesh for the domain 02pi times 0 40, with 2times 4times 4^textinitial-refinement-level cells. ","category":"page"},{"location":"tuto/","page":"Tutorial","title":"Tutorial","text":"In Trixi.jl, the P4est mesh has four labeled boundaries: ***xneg*** (left boundary), ***xpos*** (right boundary), ***yneg*** (bottom boundary), and ***ypos*** (top boundary). These labels are used to apply boundary conditions:","category":"page"},{"location":"tuto/","page":"Tutorial","title":"Tutorial","text":"bc = Dict(\n :y_neg =>boundary_condition_pressure_in,\n :y_pos => Trixi.BoundaryConditionDoNothing()\n )\n","category":"page"},{"location":"tuto/","page":"Tutorial","title":"Tutorial","text":"boundary_condition_pressure_in applies a pressure inflow condition at the bottom boundary.\nTrixi.BoundaryConditionDoNothing() specifies a \"do nothing\" boundary condition at the right boundary, meaning no flux is imposed.","category":"page"},{"location":"tuto/#Boundary-condition-implementation-2","page":"Tutorial","title":"Boundary condition implementation","text":"","category":"section"},{"location":"tuto/","page":"Tutorial","title":"Tutorial","text":"The inflow boundary condition is defined as:","category":"page"},{"location":"tuto/","page":"Tutorial","title":"Tutorial","text":"boundary_condition_pressure_in(u_inner, normal, x, t, surface_flux_function, eq::BloodFlowEquations2D)","category":"page"},{"location":"tuto/","page":"Tutorial","title":"Tutorial","text":"This function applies a time-dependent pressure inflow condition.","category":"page"},{"location":"tuto/#Parameters-5","page":"Tutorial","title":"Parameters","text":"","category":"section"},{"location":"tuto/","page":"Tutorial","title":"Tutorial","text":"u_inner: State vector inside the domain near the boundary.\nnormal: Normal of the boundary.\nx: Position vector.\nt: Time scalar.\nsurface_flux_function: Function to compute flux at the boundary.\neq: Instance of BloodFlowEquations2D.","category":"page"},{"location":"tuto/#Returns-5","page":"Tutorial","title":"Returns","text":"","category":"section"},{"location":"tuto/","page":"Tutorial","title":"Tutorial","text":"The boundary flux is computed based on the inflow pressure: $ P{\\text{in}} = \\begin{cases} 2 \\times 10^4 \\sin^2\\left(\\frac{\\pi t}{0.125}\\right) & \\text{if } t < 0.125 \\\n0 & \\text{otherwise} \\end{cases} $ This time-dependent inflow pressure mimics a pulsatile flow, typical in arterial blood flow. The inflow area A{\\text{in}}$ is determined using the inverse pressure relation, ensuring consistency with the physical model.","category":"page"},{"location":"tuto/#Numerical-flux-2","page":"Tutorial","title":"Numerical flux","text":"","category":"section"},{"location":"tuto/","page":"Tutorial","title":"Tutorial","text":"To compute fluxes at cell interfaces, we use a combination of conservative and non-conservative fluxes:","category":"page"},{"location":"tuto/","page":"Tutorial","title":"Tutorial","text":"volume_flux = (flux_lax_friedrichs, flux_nonconservative)\nsurface_flux = (flux_lax_friedrichs, flux_nonconservative)","category":"page"},{"location":"tuto/","page":"Tutorial","title":"Tutorial","text":"flux_lax_friedrichs is a standard numerical flux for hyperbolic conservation laws.\nflux_nonconservative handles the non-conservative terms in the model, particularly those related to pressure discontinuities.","category":"page"},{"location":"tuto/","page":"Tutorial","title":"Tutorial","text":"The non-conservative flux function is defined as:","category":"page"},{"location":"tuto/","page":"Tutorial","title":"Tutorial","text":"flux_nonconservative(u_ll, u_rr, normal::Integer, eq::BloodFlowEquations2D)","category":"page"},{"location":"tuto/#Parameters-6","page":"Tutorial","title":"Parameters","text":"","category":"section"},{"location":"tuto/","page":"Tutorial","title":"Tutorial","text":"u_ll: Left state vector.\nu_rr: Right state vector.\nnormal: normal vector.\neq: Instance of BloodFlowEquations2D.","category":"page"},{"location":"tuto/#Returns-6","page":"Tutorial","title":"Returns","text":"","category":"section"},{"location":"tuto/","page":"Tutorial","title":"Tutorial","text":"The function returns the non-conservative flux vector, which is essential for capturing sharp pressure changes in the simulation.","category":"page"},{"location":"tuto/#Basis-functions-and-Shock-Capturing-DG-scheme-2","page":"Tutorial","title":"Basis functions and Shock Capturing DG scheme","text":"","category":"section"},{"location":"tuto/","page":"Tutorial","title":"Tutorial","text":"To approximate the solution, we use polynomial basis functions:","category":"page"},{"location":"tuto/","page":"Tutorial","title":"Tutorial","text":"basis = LobattoLegendreBasis(2)","category":"page"},{"location":"tuto/","page":"Tutorial","title":"Tutorial","text":"This defines a Lobatto-Legendre basis of polynomial degree 2, which is commonly used in high-order methods like Discontinuous Galerkin (DG) schemes.","category":"page"},{"location":"tuto/","page":"Tutorial","title":"Tutorial","text":"We then define an indicator for shock capturing, focusing on the first variable (area perturbation $a$):","category":"page"},{"location":"tuto/","page":"Tutorial","title":"Tutorial","text":"id = IndicatorHennemannGassner(eq, basis; variable=first)","category":"page"},{"location":"tuto/","page":"Tutorial","title":"Tutorial","text":"This indicator helps detect shocks or discontinuities in the solution and applies appropriate stabilization.","category":"page"},{"location":"tuto/","page":"Tutorial","title":"Tutorial","text":"The solver is defined as:","category":"page"},{"location":"tuto/","page":"Tutorial","title":"Tutorial","text":"vol = VolumeIntegralShockCapturingHG(id, volume_flux_dg=volume_flux, volume_flux_fv=surface_flux)\nsolver = DGSEM(basis, surface_flux, vol)","category":"page"},{"location":"tuto/","page":"Tutorial","title":"Tutorial","text":"Here, DGSEM represents the Discontinuous Galerkin Spectral Element Method, a high-order accurate scheme suitable for hyperbolic problems.","category":"page"},{"location":"tuto/#Semi-discretization-2","page":"Tutorial","title":"Semi-discretization","text":"","category":"section"},{"location":"tuto/","page":"Tutorial","title":"Tutorial","text":"We are now ready to semi-discretize the problem:","category":"page"},{"location":"tuto/","page":"Tutorial","title":"Tutorial","text":"semi = SemidiscretizationHyperbolic(\n mesh,\n eq,\n initial_condition_simple,\n source_terms = source_term_simple,\n solver,\n boundary_conditions=bc\n)","category":"page"},{"location":"tuto/","page":"Tutorial","title":"Tutorial","text":"This step sets up the semi-discretized form of the PDE, which will be advanced in time using an ODE solver.","category":"page"},{"location":"tuto/#Source-term-2","page":"Tutorial","title":"Source term","text":"","category":"section"},{"location":"tuto/","page":"Tutorial","title":"Tutorial","text":"The source term accounts for additional forces acting on the blood flow, such as friction:","category":"page"},{"location":"tuto/","page":"Tutorial","title":"Tutorial","text":"source_term_simple(u, x, t, eq::BloodFlowEquations2D)","category":"page"},{"location":"tuto/#Parameters-7","page":"Tutorial","title":"Parameters","text":"","category":"section"},{"location":"tuto/","page":"Tutorial","title":"Tutorial","text":"u: State vector containing area perturbation, flow rate, elasticity modulus, and reference area.\nx: Position vector.\nt: Time scalar.\neq::BloodFlowEquations2D: Instance of the blood flow model.","category":"page"},{"location":"tuto/#Returns-7","page":"Tutorial","title":"Returns","text":"","category":"section"},{"location":"tuto/","page":"Tutorial","title":"Tutorial","text":"The source term vector is given by:","category":"page"},{"location":"tuto/","page":"Tutorial","title":"Tutorial","text":"s_1 = 0\n(no source for area perturbation).\ns_2 = frac2R3mathcalC sin theta fracQ_s^2A + frac3Rk Q_RθA\n.\ns_3 = -frac2R3mathcalC sin theta fracQ_s Q_RthetaA + fracRkQ_sA\n.","category":"page"},{"location":"tuto/","page":"Tutorial","title":"Tutorial","text":"The friction coefficient k is computed using a model-specific friction function, and the radius R is obtained from the state vector using the radius function. Also, the curvature mathcalC is computed using the curvature function and is equal to 1 here.","category":"page"},{"location":"tuto/#Initial-condition-2","page":"Tutorial","title":"Initial condition","text":"","category":"section"},{"location":"tuto/","page":"Tutorial","title":"Tutorial","text":"The initial condition specifies the starting state of the simulation:","category":"page"},{"location":"tuto/","page":"Tutorial","title":"Tutorial","text":"initial_condition_simple(x, t, eq::BloodFlowEquations2D; R0=2.0)","category":"page"},{"location":"tuto/","page":"Tutorial","title":"Tutorial","text":"This function generates a simple initial condition with a uniform radius R0.","category":"page"},{"location":"tuto/#Parameters-8","page":"Tutorial","title":"Parameters","text":"","category":"section"},{"location":"tuto/","page":"Tutorial","title":"Tutorial","text":"x: Position vector.\nt: Time scalar.\neq::BloodFlowEquations2D: Instance of the blood flow model.\nR0: Initial radius (default: 2.0).","category":"page"},{"location":"tuto/#Returns-8","page":"Tutorial","title":"Returns","text":"","category":"section"},{"location":"tuto/","page":"Tutorial","title":"Tutorial","text":"The function returns a state vector with:","category":"page"},{"location":"tuto/","page":"Tutorial","title":"Tutorial","text":"Zero initial area perturbation.\nZero initial flow rate (in theta and s directions).\nConstant elasticity modulus.\nReference area A_0 = fracR_0^22.","category":"page"},{"location":"tuto/","page":"Tutorial","title":"Tutorial","text":"This simple initial condition is suitable for testing the model without introducing complex dynamics.","category":"page"},{"location":"tuto/#Run-the-simulation-2","page":"Tutorial","title":"Run the simulation","text":"","category":"section"},{"location":"tuto/","page":"Tutorial","title":"Tutorial","text":"First, we discretize the problem in time:","category":"page"},{"location":"tuto/","page":"Tutorial","title":"Tutorial","text":"Trixi.default_analysis_integrals(::BloodFlowEquations2D) = ()\ntspan = (0.0, 0.3)\node = semidiscretize(semi, tspan)","category":"page"},{"location":"tuto/","page":"Tutorial","title":"Tutorial","text":"Here, tspan defines the time interval for the simulation.","category":"page"},{"location":"tuto/","page":"Tutorial","title":"Tutorial","text":"Next, we add some callbacks to monitor the simulation:","category":"page"},{"location":"tuto/","page":"Tutorial","title":"Tutorial","text":"summary_callback = SummaryCallback()\nanalysis_callback = AnalysisCallback(semi, interval=200)\nstepsize_callback = StepsizeCallback(; cfl=0.5)\ncallbacks = CallbackSet(summary_callback, analysis_callback, stepsize_callback)","category":"page"},{"location":"tuto/","page":"Tutorial","title":"Tutorial","text":"SummaryCallback provides a summary of the simulation progress.\nAnalysisCallback computes analysis metrics at specified intervals.\nStepsizeCallback adjusts the time step based on the CFL condition.","category":"page"},{"location":"tuto/","page":"Tutorial","title":"Tutorial","text":"Finally, we solve the problem:","category":"page"},{"location":"tuto/","page":"Tutorial","title":"Tutorial","text":"dt = stepsize_callback(ode)\nsol = solve(ode, SSPRK33(), dt=dt, dtmax=1e-4, dtmin=1e-11,\n save_everystep=false, saveat=0.003, callback=callbacks)","category":"page"},{"location":"tuto/","page":"Tutorial","title":"Tutorial","text":"Here, SSPRK33() is a third-order Strong Stability Preserving Runge-Kutta method, suitable for hyperbolic PDEs.","category":"page"},{"location":"tuto/#Plot-the-results-2","page":"Tutorial","title":"Plot the results","text":"","category":"section"},{"location":"tuto/","page":"Tutorial","title":"Tutorial","text":"The results can be visualized using the following code:","category":"page"},{"location":"tuto/","page":"Tutorial","title":"Tutorial","text":"@gif for i in eachindex(sol)\npd = PlotData2D(sol[i],semi,solution_variables=cons2prim)\nplt1 = Plots.plot(pd[\"A\"],aspect_ratio=0.2)\nplt2 = Plots.plot(pd[\"wtheta\"],aspect_ratio=0.2)\nplt3 = Plots.plot(pd[\"ws\"],aspect_ratio=0.2)\nplt4 = Plots.plot(pd[\"P\"],aspect_ratio=0.2)\nplot(plt1,plt2,plt3,plt4,layout=(2,2))\nend","category":"page"},{"location":"tuto/","page":"Tutorial","title":"Tutorial","text":"This code generates an animated GIF showing the evolution of the velocity profile over time. The velocity is computed as QA, where Q is the flow rate, and A is the cross-sectional area.","category":"page"},{"location":"tuto/#Plain-code-2","page":"Tutorial","title":"Plain code","text":"","category":"section"},{"location":"tuto/","page":"Tutorial","title":"Tutorial","text":"using Trixi\nusing BloodFlowTrixi\nusing OrdinaryDiffEq,Plots\neq = BloodFlowEquations2D(;h=0.1)\nmesh = P4estMesh(\n (2,4),\n polydeg= 2,\n coordinates_min =(0.0,0.0),\n coordinates_max = (2*pi,40.0),\n initial_refinement_level = 4,\n periodicity = (true, false)\n)\nbc = Dict(\n :y_neg =>boundary_condition_pressure_in,\n :y_pos => Trixi.BoundaryConditionDoNothing()\n )\nvolume_flux = (flux_lax_friedrichs,flux_nonconservative)\nsurface_flux = (flux_lax_friedrichs,flux_nonconservative)\nbasis = LobattoLegendreBasis(2)\nid = IndicatorHennemannGassner(eq,basis;variable=first)\nvol =VolumeIntegralShockCapturingHG(id,volume_flux_dg = surface_flux,volume_flux_fv = volume_flux) \nsolver = DGSEM(basis,surface_flux,vol)\nsemi = SemidiscretizationHyperbolic(mesh,\neq,\ninitial_condition_simple,\nsource_terms = source_term_simple,\nsolver,\nboundary_conditions = bc)\nTrixi.default_analysis_integrals(::BloodFlowEquations2D) = ()\ntspan = (0.0, 0.3)\node = semidiscretize(semi, tspan)\nsummary_callback = SummaryCallback()\nanalysis_callback = AliveCallback(analysis_interval=1000)\nstepsize_callback = StepsizeCallback(; cfl=0.5)\ncallbacks = CallbackSet(summary_callback,analysis_callback,stepsize_callback)\ndt = stepsize_callback(ode)\nsol = solve(ode, SSPRK33(),dt=dt, dtmax = 1e-4,dtmin = 1e-12,save_everystep = false,saveat = 0.003, callback = callbacks)\n@gif for i in eachindex(sol)\npd = PlotData2D(sol[i],semi,solution_variables=cons2prim)\nplt1 = Plots.plot(pd[\"A\"],aspect_ratio=0.2)\nplt2 = Plots.plot(pd[\"wtheta\"],aspect_ratio=0.2)\nplt3 = Plots.plot(pd[\"ws\"],aspect_ratio=0.2)\nplt4 = Plots.plot(pd[\"P\"],aspect_ratio=0.2)\nplot(plt1,plt2,plt3,plt4,layout=(2,2))\nend","category":"page"},{"location":"tuto/","page":"Tutorial","title":"Tutorial","text":"(Image: Alt Text)","category":"page"},{"location":"math/#Modèles-mathématiques-1D-et-2D-pour-l’écoulement-sanguin","page":"Mathematics","title":"Modèles mathématiques 1D et 2D pour l’écoulement sanguin","text":"","category":"section"},{"location":"math/#Modèle-1D","page":"Mathematics","title":"Modèle 1D","text":"","category":"section"},{"location":"math/","page":"Mathematics","title":"Mathematics","text":"Le modèle 1D repose sur une intégration par section des équations de Navier-Stokes sous l’hypothèse d’un écoulement incompressible dans des artères supposées fines. Ce modèle est particulièrement adapté pour des études globales du réseau artériel, où la géométrie est approximativement linéaire ou faiblement courbée.","category":"page"},{"location":"math/#Hypothèses-et-simplifications","page":"Mathematics","title":"Hypothèses et simplifications","text":"","category":"section"},{"location":"math/","page":"Mathematics","title":"Mathematics","text":"L’écoulement est considéré comme incompressible.\nL’artère est modélisée comme un tube cylindrique de section variable en fonction de la pression.\nUn profil de vitesse parabolique est utilisé, permettant une moyennisation sur la section transversale de l’artère.","category":"page"},{"location":"math/#Équations-principales","page":"Mathematics","title":"Équations principales","text":"","category":"section"},{"location":"math/","page":"Mathematics","title":"Mathematics","text":"Les équations dérivées sont un système d’équations hyperboliques aux dérivées partielles décrivant la conservation de la masse et de la quantité de mouvement :","category":"page"},{"location":"math/","page":"Mathematics","title":"Mathematics","text":"Conservation de la masse :","category":"page"},{"location":"math/","page":"Mathematics","title":"Mathematics","text":"_t A + _x Q = 0","category":"page"},{"location":"math/","page":"Mathematics","title":"Mathematics","text":"Conservation de la quantité de mouvement :","category":"page"},{"location":"math/","page":"Mathematics","title":"Mathematics","text":"_t Q + _x left( alpha fracQ^2A + frac1rho A P(A x) right) - partial_x left( A partial_xleft(frac Q Aright) right) = frac1rho P(A x) _x A - K fracQA","category":"page"},{"location":"math/#Énergie-et-relation-d’entropie-du-modèle-1D","page":"Mathematics","title":"Énergie et relation d’entropie du modèle 1D","text":"","category":"section"},{"location":"math/","page":"Mathematics","title":"Mathematics","text":"L’énergie associée au système est donnée par :","category":"page"},{"location":"math/","page":"Mathematics","title":"Mathematics","text":"E(t x) = fracA u_x^22 + frac1rho A P(A x) - fracbeta(x)3 rho A_0(x) A^32","category":"page"},{"location":"math/","page":"Mathematics","title":"Mathematics","text":"La relation d’entropie vérifiée par cette énergie est :","category":"page"},{"location":"math/","page":"Mathematics","title":"Mathematics","text":"_t E + _x left( left( E + fracbeta(x)3 rho A_0(x) A^32 right) u_x right) = _x left( 3 nu A _x left( fracQA right) right) u_x + frac2 pi R k1 - R k 4 nu u_x^2 0","category":"page"},{"location":"math/","page":"Mathematics","title":"Mathematics","text":"Sous des conditions aux limites nulles :","category":"page"},{"location":"math/","page":"Mathematics","title":"Mathematics","text":"_t left( int_0^L E dx right) = - 3 nu int_0^L A (_x u_x)^2 dx - frac2 pi R k1 - R k 4 nu int_0^L u_x^2 dx 0","category":"page"},{"location":"math/","page":"Mathematics","title":"Mathematics","text":"","category":"page"},{"location":"math/#Modèle-2D","page":"Mathematics","title":"Modèle 2D","text":"","category":"section"},{"location":"math/","page":"Mathematics","title":"Mathematics","text":"Le modèle 2D est dérivé à partir d’une intégration radiale des équations de Navier-Stokes, permettant de mieux représenter les effets locaux dans des configurations géométriques complexes, comme les bifurcations artérielles et les anévrismes sévères.","category":"page"},{"location":"math/#Hypothèses-et-simplifications-2","page":"Mathematics","title":"Hypothèses et simplifications","text":"","category":"section"},{"location":"math/","page":"Mathematics","title":"Mathematics","text":"L’écoulement est supposé incompressible.\nLa géométrie de l’artère est décrite à l’aide d’un système de coordonnées curvilignes (( s, \theta )).\nLe profil de vitesse est obtenu sans recourir à un ansatz spécifique.","category":"page"},{"location":"math/#Équations-principales-2","page":"Mathematics","title":"Équations principales","text":"","category":"section"},{"location":"math/","page":"Mathematics","title":"Mathematics","text":"Conservation de la masse :","category":"page"},{"location":"math/","page":"Mathematics","title":"Mathematics","text":"_t A + _θ left( fracQ_RθA right) + _s(Q_s) = 0","category":"page"},{"location":"math/","page":"Mathematics","title":"Mathematics","text":"Conservation de la quantité de mouvement (composante radiale et axiale) :","category":"page"},{"location":"math/","page":"Mathematics","title":"Mathematics","text":"_t (Q_Rθ) + _θ left( fracQ_Rθ^22 A^2 + A P right) + _s left( fracQ_Rθ Q_sA right) = frac2 R3 C sin θ fracQ_s^2A + frac2 R k Q_RθA + _θ (A P)","category":"page"},{"location":"math/","page":"Mathematics","title":"Mathematics","text":"_t (Q_s) + _θ left( fracQ_s Q_RθA^2 right) + _s left( fracQ_s^2A - fracQ_Rθ^22 A^2 + A P right) = - frac2 R3 C sin θ fracQ_Rθ Q_sA^2 + frack R Q_sA + _s (A P)","category":"page"},{"location":"math/#Énergie-et-relation-d’entropie-du-modèle-2D","page":"Mathematics","title":"Énergie et relation d’entropie du modèle 2D","text":"","category":"section"},{"location":"math/","page":"Mathematics","title":"Mathematics","text":"L’énergie associée au système est donnée par :","category":"page"},{"location":"math/","page":"Mathematics","title":"Mathematics","text":"E(t θ s) = A left( frac98 u_θ^2 + fracu_s^22 + p right) - p","category":"page"},{"location":"math/","page":"Mathematics","title":"Mathematics","text":"La relation d’entropie correspondante est :","category":"page"},{"location":"math/","page":"Mathematics","title":"Mathematics","text":"_t E + _θ left( frac32 fracu_θR left( E + p - frac916 A u_θ^2 right) right) + _s left( u_s left( E + p - frac916 A u_θ^2 right) right) = frac94 R k u_θ^2 + k R u_s^2 0","category":"page"},{"location":"math/","page":"Mathematics","title":"Mathematics","text":"Cette relation garantit que l’énergie décroît localement dans le temps, ce qui assure la stabilité du modèle.","category":"page"},{"location":"math/","page":"Mathematics","title":"Mathematics","text":"","category":"page"},{"location":"math/#Comparaison-des-modèles-1D-et-2D","page":"Mathematics","title":"Comparaison des modèles 1D et 2D","text":"","category":"section"},{"location":"math/","page":"Mathematics","title":"Mathematics","text":"Modèle 1D :\nRapide et efficace pour des simulations globales sur de grands réseaux artériels.\nBien adapté pour des géométries simples ou faiblement courbées.\nCoût de calcul très faible.\nModèle 2D :\nPlus précis pour des géométries complexes (bifurcations, anévrismes).\nPermet de mieux capturer les effets locaux et les interactions fluide-structure.\nCoût de calcul modéré par rapport aux modèles tridimensionnels (NS-FSI 3D).","category":"page"},{"location":"math/","page":"Mathematics","title":"Mathematics","text":"L’utilisation combinée de ces deux modèles permet une alternative efficace aux simulations 3D, tout en offrant un bon compromis entre précision et coût de calcul.","category":"page"},{"location":"","page":"Home","title":"Home","text":"CurrentModule = BloodFlowTrixi","category":"page"},{"location":"#BloodFlowTrixi.jl","page":"Home","title":"BloodFlowTrixi.jl","text":"","category":"section"},{"location":"","page":"Home","title":"Home","text":"BloodFlowTrixi.jl is a Julia package that implements one-dimensional (1D) and two-dimensional (2D) blood flow models for arterial circulation. These models are derived from the Navier-Stokes equations and were developed as part of my PhD research in applied mathematics, focusing on cardiovascular pathologies such as aneurysms and stenoses.","category":"page"},{"location":"#Description","page":"Home","title":"Description","text":"","category":"section"},{"location":"","page":"Home","title":"Home","text":"This package provides:","category":"page"},{"location":"","page":"Home","title":"Home","text":"1D Blood Flow Model: This model describes blood flow along compliant arteries in a single spatial dimension. It was derived under the assumption of axisymmetric flow and accounts for arterial compliance, inertia, and frictional losses.\nMore details about this model can be found in my corresponding publication: [Article 1D]","category":"page"},{"location":"","page":"Home","title":"Home","text":"leftbeginaligned\n fracpartial apartial t + fracpartialpartial x(Q) = 0 \n fracpartial Qpartial t + fracpartialpartial xleft(fracQ^2A + A P(a)right) = P(a) fracpartial Apartial x - 2 pi R k frac Q A\n P(a) = P_ext + fracEhsqrtpi1-xi^2fracsqrtA - sqrtA_0A_0 \n R = sqrtfracApi\nendalignedright","category":"page"},{"location":"","page":"Home","title":"Home","text":"2D Blood Flow Model: The 2D model extends the Navier-Stokes equations under the thin-artery assumption, allowing for simulations in complex arterial geometries using curvilinear coordinates. It captures both longitudinal and angular dynamics, making it more accurate than classical 1D models while being less computationally expensive than full 3D models.\nThis model is described in detail in: [Article 2D]","category":"page"},{"location":"","page":"Home","title":"Home","text":"leftbeginaligned\n fracpartial apartial t + fracpartialpartial thetaleft( fracQ_RthetaA right) + fracpartialpartial s(Q_s) = 0 \n fracpartial Q_Rthetapartial t + fracpartialpartial thetaleft(fracQ_Rtheta^22A^2 + A P(a)right) + fracpartialpartial sleft( fracQ_RthetaQ_sA right) = P(a) fracpartial Apartial theta - 2 R k fracQ_RthetaA + frac2R3 mathcalCsin theta fracQ_s^2A \n fracpartial Q_spartial t + fracpartialpartial thetaleft(fracQ_Rtheta Q_sA^2 right) + fracpartialpartial sleft( fracQ_s^2A - fracQ_Rtheta^22A^2 + A P(a) right) = P(a) fracpartial Apartial s - R k fracQ_sA - frac2R3 mathcalCsin theta fracQ_s Q_RthetaA^2 \n P(a) = P_ext + fracEhsqrt2left(1-xi^2right)fracsqrtA - sqrtA_0A_0 \n R = sqrt2A\nendalignedright","category":"page"},{"location":"","page":"Home","title":"Home","text":"Both models were designed to be used with Trixi.jl, a flexible and high-performance framework for solving systems of conservation laws using the Discontinuous Galerkin (DG) method.","category":"page"},{"location":"#Features","page":"Home","title":"Features","text":"","category":"section"},{"location":"","page":"Home","title":"Home","text":"1D and 2D models for arterial blood flow.\nDerived from the Navier-Stokes equations with appropriate assumptions for compliant arteries.\nTo be used with Trixi.jl for DG-based numerical simulations.\nSupport for curvilinear geometries and compliant wall dynamics.","category":"page"},{"location":"#Installation","page":"Home","title":"Installation","text":"","category":"section"},{"location":"","page":"Home","title":"Home","text":"To install BloodFlowTrixi.jl, use the following commands in Julia:","category":"page"},{"location":"","page":"Home","title":"Home","text":"julia> ]\npkg> add Trixi\npkg> add https://github.com/your-repo/BloodFlowTrixi.jl","category":"page"},{"location":"#Future-Plans","page":"Home","title":"Future Plans","text":"","category":"section"},{"location":"","page":"Home","title":"Home","text":"short term","category":"page"},{"location":"","page":"Home","title":"Home","text":"Add second order 1D model.\nDesign prim variables for 1D and 2D models.\nAdd proper tests for 1D and 2D models.\nAdd 3D representations of the solutions for 1D and 2D models.\nDesign easy to use interfaces for users to define their own initial and boundary conditions and source terms.","category":"page"},{"location":"","page":"Home","title":"Home","text":"long term","category":"page"},{"location":"","page":"Home","title":"Home","text":"Add 3D fluid-structure interaction models for complex arterial geometries.\nDesign support for artery networks and simulate vascular networks using the 2D and 1D model.\nAutodiff support for 1D and 2D models for parameter optimization.","category":"page"},{"location":"#License","page":"Home","title":"License","text":"","category":"section"},{"location":"","page":"Home","title":"Home","text":"This package is licensed under the MIT license.","category":"page"},{"location":"#Acknowledgments","page":"Home","title":"Acknowledgments","text":"","category":"section"},{"location":"","page":"Home","title":"Home","text":"This package was developed as part of my PhD research in applied mathematics, focusing on mathematical modeling and numerical simulation of blood flow in arteries. Special thanks to the developers of Trixi.jl, whose framework was invaluable in implementing and testing these models.","category":"page"},{"location":"","page":"Home","title":"Home","text":"","category":"page"},{"location":"","page":"Home","title":"Home","text":"Modules = [BloodFlowTrixi]","category":"page"},{"location":"#BloodFlowTrixi.BloodFlowTrixi","page":"Home","title":"BloodFlowTrixi.BloodFlowTrixi","text":"Package BloodFlowTrixi v0.0.1\n\nThis package implements 1D and 2D blood flow models for arterial circulation using Trixi.jl, enabling efficient numerical simulation and analysis.\n\nDocs under https://yolhan83.github.io/BloodFlowTrixi.jl\n\n\n\n\n\n","category":"module"},{"location":"#BloodFlowTrixi.BloodFlowEquations1D","page":"Home","title":"BloodFlowTrixi.BloodFlowEquations1D","text":"BloodFlowEquations1D(;h,rho=1.0,xi=0.25,nu=0.04)\n\nBlood Flow equations in one space dimension. This model describes the dynamics of blood flow along a compliant artery using one-dimensional equations derived from the Navier-Stokes equations. The equations account for conservation of mass and momentum, incorporating the effect of arterial compliance and frictional losses.\n\nThe governing equations are given by\n\nleftbeginaligned\n fracpartial apartial t + fracpartialpartial x(Q) = 0 \n fracpartial Qpartial t + fracpartialpartial xleft(fracQ^2A + A P(a)right) = P(a) fracpartial Apartial x - 2 pi R k frac Q A\n P(a) = P_ext + fracEhsqrtpi1-xi^2fracsqrtA - sqrtA_0A_0 \n R = sqrtfracApi\nendalignedright\n\n\n\n\n\n","category":"type"},{"location":"#BloodFlowTrixi.BloodFlowEquations2D","page":"Home","title":"BloodFlowTrixi.BloodFlowEquations2D","text":"BloodFlowEquations2D(;h,rho=1.0,xi=0.25)\n\nDefines the two-dimensional blood flow equations derived from the Navier-Stokes equations in curvilinear coordinates under the thin-artery assumption. This model describes the dynamics of blood flow along a compliant artery in two spatial dimensions (s, θ).\n\nParameters\n\nh::T: Wall thickness of the artery.\nrho::T: Fluid density (default 1.0).\nxi::T: Poisson's ratio (default 0.25).\nnu::T: Viscosity coefficient.\n\nThe governing equations account for conservation of mass and momentum, incorporating the effects of arterial compliance, curvature, and frictional losses.\n\nleftbeginaligned\n fracpartial apartial t + fracpartialpartial thetaleft( fracQ_RthetaA right) + fracpartialpartial s(Q_s) = 0 \n fracpartial Q_Rthetapartial t + fracpartialpartial thetaleft(fracQ_Rtheta^22A^2 + A P(a)right) + fracpartialpartial sleft( fracQ_RthetaQ_sA right) = P(a) fracpartial Apartial theta - 2 R k fracQ_RthetaA + frac2R3 mathcalCsin theta fracQ_s^2A \n fracpartial Q_spartial t + fracpartialpartial thetaleft(fracQ_Rtheta Q_sA^2 right) + fracpartialpartial sleft( fracQ_s^2A - fracQ_Rtheta^22A^2 + A P(a) right) = P(a) fracpartial Apartial s - R k fracQ_sA - frac2R3 mathcalCsin theta fracQ_s Q_RthetaA^2 \n P(a) = P_ext + fracEhsqrt2left(1-xi^2right)fracsqrtA - sqrtA_0A_0 \n R = sqrt2A\nendalignedright\n\n\n\n\n\n","category":"type"},{"location":"#Trixi.DissipationLocalLaxFriedrichs-Tuple{Any, Any, Any, BloodFlowEquations1D}","page":"Home","title":"Trixi.DissipationLocalLaxFriedrichs","text":"(dissipation::Trixi.DissipationLocalLaxFriedrichs)(u_ll, u_rr, orientation_or_normal_direction, eq::BloodFlowEquations1D)\n\nCalculates the dissipation term using the Local Lax-Friedrichs method.\n\nParameters\n\nu_ll: Left state vector.\nu_rr: Right state vector.\norientation_or_normal_direction: Orientation or normal direction.\neq: Instance of BloodFlowEquations1D.\n\nReturns\n\nDissipation vector.\n\n\n\n\n\n","category":"method"},{"location":"#BloodFlowTrixi.boundary_condition_outflow-Tuple{Any, Any, Any, Any, Any, Any, BloodFlowEquations1D}","page":"Home","title":"BloodFlowTrixi.boundary_condition_outflow","text":"boundary_condition_outflow(u_inner, orientation_or_normal, direction, x, t, surface_flux_function, eq::BloodFlowEquations1D)\n\nImplements the outflow boundary condition, assuming that there is no reflection at the boundary.\n\nParameters\n\nu_inner: State vector inside the domain near the boundary.\norientation_or_normal: Normal orientation of the boundary.\ndirection: Integer indicating the direction of the boundary.\nx: Position vector.\nt: Time.\nsurface_flux_function: Function to compute flux at the boundary.\neq: Instance of BloodFlowEquations1D.\n\nReturns\n\nComputed boundary flux.\n\n\n\n\n\n","category":"method"},{"location":"#BloodFlowTrixi.boundary_condition_pressure_in-Tuple{Any, Any, Any, Any, Any, Any, BloodFlowEquations1D}","page":"Home","title":"BloodFlowTrixi.boundary_condition_pressure_in","text":"boundary_condition_pressure_in(u_inner, orientation_or_normal, direction, x, t, surface_flux_function, eq::BloodFlowEquations1D)\n\nImplements a pressure inflow boundary condition where the inflow pressure varies with time.\n\nParameters\n\nu_inner: State vector inside the domain near the boundary.\norientation_or_normal: Normal orientation of the boundary.\ndirection: Integer indicating the boundary direction.\nx: Position vector.\nt: Time scalar.\nsurface_flux_function: Function to compute flux at the boundary.\neq: Instance of BloodFlowEquations1D.\n\nReturns\n\nComputed boundary flux with inflow pressure specified by:\n\nP_in = begincases\n2 times 10^4 sin^2(pi t 0125) textif t 0125 \n0 textotherwise\nendcases\n\nThe corresponding inflow area A_{in} is computed using the inverse pressure relation, and the boundary state is constructed accordingly.\n\n\n\n\n\n","category":"method"},{"location":"#BloodFlowTrixi.boundary_condition_slip_wall-Tuple{Any, Any, Any, Any, Any, Any, BloodFlowEquations1D}","page":"Home","title":"BloodFlowTrixi.boundary_condition_slip_wall","text":"boundary_condition_slip_wall(u_inner, orientation_or_normal, direction, x, t, surface_flux_function, eq::BloodFlowEquations1D)\n\nImplements a slip wall boundary condition where the normal component of velocity is reflected.\n\nParameters\n\nu_inner: State vector inside the domain near the boundary.\norientation_or_normal: Normal orientation of the boundary.\ndirection: Integer indicating the direction of the boundary.\nx: Position vector.\nt: Time.\nsurface_flux_function: Function to compute flux at the boundary.\neq: Instance of BloodFlowEquations1D.\n\nReturns\n\nComputed boundary flux at the slip wall.\n\n\n\n\n\n","category":"method"},{"location":"#BloodFlowTrixi.flux_nonconservative-Tuple{Any, Any, Integer, BloodFlowEquations1D}","page":"Home","title":"BloodFlowTrixi.flux_nonconservative","text":"flux_nonconservative(u_ll, u_rr, orientation::Integer, eq::BloodFlowEquations1D)\n\nComputes the non-conservative flux for the model, used for handling discontinuities in pressure.\n\nParameters\n\nu_ll: Left state vector.\nu_rr: Right state vector.\norientation::Integer: Orientation index.\neq: Instance of BloodFlowEquations1D.\n\nReturns\n\nNon-conservative flux vector.\n\n\n\n\n\n","category":"method"},{"location":"#BloodFlowTrixi.friction-Tuple{Any, Any, BloodFlowEquations1D}","page":"Home","title":"BloodFlowTrixi.friction","text":"friction(u, x, eq::BloodFlowEquations1D)\n\nCalculates the friction term for the blood flow equations, which represents viscous resistance to flow along the artery wall.\n\nParameters\n\nu: State vector containing cross-sectional area and flow rate.\nx: Position along the artery.\neq::BloodFlowEquations1D: Instance of the blood flow model.\n\nReturns\n\nFriction coefficient as a scalar.\n\n\n\n\n\n","category":"method"},{"location":"#BloodFlowTrixi.initial_condition_simple-Tuple{Any, Any, BloodFlowEquations1D}","page":"Home","title":"BloodFlowTrixi.initial_condition_simple","text":"initial_condition_simple(x, t, eq::BloodFlowEquations1D; R0=2.0)\n\nGenerates a simple initial condition with a specified initial radius R0.\n\nParameters\n\nx: Position vector.\nt: Time scalar.\neq::BloodFlowEquations1D: Instance of the blood flow model.\nR0: Initial radius (default: 2.0).\n\nReturns\n\nState vector with zero initial area perturbation, zero flow rate, constant elasticity modulus, and reference area computed as A_0 = \\pi R_0^2.\n\nThis initial condition is suitable for basic tests without complex dynamics.\n\n\n\n\n\n","category":"method"},{"location":"#BloodFlowTrixi.inv_pressure-Tuple{Any, Any, BloodFlowEquations1D}","page":"Home","title":"BloodFlowTrixi.inv_pressure","text":"inv_pressure(p, u, eq::BloodFlowEquations1D)\n\nComputes the inverse relation of pressure to cross-sectional area.\n\nParameters\n\np: Pressure.\nu: State vector.\neq: Instance of BloodFlowEquations1D.\n\nReturns\n\nCross-sectional area corresponding to the given pressure.\n\n\n\n\n\n","category":"method"},{"location":"#BloodFlowTrixi.pressure-Tuple{Any, BloodFlowEquations1D}","page":"Home","title":"BloodFlowTrixi.pressure","text":"pressure(u, eq::BloodFlowEquations1D)\n\nComputes the pressure given the state vector based on the compliance of the artery.\n\nParameters\n\nu: State vector.\neq: Instance of BloodFlowEquations1D.\n\nReturns\n\nPressure as a scalar.\n\n\n\n\n\n","category":"method"},{"location":"#BloodFlowTrixi.pressure_der-Tuple{Any, BloodFlowEquations1D}","page":"Home","title":"BloodFlowTrixi.pressure_der","text":"pressure_der(u, eq::BloodFlowEquations1D)\n\nComputes the derivative of pressure with respect to cross-sectional area.\n\nParameters\n\nu: State vector.\neq: Instance of BloodFlowEquations1D.\n\nReturns\n\nDerivative of pressure.\n\n\n\n\n\n","category":"method"},{"location":"#BloodFlowTrixi.radius-Tuple{Any, BloodFlowEquations1D}","page":"Home","title":"BloodFlowTrixi.radius","text":"radius(u, eq::BloodFlowEquations1D)\n\nComputes the radius of the artery based on the cross-sectional area.\n\nParameters\n\nu: State vector.\neq: Instance of BloodFlowEquations1D.\n\nReturns\n\nRadius as a scalar.\n\n\n\n\n\n","category":"method"},{"location":"#BloodFlowTrixi.source_term_simple-Tuple{Any, Any, Any, BloodFlowEquations1D}","page":"Home","title":"BloodFlowTrixi.source_term_simple","text":"source_term_simple(u, x, t, eq::BloodFlowEquations1D)\n\nComputes a simple source term for the blood flow model, focusing on frictional effects.\n\nParameters\n\nu: State vector containing area perturbation, flow rate, elasticity modulus, and reference area.\nx: Position vector.\nt: Time scalar.\neq::BloodFlowEquations1D: Instance of the blood flow model.\n\nReturns\n\nSource terms vector where:\n\ns_1 = 0 (no source for area perturbation).\ns_2 represents the friction term given by s_2 = \\frac{2 \\pi k Q}{R A}.\n\nFriction coefficient k is computed using the friction function, and the radius R is obtained using the radius function.\n\n\n\n\n\n","category":"method"},{"location":"#Trixi.cons2prim-Tuple{Any, BloodFlowEquations1D}","page":"Home","title":"Trixi.cons2prim","text":"Trixi.cons2prim(u, eq::BloodFlowEquations1D)\n\nConverts the conserved variables to primitive variables.\n\nParameters\n\nu: State vector.\neq: Instance of BloodFlowEquations1D.\n\nReturns\n\nPrimitive variable vector.\n\n\n\n\n\n","category":"method"},{"location":"#Trixi.flux-Tuple{Any, Integer, BloodFlowEquations1D}","page":"Home","title":"Trixi.flux","text":"Trixi.flux(u, orientation::Integer, eq::BloodFlowEquations1D)\n\nComputes the flux vector for the conservation laws of the blood flow model.\n\nParameters\n\nu: State vector.\norientation::Integer: Orientation index for flux computation.\neq: Instance of BloodFlowEquations1D.\n\nReturns\n\nFlux vector as an SVector.\n\n\n\n\n\n","category":"method"},{"location":"#Trixi.initial_condition_convergence_test-Tuple{Any, Any, BloodFlowEquations1D}","page":"Home","title":"Trixi.initial_condition_convergence_test","text":"initial_condition_convergence_test(x, t, eq::BloodFlowEquations1D)\n\nGenerates a smooth initial condition for convergence tests of the blood flow equations.\n\nParameters\n\nx: Position vector.\nt: Time scalar.\neq::BloodFlowEquations1D: Instance of the blood flow model.\n\nReturns\n\nInitial condition state vector with zero initial area perturbation, sinusoidal flow rate, a constant elasticity modulus, and reference area.\n\nDetails\n\nThe returned initial condition has:\n\nZero perturbation in area (a = 0).\nA sinusoidal flow rate given by Q = sin(\\pi x t).\nA constant elasticity modulus E.\nA reference cross-sectional area A_0 = \\pi R_0^2 for R_0 = 1.\n\nThis initial condition can be used to verify the accuracy and stability of numerical solvers.\n\n\n\n\n\n","category":"method"},{"location":"#Trixi.max_abs_speed_naive-Tuple{Any, Any, Integer, BloodFlowEquations1D}","page":"Home","title":"Trixi.max_abs_speed_naive","text":"Trixi.max_abs_speed_naive(u_ll, u_rr, orientation::Integer, eq::BloodFlowEquations1D)\n\nCalculates the maximum absolute speed for wave propagation in the blood flow model using a naive approach.\n\nParameters\n\nu_ll: Left state vector.\nu_rr: Right state vector.\norientation::Integer: Orientation index.\neq: Instance of BloodFlowEquations1D.\n\nReturns\n\nMaximum absolute speed.\n\n\n\n\n\n","category":"method"},{"location":"#Trixi.prim2cons-Tuple{Any, BloodFlowEquations1D}","page":"Home","title":"Trixi.prim2cons","text":"Trixi.prim2cons(u, eq::BloodFlowEquations1D)\n\nConverts the primitive variables to conserved variables.\n\nParameters\n\nu: Primitive variable vector.\neq: Instance of BloodFlowEquations1D.\n\nReturns\n\nConserved variable vector.\n\n\n\n\n\n","category":"method"},{"location":"#Trixi.source_terms_convergence_test-Tuple{Any, Any, Any, BloodFlowEquations1D}","page":"Home","title":"Trixi.source_terms_convergence_test","text":"source_terms_convergence_test(u, x, t, eq::BloodFlowEquations1D)\n\nComputes the source terms for convergence tests of the blood flow equations.\n\nParameters\n\nu: State vector containing area perturbation, flow rate, elasticity modulus, and reference area.\nx: Position vector.\nt: Time scalar.\neq::BloodFlowEquations1D: Instance of the blood flow model.\n\nReturns\n\nSource terms vector.\n\nDetails\n\nThe source terms are derived based on the smooth initial condition and friction effects:\n\ns_1 represents the source term for area perturbation and is given by s_1 = \\pi t \\cos(\\pi x t).\ns_2 represents the source term for the flow rate and includes contributions from spatial and temporal variations as well as friction effects.\n\nThe radius R is computed using the radius function, and the friction coefficient k is obtained using the friction function.\n\nThis function is useful for evaluating the correctness of source term handling in numerical solvers.\n\n\n\n\n\n","category":"method"}] } diff --git a/dev/tuto/index.html b/dev/tuto/index.html index 0603f9d..843fd9c 100644 --- a/dev/tuto/index.html +++ b/dev/tuto/index.html @@ -148,4 +148,4 @@ plt3 = Plots.plot(pd["ws"],aspect_ratio=0.2) plt4 = Plots.plot(pd["P"],aspect_ratio=0.2) plot(plt1,plt2,plt3,plt4,layout=(2,2)) -end

Alt Text

+end

Alt Text