Add second order 1D model

README.md

@@ -61,7 +61,6 @@ pkg> add https://github.com/your-repo/BloodFlowTrixi.jl
 
 **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.

docs/src/math.md

@@ -56,10 +56,10 @@ Le modèle 2D est dérivé à partir d’une **intégration radiale des équatio
 
 2. **Conservation de la quantité de mouvement (composante radiale et axiale)** :
 ```math
-∂_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_{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} + P∂_θ (A)
 ```
 ```math
-∂_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)
+∂_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} + P∂_s (A)
 ```
 
 ### Énergie et relation d’entropie du modèle 2D

src/1DModel/1dmodel.jl

@@ -148,4 +148,5 @@ end
 include("./variables.jl")
 include("./bc1d.jl")
 include("./Test_Cases/pressure_in.jl")
-include("./Test_Cases/convergence_test.jl")
+include("./Test_Cases/convergence_test.jl")
+include("Ord2/1dmodelord2.jl")

src/1DModel/Ord2/1dmodelord2.jl

@@ -0,0 +1,68 @@
+struct BloodFlowEquations1DOrd2{T <:Real,E} <: Trixi.AbstractEquationsParabolic{1, 4, GradientVariablesConservative}
+    nu ::T
+    model1d ::E
+end
+Trixi.varnames(mapin,eq::BloodFlowTrixi.BloodFlowEquations1DOrd2) = Trixi.varnames(mapin,eq.model1d)
+
+function Trixi.flux(u,gradients,orientation::Int,eq_parab ::BloodFlowEquations1DOrd2)
+    dudx = gradients
+    a,Q,_,A0 = u
+    A = a+A0
+    val = 3*eq_parab.nu * (-(dudx[1] + dudx[4])*Q/A + dudx[2])
+    return SVector(0.0,val,0,0)
+end
+
+function source_term_simple(u, x, t, eq::BloodFlowEquations1DOrd2)
+    res = source_term_simple(u, x, t, eq.model1d)
+    k = friction(u,x,eq.model1d)
+    R = radius(u,eq.model1d)
+    return SVector(res[1],res[2]/(1-R*k/(4*eq.nu)),res[3],res[4])
+end
+
+@inline function boundary_condition_pressure_in(flux_inner, u_inner,
+    orientation_or_normal,direction,
+    x, t,
+    operator_type::Trixi.Gradient,
+    equations_parabolic::BloodFlowEquations1DOrd2)
+return boundary_condition_pressure_in(u_inner,orientation_or_normal,direction,x,t,flux_lax_friedrichs,equations_parabolic.model1d)
+end
+@inline function boundary_condition_pressure_in(flux_inner, u_inner,
+    orientation_or_normal,direction,
+    x, t,
+    operator_type::Trixi.Divergence,
+    equations_parabolic::BloodFlowEquations1DOrd2)
+return flux_inner
+end
+# Dirichlet and Neumann boundary conditions for use with parabolic solvers in weak form.
+# Note that these are general, so they apply to LaplaceDiffusion in any spatial dimension.
+@inline function (boundary_condition::Trixi.BoundaryConditionDirichlet)(flux_inner, u_inner,
+    normal::AbstractVector,
+    x, t,
+    operator_type::Trixi.Gradient,
+    equations_parabolic::BloodFlowEquations1DOrd2)
+return boundary_condition.boundary_value_function(x, t, equations_parabolic)
+end
+
+@inline function (boundary_condition::Trixi.BoundaryConditionDirichlet)(flux_inner, u_inner,
+    normal::AbstractVector,
+    x, t,
+    operator_type::Trixi.Divergence,
+    equations_parabolic::BloodFlowEquations1DOrd2)
+return flux_inner
+end
+
+@inline function (boundary_condition::Trixi.BoundaryConditionNeumann)(flux_inner, u_inner,
+  normal::AbstractVector,
+  x, t,
+  operator_type::Trixi.Divergence,
+  equations_parabolic::BloodFlowEquations1DOrd2)
+return boundary_condition.boundary_normal_flux_function(x, t, equations_parabolic)
+end
+
+@inline function (boundary_condition::Trixi.BoundaryConditionNeumann)(flux_inner, u_inner,
+  normal::AbstractVector,
+  x, t,
+  operator_type::Trixi.Gradient,
+  equations_parabolic::BloodFlowEquations1DOrd2)
+return flux_inner
+end