diff --git a/.JuliaFormatter.toml b/.JuliaFormatter.toml
new file mode 100644
index 0000000..323237b
--- /dev/null
+++ b/.JuliaFormatter.toml
@@ -0,0 +1 @@
+style = "blue"
diff --git a/.github/workflows/CI.yml b/.github/workflows/CI.yml
index deb0e65..d9c5f91 100644
--- a/.github/workflows/CI.yml
+++ b/.github/workflows/CI.yml
@@ -39,6 +39,13 @@ jobs:
       - uses: julia-actions/cache@v2
       - uses: julia-actions/julia-buildpkg@v1
       - uses: julia-actions/julia-runtest@v1
+      - uses: julia-actions/julia-processcoverage@v1
+        if: matrix.version == '1.10' && matrix.os == 'ubuntu-latest'
+      - uses: codecov/codecov-action@v5
+        if: matrix.version == '1.10' && matrix.os == 'ubuntu-latest'
+        with:
+          token: ${{ secrets.CODECOV_TOKEN }}
+          fail_ci_if_error: false
   docs:
     name: Documentation
     runs-on: ubuntu-latest
diff --git a/CODE_OF_CONDUCT.md b/CODE_OF_CONDUCT.md
new file mode 100644
index 0000000..06729bd
--- /dev/null
+++ b/CODE_OF_CONDUCT.md
@@ -0,0 +1,6 @@
+JuliaHealth Community Standards
+=========================
+
+The JuliaHealth community abides by the broader Julia Community Code of Conduct (CoC). You can find that CoC listed here: https://julialang.org/community/standards/
+
+If you have a conflict or concern that requires resolution, please contact the [Julia Community Stewards](https://julialang.org/community/stewards/).
diff --git a/CONTRIBUTING.md b/CONTRIBUTING.md
new file mode 100644
index 0000000..72407fb
--- /dev/null
+++ b/CONTRIBUTING.md
@@ -0,0 +1,41 @@
+# Contributing Guidelines
+
+Thanks for taking the time to contribute to BloodFlowTrixi.jl! ❤️
+
+We appreciate your interest in improving this project. Before you start contributing, please take a moment to review the following guidelines.
+
+
+## Reporting Issues
+
+First off, we assume that you have read the available [Documentation](https://juliahealth.org/BloodFlowTrixi.jl/dev/).
+
+Before you report an issue, it is best to search for existing [Issues](https://github.com/JuliaHealth/BloodFlowTrixi.jl/issues) that might help you. In case you have found a suitable issue and still need clarification, you can write your question in this issue.
+
+If you then still feel the need to report an issue, we recommend the following:
+
+- Open an [Issue on GitHub](https://github.com/JuliaHealth/BloodFlowTrixi.jl/issues/new).
+- Provide as much context as you can about what you're running into.
+- Provide project and platform versions, depending on what seems relevant.
+
+We will then take care of the issue as soon as possible.
+
+
+## How to Contribute
+
+- **Implement Code Changes**: Introduce your modifications, ensuring adherence to the [Julia Blue Style Guidelines](https://github.com/invenia/BlueStyle). For new features, include informative comments, docstrings, and consider enriching the documentation with relevant examples.
+
+- **Format Code**: This project uses JuliaFormatter with Blue style. Run formatting from the repository root using these steps:
+	- `julia --project=dev`
+	- `using JuliaFormatter`
+	- `format(".")`
+
+- **Validate Changes with Tests**: Execute existing tests to verify the compatibility of your alterations with the current functionality. If applicable, incorporate additional tests to validate your new contributions.
+
+- **Undergo Code Review**: Subject your code to review by maintainers, who will provide feedback and may request further adjustments before merging.
+
+
+## License
+
+By contributing to BloodFlowTrixi.jl, you agree that your contributions will be licensed under the [MIT License](LICENSE).
+
+Thank you for contributing to BloodFlowTrixi.jl! 🌟
diff --git a/Project.toml b/Project.toml
index cddd7d6..6c5d242 100644
--- a/Project.toml
+++ b/Project.toml
@@ -1,7 +1,7 @@
 name = "BloodFlowTrixi"
 uuid = "7b46d5ad-635f-4b66-9393-904f183e39b5"
-version = "0.1.7"
 authors = ["yolhan83 <yolhan@laposte.net>"]
+version = "0.1.7"
 
 [deps]
 ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
diff --git a/README.md b/README.md
index 691ec58..d3c71da 100644
--- a/README.md
+++ b/README.md
@@ -3,6 +3,7 @@
 [![Stable](https://img.shields.io/badge/docs-stable-blue.svg)](https://juliahealth.org/BloodFlowTrixi.jl/stable/)
 [![Dev](https://img.shields.io/badge/docs-dev-blue.svg)](https://juliahealth.org/BloodFlowTrixi.jl/dev/)
 [![Build Status](https://github.com/yolhan83/BloodFlowTrixi.jl/actions/workflows/CI.yml/badge.svg?branch=master)](https://github.com/yolhan83/BloodFlowTrixi.jl/actions/workflows/CI.yml?query=branch%3Amaster)
+[![Coverage](https://codecov.io/gh/JuliaHealth/BloodFlowTrixi.jl/branch/master/graph/badge.svg)](https://codecov.io/gh/JuliaHealth/BloodFlowTrixi.jl)
 [![Aqua](https://raw.githubusercontent.com/JuliaTesting/Aqua.jl/master/badge.svg)](https://github.com/JuliaTesting/Aqua.jl)
 
 **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.
@@ -57,6 +58,10 @@ pkg> add Trixi
 pkg> add BloodFlowTrixi
 ```
 
+## Usage Examples
+
+For detailed usage examples and tutorials on how to use BloodFlowTrixi.jl, please visit the [Tutorial Documentation](https://juliahealth.org/BloodFlowTrixi.jl/stable/tuto/).
+
 # Simulations
 ![Alt Text](./docs/src/graph.gif)
 ![Alt Text](./docs/src/graph2d.gif)
@@ -80,6 +85,14 @@ pkg> add BloodFlowTrixi
 
 This package is licensed under the MIT license.
 
+## Contributing
+
+We welcome contributions from the community! If you find a bug, have a feature request, or want to improve the package, please:
+
+1. Check the [Contributing Guidelines](CONTRIBUTING.md)
+2. See the [Code of Conduct](CODE_OF_CONDUCT.md)
+3. Open a [new issue](https://github.com/JuliaHealth/BloodFlowTrixi.jl/issues) or [pull request](https://github.com/JuliaHealth/BloodFlowTrixi.jl/pulls)
+
 ## Acknowledgments
 
 This package was developed as part of my PhD research in applied mathematics in the IMATH laboratory and founded by the University of Toulon in France, focusing on mathematical modeling and numerical simulation of blood flow in arteries. 
diff --git a/dev/.gitignore b/dev/.gitignore
new file mode 100644
index 0000000..f3ca75e
--- /dev/null
+++ b/dev/.gitignore
@@ -0,0 +1 @@
+/Manifest.toml
\ No newline at end of file
diff --git a/dev/Project.toml b/dev/Project.toml
new file mode 100644
index 0000000..7fbe771
--- /dev/null
+++ b/dev/Project.toml
@@ -0,0 +1,5 @@
+[deps]
+JuliaFormatter = "98e50ef6-434e-11e9-1051-2b60c6c9e899"
+
+[compat]
+JuliaFormatter = "2.3.0"
diff --git a/docs/make.jl b/docs/make.jl
index 1d7e087..7d1e53f 100644
--- a/docs/make.jl
+++ b/docs/make.jl
@@ -12,14 +12,7 @@ makedocs(;
         edit_link="master",
         assets=String[],
     ),
-    pages=[
-        "Home" => "index.md",
-        "Tutorial" => "tuto.md",
-        "Mathematics" => "math.md"
-    ],
+    pages=["Home" => "index.md", "Tutorial" => "tuto.md", "Mathematics" => "math.md"],
 )
 
-deploydocs(;
-    repo="github.com/JuliaHealth/BloodFlowTrixi.jl",
-    devbranch="master",
-)
+deploydocs(; repo="github.com/JuliaHealth/BloodFlowTrixi.jl", devbranch="master")
diff --git a/exemples/Model1D/exemple.jl b/exemples/Model1D/exemple.jl
index 7af0752..53df303 100644
--- a/exemples/Model1D/exemple.jl
+++ b/exemples/Model1D/exemple.jl
@@ -2,42 +2,34 @@ using Trixi
 using BloodFlowTrixi
 using OrdinaryDiffEq
 
-eq = BloodFlowEquations1D(; h = 0.1)
+eq = BloodFlowEquations1D(; h=0.1)
 
-mesh = TreeMesh(0.0, 40.0, 
-    initial_refinement_level = 4, 
-    n_cells_max = 10^4, 
-    periodicity = false)
+mesh = TreeMesh(0.0, 40.0; initial_refinement_level=4, n_cells_max=10^4, periodicity=false)
 
-bc = (;
-    x_neg = boundary_condition_pressure_in,
-    x_pos = Trixi.BoundaryConditionDoNothing()
-)
+bc = (; x_neg=boundary_condition_pressure_in, x_pos=Trixi.BoundaryConditionDoNothing())
 
-solver = DGSEM(polydeg = 2,
-    surface_flux = (flux_lax_friedrichs,flux_nonconservative),
-    volume_integral = VolumeIntegralFluxDifferencing((flux_lax_friedrichs,flux_nonconservative))
-    )
+solver = DGSEM(;
+    polydeg=2,
+    surface_flux=(flux_lax_friedrichs, flux_nonconservative),
+    volume_integral=VolumeIntegralFluxDifferencing((
+        flux_lax_friedrichs, flux_nonconservative
+    )),
+)
 
 semi = SemidiscretizationHyperbolic(
     mesh,
     eq,
     initial_condition_simple,
-    source_terms = source_term_simple,
-    solver,
-    boundary_conditions = bc
+    solver;
+    source_terms=source_term_simple,
+    boundary_conditions=bc,
 )
 
 tspan = (0.0, 0.3)
 ode = semidiscretize(semi, tspan)
 
-dt_adapt = StepsizeCallback(;cfl=0.5)
-analyse = AliveCallback(
-    alive_interval = 10,
-    analysis_interval = 100
-)
-cb = CallbackSet(
-    dt_adapt,analyse
-)
+dt_adapt = StepsizeCallback(; cfl=0.5)
+analyse = AliveCallback(; alive_interval=10, analysis_interval=100)
+cb = CallbackSet(dt_adapt, analyse)
 
-sol = solve(ode, SSPRK33(),dt = dt_adapt(ode),callback= cb)
\ No newline at end of file
+sol = solve(ode, SSPRK33(); dt=dt_adapt(ode), callback=cb)
diff --git a/exemples/Model1DOrd2/exemple.jl b/exemples/Model1DOrd2/exemple.jl
index 712b685..a2f8b0b 100644
--- a/exemples/Model1DOrd2/exemple.jl
+++ b/exemples/Model1DOrd2/exemple.jl
@@ -2,48 +2,40 @@ using Trixi
 using BloodFlowTrixi
 using OrdinaryDiffEq
 
-eq = BloodFlowEquations1D(; h = 0.1)
+eq = BloodFlowEquations1D(; h=0.1)
 eq_ord2 = BloodFlowEquations1DOrd2(eq)
 
-mesh = TreeMesh(0.0, 40.0, 
-    initial_refinement_level = 4, 
-    n_cells_max = 10^4, 
-    periodicity = false)
+mesh = TreeMesh(0.0, 40.0; initial_refinement_level=4, n_cells_max=10^4, periodicity=false)
 
-bc_hypo = (;
-    x_neg = boundary_condition_pressure_in,
-    x_pos = Trixi.BoundaryConditionDoNothing()
-)
+bc_hypo = (; x_neg=boundary_condition_pressure_in, x_pos=Trixi.BoundaryConditionDoNothing())
 
 bc_parab = (;
-    x_neg = BoundaryConditionNeumann((x,t,eq) -> SVector(0.0,0,0,0,0)),
-    x_pos = BoundaryConditionNeumann((x,t,eq) -> SVector(0.0,0,0,0,0))
+    x_neg=BoundaryConditionNeumann((x, t, eq) -> SVector(0.0, 0, 0, 0, 0)),
+    x_pos=BoundaryConditionNeumann((x, t, eq) -> SVector(0.0, 0, 0, 0, 0)),
 )
 
-solver = DGSEM(polydeg = 2,
-    surface_flux = (flux_lax_friedrichs,flux_nonconservative),
-    volume_integral = VolumeIntegralFluxDifferencing((flux_lax_friedrichs,flux_nonconservative))
-    )
+solver = DGSEM(;
+    polydeg=2,
+    surface_flux=(flux_lax_friedrichs, flux_nonconservative),
+    volume_integral=VolumeIntegralFluxDifferencing((
+        flux_lax_friedrichs, flux_nonconservative
+    )),
+)
 
 semi = SemidiscretizationHyperbolicParabolic(
     mesh,
-    (eq,eq_ord2),
-    initial_condition_simple,
-    source_terms = source_term_simple_ord2,
+    (eq, eq_ord2),
+    initial_condition_simple;
+    source_terms=source_term_simple_ord2,
     solver,
-    boundary_conditions = (bc_hypo,bc_parab)
+    boundary_conditions=(bc_hypo, bc_parab),
 )
 
 tspan = (0.0, 0.3)
 ode = semidiscretize(semi, tspan)
 
-dt_adapt = StepsizeCallback(;cfl=0.5)
-analyse = AliveCallback(
-    alive_interval = 10,
-    analysis_interval = 100
-)
-cb = CallbackSet(
-    dt_adapt,analyse
-)
+dt_adapt = StepsizeCallback(; cfl=0.5)
+analyse = AliveCallback(; alive_interval=10, analysis_interval=100)
+cb = CallbackSet(dt_adapt, analyse)
 
-sol = solve(ode, SSPRK33(),dt = dt_adapt(ode),callback= cb)
\ No newline at end of file
+sol = solve(ode, SSPRK33(); dt=dt_adapt(ode), callback=cb)
diff --git a/exemples/Model2D/diexemple.jl b/exemples/Model2D/diexemple.jl
index 75d3435..640405d 100644
--- a/exemples/Model2D/diexemple.jl
+++ b/exemples/Model2D/diexemple.jl
@@ -6,56 +6,53 @@ using StaticArrays, LinearAlgebra
 using QuadGK
 using LinearAlgebra
 
-eq = BloodFlowEquations2D(; h = 0.1)
-xyz_data = [SA[cos(0.2*si),sin(0.2*si),si] for si in range(0,40,100)]
+eq = BloodFlowEquations2D(; h=0.1)
+xyz_data = [SA[cos(0.2*si), sin(0.2*si), si] for si in range(0, 40, 100)]
 
 curve = interpolate_curve(xyz_data)
-L = curve.t[end-1]
+L = curve.t[end - 1]
 println("curve length : $L")
-BloodFlowTrixi.curvature(s) = norm(DataInterpolations.derivative(curve,s,2))
+BloodFlowTrixi.curvature(s) = norm(DataInterpolations.derivative(curve, s, 2))
 
 mesh = P4estMesh(
-    (1,2),
-    polydeg = 1,
-    periodicity = (true,false),
-    coordinates_min = (0.0,0.0),
-    coordinates_max = (2*pi,L),
-    initial_refinement_level = 4
+    (1, 2);
+    polydeg=1,
+    periodicity=(true, false),
+    coordinates_min=(0.0, 0.0),
+    coordinates_max=(2*pi, L),
+    initial_refinement_level=4,
 )
 
-bc = (;
-    y_neg = boundary_condition_pressure_in,
-    y_pos = Trixi.BoundaryConditionDoNothing()
-)
+bc = (; y_neg=boundary_condition_pressure_in, y_pos=Trixi.BoundaryConditionDoNothing())
 
-solver = DGSEM(polydeg = 1,
-    surface_flux = (flux_lax_friedrichs,flux_nonconservative),
-    volume_integral = VolumeIntegralFluxDifferencing((flux_lax_friedrichs,flux_nonconservative))
-    )
+solver = DGSEM(;
+    polydeg=1,
+    surface_flux=(flux_lax_friedrichs, flux_nonconservative),
+    volume_integral=VolumeIntegralFluxDifferencing((
+        flux_lax_friedrichs, flux_nonconservative
+    )),
+)
 
 semi = SemidiscretizationHyperbolic(
     mesh,
     eq,
     initial_condition_simple,
-    source_terms = source_term_simple,
-    solver,
-    boundary_conditions = bc
+    solver;
+    source_terms=source_term_simple,
+    boundary_conditions=bc,
 )
 
 tspan = (0.0, 0.3)
 ode = semidiscretize(semi, tspan)
 
-dt_adapt = StepsizeCallback(;cfl=0.5)
-analyse = AliveCallback(
-    alive_interval = 10,
-    analysis_interval = 100
-)
-cb = CallbackSet(
-    dt_adapt,analyse
-)
+dt_adapt = StepsizeCallback(; cfl=0.5)
+analyse = AliveCallback(; alive_interval=10, analysis_interval=100)
+cb = CallbackSet(dt_adapt, analyse)
 
-sol = solve(ode, SSPRK33(),dt = dt_adapt(ode),callback= cb,saveat = 0.03,save_everystep = false)
+sol = solve(
+    ode, SSPRK33(); dt=dt_adapt(ode), callback=cb, saveat=0.03, save_everystep=false
+)
 
 # artery center-line
 
-res = get3DData(eq,xyz_data,semi,sol,1)
+res = get3DData(eq, xyz_data, semi, sol, 1)
diff --git a/exemples/Model2D/exemple.jl b/exemples/Model2D/exemple.jl
index 7f5b504..16107c2 100644
--- a/exemples/Model2D/exemple.jl
+++ b/exemples/Model2D/exemple.jl
@@ -2,46 +2,41 @@ using Trixi
 using BloodFlowTrixi
 using OrdinaryDiffEq
 
-eq = BloodFlowEquations2D(; h = 0.1)
+eq = BloodFlowEquations2D(; h=0.1)
 
 mesh = P4estMesh(
-    (1,2),
-    polydeg = 1,
-    periodicity = (true,false),
-    coordinates_min = (0.0,0.0),
-    coordinates_max = (2*pi,40.0),
-    initial_refinement_level = 4
+    (1, 2);
+    polydeg=1,
+    periodicity=(true, false),
+    coordinates_min=(0.0, 0.0),
+    coordinates_max=(2*pi, 40.0),
+    initial_refinement_level=4,
 )
 
-bc = (;
-    y_neg = boundary_condition_pressure_in,
-    y_pos = Trixi.BoundaryConditionDoNothing()
-)
+bc = (; y_neg=boundary_condition_pressure_in, y_pos=Trixi.BoundaryConditionDoNothing())
 
-solver = DGSEM(polydeg = 1,
-    surface_flux = (flux_lax_friedrichs,flux_nonconservative),
-    volume_integral = VolumeIntegralFluxDifferencing((flux_lax_friedrichs,flux_nonconservative))
-    )
+solver = DGSEM(;
+    polydeg=1,
+    surface_flux=(flux_lax_friedrichs, flux_nonconservative),
+    volume_integral=VolumeIntegralFluxDifferencing((
+        flux_lax_friedrichs, flux_nonconservative
+    )),
+)
 
 semi = SemidiscretizationHyperbolic(
     mesh,
     eq,
     initial_condition_simple,
-    source_terms = source_term_simple,
-    solver,
-    boundary_conditions = bc
+    solver;
+    source_terms=source_term_simple,
+    boundary_conditions=bc,
 )
 
 tspan = (0.0, 0.3)
 ode = semidiscretize(semi, tspan)
 
-dt_adapt = StepsizeCallback(;cfl=0.5)
-analyse = AliveCallback(
-    alive_interval = 10,
-    analysis_interval = 100
-)
-cb = CallbackSet(
-    dt_adapt,analyse
-)
+dt_adapt = StepsizeCallback(; cfl=0.5)
+analyse = AliveCallback(; alive_interval=10, analysis_interval=100)
+cb = CallbackSet(dt_adapt, analyse)
 
-sol = solve(ode, SSPRK33(),dt = dt_adapt(ode),callback= cb)
+sol = solve(ode, SSPRK33(); dt=dt_adapt(ode), callback=cb)
diff --git a/ext/BloodFlowTrixiDataInterpolationsExt.jl b/ext/BloodFlowTrixiDataInterpolationsExt.jl
index b587dff..f85c125 100644
--- a/ext/BloodFlowTrixiDataInterpolationsExt.jl
+++ b/ext/BloodFlowTrixiDataInterpolationsExt.jl
@@ -1,41 +1,53 @@
 module BloodFlowTrixiDataInterpolationsExt
-    if isdefined(Base, :get_extension) 
-        using BloodFlowTrixi
-        using DataInterpolations
-    else
-        using ..BloodFlowTrixi
-        using ..DataInterpolations
-    end
-    using StaticArrays, LinearAlgebra
-    using ForwardDiff
-    function BloodFlowTrixi.interpolate_curve(curve_data::AbstractArray)
-        N = length(curve_data)
-        quadinterp = QuadraticSpline(curve_data,range(0,1,N))
-        curve = SmoothArcLengthInterpolation(quadinterp;m=N,in_place=false)
-        return curve
-    end
-    function BloodFlowTrixi.get3DData(eq::BloodFlowEquations2D,curve_data::AbstractArray,semi,sol,time_index ::Int = 1;vtk ::Bool=false,out ::T="./datas") where {T<:AbstractString}
-        curve = interpolate_curve(curve_data)
-        tanj(s) = ForwardDiff.derivative(curve,s)
-        function nor(s) 
-            res= ForwardDiff.derivative(tanj,s)
-            n = norm(res)
-            if n ≈ 0
-                a,b,c = tanj(s)
-                # return a any normal vector
-                if a != 0
-                    return SA[-b,a,0]/sqrt(a^2+b^2)
-                elseif b != 0
-                    return SA[-b,a,0]/sqrt(a^2+b^2)
-                else
-                    return SA[-sign(c),0,0]
-                end
+if isdefined(Base, :get_extension)
+    using BloodFlowTrixi
+    using DataInterpolations
+else
+    using ..BloodFlowTrixi
+    using ..DataInterpolations
+end
+using StaticArrays, LinearAlgebra
+using ForwardDiff
+function BloodFlowTrixi.interpolate_curve(curve_data::AbstractArray)
+    N = length(curve_data)
+    quadinterp = QuadraticSpline(curve_data, range(0, 1, N))
+    curve = SmoothArcLengthInterpolation(quadinterp; m=N, in_place=false)
+    return curve
+end
+function BloodFlowTrixi.get3DData(
+    eq::BloodFlowEquations2D,
+    curve_data::AbstractArray,
+    semi,
+    sol,
+    time_index::Int=1;
+    vtk::Bool=false,
+    out::T="./datas",
+) where {T<:AbstractString}
+    curve = interpolate_curve(curve_data)
+    tanj(s) = ForwardDiff.derivative(curve, s)
+    function nor(s)
+        res = ForwardDiff.derivative(tanj, s)
+        n = norm(res)
+        if n ≈ 0
+            a, b, c = tanj(s)
+            # return a any normal vector
+            if a != 0
+                return SA[-b, a, 0]/sqrt(a^2+b^2)
+            elseif b != 0
+                return SA[-b, a, 0]/sqrt(a^2+b^2)
+            else
+                return SA[-sign(c), 0, 0]
             end
-            return res/n
         end
-        ∧(v,w) = SA[v[2]*w[3]-v[3]*w[2],v[3]*w[1]-v[1]*w[3],v[1]*w[2]-v[2]*w[1]]
-        er(theta,s) = cos(theta).*nor(s) .+ sin(theta).*∧(tanj(s),nor(s))
-        return BloodFlowTrixi.get3DData(eq,s->curve(s),er,semi,sol,time_index;vtk=vtk,out=out)
+        return res/n
     end
-    export get3DData,interpolate_curve
-end
\ No newline at end of file
+    ∧(v, w) = SA[
+        v[2] * w[3] - v[3] * w[2], v[3] * w[1] - v[1] * w[3], v[1] * w[2] - v[2] * w[1]
+    ]
+    er(theta, s) = cos(theta) .* nor(s) .+ sin(theta) .* ∧(tanj(s), nor(s))
+    return BloodFlowTrixi.get3DData(
+        eq, s->curve(s), er, semi, sol, time_index; vtk=vtk, out=out
+    )
+end
+export get3DData, interpolate_curve
+end
diff --git a/src/1DModel/1dmodel.jl b/src/1DModel/1dmodel.jl
index 64763c3..f106f8e 100644
--- a/src/1DModel/1dmodel.jl
+++ b/src/1DModel/1dmodel.jl
@@ -15,14 +15,14 @@ The governing equations are given by
 """
 struct BloodFlowEquations1D{T<:Real} <: AbstractBloodFlowEquations{1,4}
     # constant coefficients
-    h ::T      # Wall thickness
+    h::T      # Wall thickness
     rho::T     # Fluid density
     xi::T      # Poisson's ratio
     nu::T      # Viscosity coefficient
 end
 
-function BloodFlowEquations1D(;h,rho=1.0,xi=0.25,nu=0.04)
-    return BloodFlowEquations1D(h,rho,xi,nu)
+function BloodFlowEquations1D(; h, rho=1.0, xi=0.25, nu=0.04)
+    return BloodFlowEquations1D(h, rho, xi, nu)
 end
 
 Trixi.have_nonconservative_terms(::BloodFlowEquations1D) = Trixi.True()
@@ -40,14 +40,14 @@ Computes the flux vector for the conservation laws of the blood flow model.
 ### Returns
 Flux vector as an `SVector`.
 """
-function Trixi.flux(u, orientation::Integer,eq::BloodFlowEquations1D)
+function Trixi.flux(u, orientation::Integer, eq::BloodFlowEquations1D)
     # up = cons2prim(u,eq)
-    P = pressure(u,eq)
-    a,Q,E,A0 = u 
+    P = pressure(u, eq)
+    a, Q, E, A0 = u
     A = a+A0
     f1 = Q
     f2 = Q^2/A+A*P
-    return SVector(f1,f2,0,0)
+    return SVector(f1, f2, 0, 0)
 end
 
 @doc raw"""
@@ -64,17 +64,17 @@ Computes the non-conservative flux for the model, used for handling discontinuit
 ### Returns
 Non-conservative flux vector.
 """
-function flux_nonconservative(u_ll,u_rr,orientation::Integer,eq::BloodFlowEquations1D)
+function flux_nonconservative(u_ll, u_rr, orientation::Integer, eq::BloodFlowEquations1D)
     T = eltype(u_ll)
-    p_ll = pressure(u_ll,eq)
-    p_rr = pressure(u_rr,eq)
+    p_ll = pressure(u_ll, eq)
+    p_rr = pressure(u_rr, eq)
     pmean = (p_ll+p_rr)/2
-    a_ll,_,_,A0_ll = u_ll
-    a_rr,_,_,A0_rr = u_rr
+    a_ll, _, _, A0_ll = u_ll
+    a_rr, _, _, A0_rr = u_rr
     A_ll = a_ll + A0_ll
     A_rr = a_rr + A0_rr
     Ajump = A_rr - A_ll
-    return SVector(zero(T),-pmean*Ajump,0,0)
+    return SVector(zero(T), -pmean*Ajump, 0, 0)
 end
 
 @doc raw"""
@@ -91,16 +91,18 @@ Calculates the maximum absolute speed for wave propagation in the blood flow mod
 ### Returns
 Maximum absolute speed.
 """
-function Trixi.max_abs_speed_naive(u_ll,u_rr,orientation::Integer,eq ::BloodFlowEquations1D)
-    a_ll,Q_ll,E_ll,A0_ll = u_ll
-    a_rr,Q_rr,E_rr,A0_rr = u_rr
+function Trixi.max_abs_speed_naive(
+    u_ll, u_rr, orientation::Integer, eq::BloodFlowEquations1D
+)
+    a_ll, Q_ll, E_ll, A0_ll = u_ll
+    a_rr, Q_rr, E_rr, A0_rr = u_rr
     A_ll = a_ll + A0_ll
     A_rr = a_rr + A0_rr
-    pp_ll = pressure_der(u_ll,eq)
-    pp_rr = pressure_der(u_rr,eq)
+    pp_ll = pressure_der(u_ll, eq)
+    pp_rr = pressure_der(u_rr, eq)
     w_ll = Q_ll/A_ll
     w_rr = Q_rr/A_rr
-    return max(abs(w_ll),abs(w_rr))+max(sqrt(A_ll*pp_ll),sqrt(A_rr*pp_rr))
+    return max(abs(w_ll), abs(w_rr))+max(sqrt(A_ll*pp_ll), sqrt(A_rr*pp_rr))
 end
 
 @doc raw"""
@@ -114,14 +116,13 @@ Computes the maximum absolute speed for wave propagation in the model.
 Maximum absolute speed as a scalar value.
 """
 
-function Trixi.max_abs_speeds(u,eq::BloodFlowEquations1D)
-    a,Q,E,A0 = u 
+function Trixi.max_abs_speeds(u, eq::BloodFlowEquations1D)
+    a, Q, E, A0 = u
     A = a+A0
-    pp= pressure_der(u,eq)
+    pp = pressure_der(u, eq)
     return abs(Q/A) + sqrt(A*pp)
 end
 
-
 @doc raw"""
     (dissipation::Trixi.DissipationLocalLaxFriedrichs)(u_ll, u_rr, orientation_or_normal_direction, eq::BloodFlowEquations1D)
 
@@ -136,13 +137,12 @@ Calculates the dissipation term using the Local Lax-Friedrichs method.
 ### Returns
 Dissipation vector.
 """
-function (dissipation::Trixi.DissipationLocalLaxFriedrichs)(u_ll, u_rr,
-    orientation_or_normal_direction,
-    eq::BloodFlowEquations1D)
-    λ = dissipation.max_abs_speed(u_ll, u_rr, orientation_or_normal_direction,
-    eq)
+function (dissipation::Trixi.DissipationLocalLaxFriedrichs)(
+    u_ll, u_rr, orientation_or_normal_direction, eq::BloodFlowEquations1D
+)
+    λ = dissipation.max_abs_speed(u_ll, u_rr, orientation_or_normal_direction, eq)
     diss = -0.5f0 .* abs(λ) .* (u_rr .- u_ll)
-    return SVector(diss[1], diss[2],0,0)
+    return SVector(diss[1], diss[2], 0, 0)
 end
 
 include("./variables.jl")
@@ -150,4 +150,4 @@ include("./bc1d.jl")
 include("./Test_Cases/pressure_in.jl")
 include("./Test_Cases/convergence_test.jl")
 include("./Ord2/1dmodelord2.jl")
-include("./viz.jl")
\ No newline at end of file
+include("./viz.jl")
diff --git a/src/1DModel/Ord2/1dmodelord2.jl b/src/1DModel/Ord2/1dmodelord2.jl
index 2b491cd..ccf3bfc 100644
--- a/src/1DModel/Ord2/1dmodelord2.jl
+++ b/src/1DModel/Ord2/1dmodelord2.jl
@@ -1,67 +1,108 @@
-struct BloodFlowEquations1DOrd2{E} <: Trixi.AbstractEquationsParabolic{1, 4, GradientVariablesConservative}
-    model1d ::E
+struct BloodFlowEquations1DOrd2{E} <:
+       Trixi.AbstractEquationsParabolic{1,4,GradientVariablesConservative}
+    model1d::E
+end
+function Trixi.varnames(mapin, eq::BloodFlowTrixi.BloodFlowEquations1DOrd2)
+    Trixi.varnames(mapin, eq.model1d)
 end
-Trixi.varnames(mapin,eq::BloodFlowTrixi.BloodFlowEquations1DOrd2) = Trixi.varnames(mapin,eq.model1d)
 
-function Trixi.flux(u,gradients,orientation::Int,eq_parab ::BloodFlowEquations1DOrd2)
+function Trixi.flux(u, gradients, orientation::Int, eq_parab::BloodFlowEquations1DOrd2)
     dudx = gradients
-    a,Q,_,A0 = u
+    a, Q, _, A0 = u
     A = a+A0
-    val = -3*eq_parab.model1d.nu * (-(dudx[1] + dudx[4])*Q/A + dudx[2])
-    return SVector(0.0,val,0,0)
+    val = -3 * eq_parab.model1d.nu * (-(dudx[1] + dudx[4])*Q/A + dudx[2])
+    return SVector(0.0, val, 0, 0)
 end
 
 function source_term_simple_ord2(u, x, t, eq::BloodFlowEquations1D)
     res = source_term_simple(u, x, t, eq)
-    k = friction(u,x,eq)
-    R = radius(u,eq)
-    return SVector(res[1],res[2]/(1-R*k/(4*eq.nu)),res[3],res[4])
+    k = friction(u, x, eq)
+    R = radius(u, eq)
+    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,
+@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)
+    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,
+@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
+    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,
-    orientation_or_normal,direction,
-    x, t,
+@inline function (boundary_condition::Trixi.BoundaryConditionDirichlet)(
+    flux_inner,
+    u_inner,
+    orientation_or_normal,
+    direction,
+    x,
+    t,
     operator_type::Trixi.Gradient,
-    equations_parabolic::BloodFlowEquations1DOrd2)
-return boundary_condition.boundary_value_function(x, t, equations_parabolic)
+    equations_parabolic::BloodFlowEquations1DOrd2,
+)
+    return boundary_condition.boundary_value_function(x, t, equations_parabolic)
 end
 
-@inline function (boundary_condition::Trixi.BoundaryConditionDirichlet)(flux_inner, u_inner,
-    orientation_or_normal,direction,
-    x, t,
+@inline function (boundary_condition::Trixi.BoundaryConditionDirichlet)(
+    flux_inner,
+    u_inner,
+    orientation_or_normal,
+    direction,
+    x,
+    t,
     operator_type::Trixi.Divergence,
-    equations_parabolic::BloodFlowEquations1DOrd2)
-return flux_inner
+    equations_parabolic::BloodFlowEquations1DOrd2,
+)
+    return flux_inner
 end
 
-@inline function (boundary_condition::Trixi.BoundaryConditionNeumann)(flux_inner, u_inner,
-  orientation_or_normal,direction,
-  x, t,
-  operator_type::Trixi.Divergence,
-  equations_parabolic::BloodFlowEquations1DOrd2)
-return boundary_condition.boundary_normal_flux_function(x, t, equations_parabolic)
+@inline function (boundary_condition::Trixi.BoundaryConditionNeumann)(
+    flux_inner,
+    u_inner,
+    orientation_or_normal,
+    direction,
+    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,
-  orientation_or_normal,direction,
-  x, t,
-  operator_type::Trixi.Gradient,
-  equations_parabolic::BloodFlowEquations1DOrd2)
-return flux_inner
-end
\ No newline at end of file
+@inline function (boundary_condition::Trixi.BoundaryConditionNeumann)(
+    flux_inner,
+    u_inner,
+    orientation_or_normal,
+    direction,
+    x,
+    t,
+    operator_type::Trixi.Gradient,
+    equations_parabolic::BloodFlowEquations1DOrd2,
+)
+    return flux_inner
+end
diff --git a/src/1DModel/Test_Cases/convergence_test.jl b/src/1DModel/Test_Cases/convergence_test.jl
index a8b438c..0e3a43f 100644
--- a/src/1DModel/Test_Cases/convergence_test.jl
+++ b/src/1DModel/Test_Cases/convergence_test.jl
@@ -55,7 +55,7 @@ This function is useful for evaluating the correctness of source term handling i
 function Trixi.source_terms_convergence_test(u, x, t, eq::BloodFlowEquations1D)
     T = eltype(u)
     A0 = u[4]
-    s1 = pi * t * cospi(x[1] * t) |> T
+    s1 = T(pi * t * cospi(x[1] * t))
     # k = friction(u, x, eq)
     # R = radius(u, eq)
     s2 = pi * x[1] * cospi(x[1] * t) + pi * t * cospi(x[1] * t) * sinpi(x[1] * t) / A0
diff --git a/src/1DModel/Test_Cases/pressure_in.jl b/src/1DModel/Test_Cases/pressure_in.jl
index e11d367..947ef1e 100644
--- a/src/1DModel/Test_Cases/pressure_in.jl
+++ b/src/1DModel/Test_Cases/pressure_in.jl
@@ -75,31 +75,28 @@ P_{in} = \begin{cases}
 ```
 The corresponding inflow area `A_{in}` is computed using the inverse pressure relation, and the boundary state is constructed accordingly.
 """
-function boundary_condition_pressure_in(u_inner, orientation_or_normal,
-direction,
-x, t,
-surface_flux_function,
-eq::BloodFlowEquations1D)
+function boundary_condition_pressure_in(
+    u_inner,
+    orientation_or_normal,
+    direction,
+    x,
+    t,
+    surface_flux_function,
+    eq::BloodFlowEquations1D,
+)
     Pin = ifelse(t < 0.125, 2e4 * sinpi(t / 0.125)^2, 0.0)
     Ain = inv_pressure(Pin, u_inner, eq)
     A0in = u_inner[4]
     ain = Ain - A0in
-    u_boundary =  SVector(
-        ain,
-        u_inner[2],
-        u_inner[3],
-        u_inner[4]
-    )
+    u_boundary = SVector(ain, u_inner[2], u_inner[3], u_inner[4])
     # calculate the boundary flux
     if iseven(direction) # u_inner is "left" of boundary, u_boundary is "right" of boundary
-        flux1 = surface_flux_function[1](u_inner, u_boundary, orientation_or_normal,
-        eq)
-        flux2 = surface_flux_function[2](u_inner, u_boundary, orientation_or_normal,
-        eq)
+        flux1 = surface_flux_function[1](u_inner, u_boundary, orientation_or_normal, eq)
+        flux2 = surface_flux_function[2](u_inner, u_boundary, orientation_or_normal, eq)
     else # u_boundary is "left" of boundary, u_inner is "right" of boundary
-        flux1 = surface_flux_function[1](u_boundary, u_inner, orientation_or_normal,eq)
-        flux2 = surface_flux_function[2](u_boundary, u_inner, orientation_or_normal,eq)
+        flux1 = surface_flux_function[1](u_boundary, u_inner, orientation_or_normal, eq)
+        flux2 = surface_flux_function[2](u_boundary, u_inner, orientation_or_normal, eq)
     end
 
-    return flux1,flux2
+    return flux1, flux2
 end
diff --git a/src/1DModel/bc1d.jl b/src/1DModel/bc1d.jl
index c0116d5..920af51 100644
--- a/src/1DModel/bc1d.jl
+++ b/src/1DModel/bc1d.jl
@@ -15,24 +15,24 @@ Implements the outflow boundary condition, assuming that there is no reflection
 ### Returns
 Computed boundary flux.
 """
-function boundary_condition_outflow(u_inner, orientation_or_normal,       
+function boundary_condition_outflow(
+    u_inner,
+    orientation_or_normal,
     direction,
-    x, t,
+    x,
+    t,
     surface_flux_function,
-    eq::BloodFlowEquations1D)
-        # calculate the boundary flux
-        if iseven(direction) # u_inner is "left" of boundary, u_boundary is "right" of boundary
-            flux1 = surface_flux_function[1](u_inner, u_inner, orientation_or_normal,
-            eq)
-            flux2 = surface_flux_function[2](u_inner, u_inner, orientation_or_normal,
-            eq)
-            else # u_inner is "left" of boundary, u_inner is "right" of boundary
-            flux1 = surface_flux_function[1](u_inner, u_inner, orientation_or_normal,
-            eq)
-            flux2 = surface_flux_function[2](u_inner, u_inner, orientation_or_normal,
-            eq)
-            end
-    return flux1,flux2
+    eq::BloodFlowEquations1D,
+)
+    # calculate the boundary flux
+    if iseven(direction) # u_inner is "left" of boundary, u_boundary is "right" of boundary
+        flux1 = surface_flux_function[1](u_inner, u_inner, orientation_or_normal, eq)
+        flux2 = surface_flux_function[2](u_inner, u_inner, orientation_or_normal, eq)
+    else # u_inner is "left" of boundary, u_inner is "right" of boundary
+        flux1 = surface_flux_function[1](u_inner, u_inner, orientation_or_normal, eq)
+        flux2 = surface_flux_function[2](u_inner, u_inner, orientation_or_normal, eq)
+    end
+    return flux1, flux2
 end
 
 @doc raw"""
@@ -52,29 +52,26 @@ Implements a slip wall boundary condition where the normal component of velocity
 ### Returns
 Computed boundary flux at the slip wall.
 """
-function boundary_condition_slip_wall(u_inner, orientation_or_normal,       
+function boundary_condition_slip_wall(
+    u_inner,
+    orientation_or_normal,
     direction,
-    x, t,
+    x,
+    t,
     surface_flux_function,
-    eq::BloodFlowEquations1D)
+    eq::BloodFlowEquations1D,
+)
     # create the "external" boundary solution state
-    u_boundary = SVector(u_inner[1],
-    -u_inner[2],
-    u_inner[3],
-    u_inner[4])
+    u_boundary = SVector(u_inner[1], -u_inner[2], u_inner[3], u_inner[4])
 
     # calculate the boundary flux
     if iseven(direction) # u_inner is "left" of boundary, u_boundary is "right" of boundary
-    flux1 = surface_flux_function[1](u_inner, u_boundary, orientation_or_normal,
-    eq)
-    flux2 = surface_flux_function[2](u_inner, u_boundary, orientation_or_normal,
-    eq)
+        flux1 = surface_flux_function[1](u_inner, u_boundary, orientation_or_normal, eq)
+        flux2 = surface_flux_function[2](u_inner, u_boundary, orientation_or_normal, eq)
     else # u_boundary is "left" of boundary, u_inner is "right" of boundary
-    flux1 = surface_flux_function[1](u_boundary, u_inner, orientation_or_normal,
-    eq)
-    flux2 = surface_flux_function[2](u_boundary, u_inner, orientation_or_normal,
-    eq)
+        flux1 = surface_flux_function[1](u_boundary, u_inner, orientation_or_normal, eq)
+        flux2 = surface_flux_function[2](u_boundary, u_inner, orientation_or_normal, eq)
     end
 
-    return flux1,flux2
-end
\ No newline at end of file
+    return flux1, flux2
+end
diff --git a/src/1DModel/variables.jl b/src/1DModel/variables.jl
index 10810c0..2dd9a05 100644
--- a/src/1DModel/variables.jl
+++ b/src/1DModel/variables.jl
@@ -9,7 +9,7 @@ Returns the variable names corresponding to the conserved variables in the blood
 ### Returns
 A tuple of variable names: `("a", "Q", "E", "A0")`.
 """
-Trixi.varnames(::typeof(cons2cons),::BloodFlowEquations1D) = ("a","Q","E","A0")
+Trixi.varnames(::typeof(cons2cons), ::BloodFlowEquations1D) = ("a", "Q", "E", "A0")
 
 @doc raw"""
     Trixi.varnames(::typeof(cons2prim), ::BloodFlowEquations1D)
@@ -22,7 +22,7 @@ Returns the variable names corresponding to the primitive variables in the blood
 ### Returns
 A tuple of variable names: `("A", "w", "P", "A0", "P")`.
 """
-Trixi.varnames(::typeof(cons2prim),::BloodFlowEquations1D) = ("A","w","P","A0","P")
+Trixi.varnames(::typeof(cons2prim), ::BloodFlowEquations1D) = ("A", "w", "P", "A0", "P")
 
 @doc raw"""
     Trixi.varnames(::typeof(cons2entropy), ::BloodFlowEquations1D)
@@ -35,7 +35,7 @@ Returns the variable names corresponding to the entropy variables in the blood f
 ### Returns
 A tuple of variable names: `("A", "w", "En", "A0", "P")`.
 """
-Trixi.varnames(::typeof(cons2entropy),::BloodFlowEquations1D) = ("A","w","En","A0","P")
+Trixi.varnames(::typeof(cons2entropy), ::BloodFlowEquations1D) = ("A", "w", "En", "A0", "P")
 
 @doc raw"""
     Trixi.cons2prim(u, eq::BloodFlowEquations1D)
@@ -49,12 +49,12 @@ Converts the conserved variables to primitive variables.
 ### Returns
 Primitive variable vector.
 """
-function Trixi.cons2prim(u,eq::BloodFlowEquations1D)
-    a,Q,E,A0 = u
-    P = pressure(u,eq)
+function Trixi.cons2prim(u, eq::BloodFlowEquations1D)
+    a, Q, E, A0 = u
+    P = pressure(u, eq)
     A = a+A0
     w = Q/A
-    return SVector(A,w,P,A0)
+    return SVector(A, w, P, A0)
 end
 
 @doc raw"""
@@ -68,13 +68,13 @@ Converts the conserved variables to entropy variables.
 ### Returns
 Entropy variable vector.
 """
-function Trixi.cons2entropy(u,eq::BloodFlowEquations1D)
-    a,Q,E,A0 = u
-    P = pressure(u,eq)
+function Trixi.cons2entropy(u, eq::BloodFlowEquations1D)
+    a, Q, E, A0 = u
+    P = pressure(u, eq)
     A = a+A0
     w = Q/A
-    En = entropy(u,eq)
-    return SVector(A,w,En,A0)
+    En = entropy(u, eq)
+    return SVector(A, w, En, A0)
 end
 
 @doc raw"""
@@ -89,15 +89,14 @@ Converts the primitive variables to conserved variables.
 ### Returns
 Conserved variable vector.
 """
-function Trixi.prim2cons(u,eq::BloodFlowEquations1D)
-    A,w,P,A0 = u
+function Trixi.prim2cons(u, eq::BloodFlowEquations1D)
+    A, w, P, A0 = u
     a = A-A0
     Q = w*A
     E = P/sqrt(pi)*A0/(sqrt(A)-sqrt(A0))*(1-eq.xi^2)/eq.h
-    return SVector(a,Q,E,A0)
+    return SVector(a, Q, E, A0)
 end
 
-
 @doc raw"""
     friction(u, x, eq::BloodFlowEquations1D)
 
@@ -111,12 +110,11 @@ Calculates the friction term for the blood flow equations, which represents visc
 ### Returns
 Friction coefficient as a scalar.
 """
-function friction(u,x,eq::BloodFlowEquations1D)
-    R=radius(u,eq) 
+function friction(u, x, eq::BloodFlowEquations1D)
+    R=radius(u, eq)
     return eltype(u)(-11*eq.nu/R)
 end
 
-
 @doc raw"""
     pressure(u, eq::BloodFlowEquations1D)
 
@@ -129,7 +127,7 @@ Computes the pressure given the state vector based on the compliance of the arte
 ### Returns
 Pressure as a scalar.
 """
-function pressure(u,eq::BloodFlowEquations1D)
+function pressure(u, eq::BloodFlowEquations1D)
     T = eltype(u)
     A = u[1]+u[4]
     E = u[3]
@@ -152,7 +150,7 @@ Computes the radius of the artery based on the cross-sectional area.
 ### Returns
 Radius as a scalar.
 """
-function radius(u,eq::BloodFlowEquations1D)
+function radius(u, eq::BloodFlowEquations1D)
     return sqrt((u[1]+u[4])/pi)
 end
 
@@ -169,7 +167,7 @@ Computes the inverse relation of pressure to cross-sectional area.
 ### Returns
 Cross-sectional area corresponding to the given pressure.
 """
-function inv_pressure(p,u,eq::BloodFlowEquations1D)
+function inv_pressure(p, u, eq::BloodFlowEquations1D)
     T = eltype(u)
     E = u[3]
     A0 = u[4]
@@ -192,7 +190,7 @@ Computes the derivative of pressure with respect to cross-sectional area.
 ### Returns
 Derivative of pressure.
 """
-function pressure_der(u,eq::BloodFlowEquations1D)
+function pressure_der(u, eq::BloodFlowEquations1D)
     T = eltype(u)
     A = u[1]+u[4]
     E = u[3]
@@ -213,12 +211,12 @@ Computes the entropy of the system for the given state vector.
 ### Returns
 Entropy as a scalar value.
 """
-function Trixi.entropy(u,eq::BloodFlowEquations1D)
-    up = cons2prim(u,eq)
-    _,_,E,_ = u
-    A,w,P,A0 = up
+function Trixi.entropy(u, eq::BloodFlowEquations1D)
+    up = cons2prim(u, eq)
+    _, _, E, _ = u
+    A, w, P, A0 = up
     psi = w^2/2+P
     b = E*eq.h/(1-eq.xi^2)
     pt = b*sqrt(pi)/(3*A0)*(A^(3/2)-A0^(3/2))
     return A*psi - pt
-end
\ No newline at end of file
+end
diff --git a/src/1DModel/viz.jl b/src/1DModel/viz.jl
index dd85b75..f11ef29 100644
--- a/src/1DModel/viz.jl
+++ b/src/1DModel/viz.jl
@@ -31,37 +31,45 @@ Named tuple containing:
 - A cylindrical coordinate transformation is applied to represent the vessel cross-section.
 - If `vtk` is `true`, the function writes the data to VTK format using `vtk_grid`.
 """
-function get3DData(eq::BloodFlowEquations1D,semi,sol,time_index ::Int = 1;theta_disc ::Int = 32,vtk ::Bool=false,out ::T="./datas") where T<:AbstractString
+function get3DData(
+    eq::BloodFlowEquations1D,
+    semi,
+    sol,
+    time_index::Int=1;
+    theta_disc::Int=32,
+    vtk::Bool=false,
+    out::T="./datas",
+) where {T<:AbstractString}
     xval_not = semi.cache.elements.node_coordinates[:]
     # Get unique values
     unique_values = unique(xval_not)
     indices = [findfirst(==(val), xval_not) for val in unique_values]
     xval = xval_not[indices]
     soltime = sol[time_index]
-    aval = @view( @view(soltime[1:4:end])[indices] )
-    Qval = @view( @view(soltime[2:4:end])[indices])
-    Eval = @view( @view(soltime[3:4:end])[indices]) 
-    A0val =  @view( @view(soltime[4:4:end])[indices])
-    Pval = map((a,Q,E,A0)->pressure(SA[a,Q,E,A0],eq),aval,Qval,Eval,A0val)
-    theta = range(0,2π,theta_disc)[1:end-1]
-    M(theta,x,R) = (x,R*cos(theta),R*sin(theta))
+    aval = @view(@view(soltime[1:4:end])[indices])
+    Qval = @view(@view(soltime[2:4:end])[indices])
+    Eval = @view(@view(soltime[3:4:end])[indices])
+    A0val = @view(@view(soltime[4:4:end])[indices])
+    Pval = map((a, Q, E, A0)->pressure(SA[a, Q, E, A0], eq), aval, Qval, Eval, A0val)
+    theta = range(0, 2π, theta_disc)[1:(end - 1)]
+    M(theta, x, R) = (x, R*cos(theta), R*sin(theta))
     Typ = eltype(xval)
     s = (theta_disc-1)*length(xval)
-    x = zeros(Typ,s)
-    y = zeros(Typ,s)
-    z = zeros(Typ,s)
-    A = zeros(Typ,s)
-    w = zeros(Typ,s)
-    P = zeros(Typ,s)
+    x = zeros(Typ, s)
+    y = zeros(Typ, s)
+    z = zeros(Typ, s)
+    A = zeros(Typ, s)
+    w = zeros(Typ, s)
+    P = zeros(Typ, s)
     c=1
-    for i in eachindex(xval,aval,Qval,A0val,Pval)
+    for i in eachindex(xval, aval, Qval, A0val, Pval)
         xvali = xval[i]
         Avali = aval[i] + A0val[i]
-        wvali =Qval[i]/Avali 
+        wvali = Qval[i]/Avali
         Pvali = Pval[i]
         Rvali = sqrt(Avali/pi)
         for thetaj in theta
-            xi,yi,zi = M(thetaj,xvali,Rvali)
+            xi, yi, zi = M(thetaj, xvali, Rvali)
             x[c] = xi
             y[c] = yi
             z[c] = zi
@@ -73,12 +81,12 @@ function get3DData(eq::BloodFlowEquations1D,semi,sol,time_index ::Int = 1;theta_
     end
     if vtk
         npoints = length(x)
-        cells = [MeshCell(VTKCellTypes.VTK_VERTEX, (i, )) for i = 1:npoints]
-        vtk_grid(joinpath(out,"./points$time_index"), x, y, z, cells) do vtk
-            vtk["Area", VTKPointData()] = A 
+        cells = [MeshCell(VTKCellTypes.VTK_VERTEX, (i,)) for i in 1:npoints]
+        vtk_grid(joinpath(out, "./points$time_index"), x, y, z, cells) do vtk
+            vtk["Area", VTKPointData()] = A
             vtk["Speed", VTKPointData()] = w
             vtk["Pressure", VTKPointData()] = P
-        end     
+        end
     end
-    return (;x=x,y=y,z=z,A=A,w=w)
-end
\ No newline at end of file
+    return (; x=x, y=y, z=z, A=A, w=w)
+end
diff --git a/src/2DModel/2dmodel.jl b/src/2DModel/2dmodel.jl
index 7f5f752..d3b522f 100644
--- a/src/2DModel/2dmodel.jl
+++ b/src/2DModel/2dmodel.jl
@@ -23,14 +23,14 @@ The governing equations account for conservation of mass and momentum, incorpora
 
 """
 struct BloodFlowEquations2D{T<:Real} <: AbstractBloodFlowEquations{2,5}
-    h ::T      # Wall thickness
+    h::T      # Wall thickness
     rho::T     # Fluid density
     xi::T      # Poisson's ratio
     nu::T      # Viscosity coefficient
 end
 
-function BloodFlowEquations2D(;h,rho=1.0,xi=0.25,nu=0.04)
-    return BloodFlowEquations2D(h,rho,xi,nu)
+function BloodFlowEquations2D(; h, rho=1.0, xi=0.25, nu=0.04)
+    return BloodFlowEquations2D(h, rho, xi, nu)
 end
 
 Trixi.have_nonconservative_terms(::BloodFlowEquations2D) = Trixi.True()
@@ -49,7 +49,7 @@ Flux vector as an `SVector`.
 """
 function Trixi.flux(u, orientation::Integer, eq::BloodFlowEquations2D)
     P = pressure(u, eq) # Compute pressure from state vector
-    a, QRθ, Qs, E, A0 = u 
+    a, QRθ, Qs, E, A0 = u
     QRθ = 0.0
     A = a + A0 # Total cross-sectional area
     if orientation == 1 # Flux in θ-direction
@@ -60,7 +60,7 @@ function Trixi.flux(u, orientation::Integer, eq::BloodFlowEquations2D)
     else # Flux in s-direction
         f1 = Qs
         f2 = QRθ * Qs / A
-        f3 = Qs^2 / A - QRθ^2/(2*A^2)+ A * P
+        f3 = Qs^2 / A - QRθ^2/(2*A^2)+A * P
         return SVector(f1, f2, f3, 0, 0)
     end
 end
@@ -78,18 +78,18 @@ Flux vector as an `SVector`.
 """
 function Trixi.flux(u, normal, eq::BloodFlowEquations2D)
     P = pressure(u, eq) # Compute pressure from state vector
-    a, QRθ, Qs, E, A0 = u 
+    a, QRθ, Qs, E, A0 = u
     A = a + A0 # Total cross-sectional area
     # if normal == 1 # Flux in θ-direction
-        f1 = QRθ / A
-        f2 = QRθ^2 / (2 * A^2) + A * P
-        f3 = QRθ * Qs / A^2
-        fl1 =  SVector(f1, f2, f3, 0, 0)
+    f1 = QRθ / A
+    f2 = QRθ^2 / (2 * A^2) + A * P
+    f3 = QRθ * Qs / A^2
+    fl1 = SVector(f1, f2, f3, 0, 0)
     # else # Flux in s-direction
-        f1 = Qs
-        f2 = QRθ * Qs / A
-        f3 = Qs^2 / A - QRθ^2/(2*A^2)+ A * P
-        fl2 = SVector(f1, f2, f3, 0, 0)
+    f1 = Qs
+    f2 = QRθ * Qs / A
+    f3 = Qs^2 / A - QRθ^2/(2*A^2)+A * P
+    fl2 = SVector(f1, f2, f3, 0, 0)
     return fl1 .* normal[1] .+ fl2 .* normal[2]
 end
 
@@ -149,9 +149,9 @@ function flux_nonconservative(u_ll, u_rr, normal, eq::BloodFlowEquations2D)
     A_rr = a_rr + A0_rr
     Ajump = A_rr - A_ll # Compute jump in area
     # if orientation == 1
-        fn1 =  SVector(zero(T), -pmean * Ajump, 0, 0, 0)
+    fn1 = SVector(zero(T), -pmean * Ajump, 0, 0, 0)
     # else
-        fn2 = SVector(zero(T), 0, -pmean * Ajump, 0, 0)
+    fn2 = SVector(zero(T), 0, -pmean * Ajump, 0, 0)
     # end
     return @. fn1*normal[1] + fn2*normal[2]
 end
@@ -170,7 +170,9 @@ Computes the maximum absolute speed for wave propagation in the 2D blood flow mo
 Maximum absolute speed.
 
 """
-function Trixi.max_abs_speed_naive(u_ll, u_rr, orientation::Integer, eq::BloodFlowEquations2D)
+function Trixi.max_abs_speed_naive(
+    u_ll, u_rr, orientation::Integer, eq::BloodFlowEquations2D
+)
     a_ll, QRθ_ll, Qs_ll, _, A0_ll = u_ll
     a_rr, QRθ_rr, Qs_rr, _, A0_rr = u_rr
     A_ll = a_ll + A0_ll
@@ -178,7 +180,9 @@ function Trixi.max_abs_speed_naive(u_ll, u_rr, orientation::Integer, eq::BloodFl
     pp_ll = pressure_der(u_ll, eq)
     pp_rr = pressure_der(u_rr, eq)
     if orientation == 1
-        return max(max(abs(QRθ_ll)/A_ll^2, abs(QRθ_rr)/A_rr^2), max(sqrt(pp_ll), sqrt(pp_rr)))
+        return max(
+            max(abs(QRθ_ll)/A_ll^2, abs(QRθ_rr)/A_rr^2), max(sqrt(pp_ll), sqrt(pp_rr))
+        )
     else
         ws_ll = Qs_ll / A_ll
         ws_rr = Qs_rr / A_rr
@@ -207,13 +211,13 @@ function Trixi.max_abs_speed_naive(u_ll, u_rr, normal, eq::BloodFlowEquations2D)
     A_rr = a_rr + A0_rr
     pp_ll = pressure_der(u_ll, eq)
     pp_rr = pressure_der(u_rr, eq)
-        ws_ll = Qs_ll / A_ll
-        ws_rr = Qs_rr / A_rr
+    ws_ll = Qs_ll / A_ll
+    ws_rr = Qs_rr / A_rr
     return max(
         abs(ws_ll*normal[2] + sqrt(A_ll*pp_ll)*sqrt(normal[1]^2/A_ll + normal[2]^2)),
         abs(ws_rr*normal[2] + sqrt(A_rr*pp_rr)*sqrt(normal[1]^2/A_rr + normal[2]^2)),
         abs(ws_ll*normal[2] + QRθ_ll/A_ll^2*normal[1]),
-        abs(ws_rr*normal[2] + QRθ_rr/A_rr^2*normal[1])
+        abs(ws_rr*normal[2] + QRθ_rr/A_rr^2*normal[1]),
     )
 end
 
@@ -230,11 +234,11 @@ Tuple containing the maximum absolute speeds in the \( \theta \)- and \( s \)-di
 
 
 """
-function Trixi.max_abs_speeds(u,eq::BloodFlowEquations2D)
-    a,QRθ,Qs,E,A0 = u 
+function Trixi.max_abs_speeds(u, eq::BloodFlowEquations2D)
+    a, QRθ, Qs, E, A0 = u
     A = a+A0
-    pp= pressure_der(u,eq)
-    return max(abs(QRθ/A^2),sqrt(pp)),abs(Qs/A) + sqrt(A*pp)
+    pp = pressure_der(u, eq)
+    return max(abs(QRθ/A^2), sqrt(pp)), abs(Qs/A) + sqrt(A*pp)
 end
 
 @doc raw"""
@@ -251,15 +255,16 @@ Computes the dissipation term using the Local Lax-Friedrichs method for the 2D b
 Dissipation vector.
 """
 
-function (dissipation::Trixi.DissipationLocalLaxFriedrichs)(u_ll, u_rr, orientation_or_normal_direction, eq::BloodFlowEquations2D)
+function (dissipation::Trixi.DissipationLocalLaxFriedrichs)(
+    u_ll, u_rr, orientation_or_normal_direction, eq::BloodFlowEquations2D
+)
     λ = dissipation.max_abs_speed(u_ll, u_rr, orientation_or_normal_direction, eq)
     diss = -0.5 .* abs(λ) .* (u_rr .- u_ll) # Compute dissipation term
     return SVector(diss[1], diss[2], diss[3], 0, 0)
 end
 
-
 include("./variables.jl")
 include("./bc2d.jl")
 include("./Test_Cases/pressure_in.jl")
 include("./Test_Cases/convergence_test.jl")
-include("./viz.jl")
\ No newline at end of file
+include("./viz.jl")
diff --git a/src/2DModel/Test_Cases/convergence_test.jl b/src/2DModel/Test_Cases/convergence_test.jl
index 99a07b3..d7828c8 100644
--- a/src/2DModel/Test_Cases/convergence_test.jl
+++ b/src/2DModel/Test_Cases/convergence_test.jl
@@ -1,21 +1,17 @@
 
-
-
 function Trixi.initial_condition_convergence_test(x, t, eq::BloodFlowEquations2D)
     T = eltype(x)
     R0 = T(1.0)
     A0 = T(R0^2/2)
     E = T(1e7)
-    QRθ = 
-    Qs = T(sinpi(x[1] * t))
-    return SVector(zero(T), QRθ,Qs, E, A0)
+    QRθ = Qs = T(sinpi(x[1] * t))
+    return SVector(zero(T), QRθ, Qs, E, A0)
 end
 
-
 function Trixi.source_terms_convergence_test(u, x, t, eq::BloodFlowEquations2D)
     T = eltype(u)
     A0 = u[4]
-    s1 = pi * t * cospi(x[1] * t) |> T
+    s1 = T(pi * t * cospi(x[1] * t))
     # k = friction(u, x, eq)
     # R = radius(u, eq)
     s2 = pi * x[1] * cospi(x[1] * t) + pi * t * cospi(x[1] * t) * sinpi(x[1] * t) / A0
diff --git a/src/2DModel/Test_Cases/pressure_in.jl b/src/2DModel/Test_Cases/pressure_in.jl
index 6b22d9f..74c4c88 100644
--- a/src/2DModel/Test_Cases/pressure_in.jl
+++ b/src/2DModel/Test_Cases/pressure_in.jl
@@ -17,11 +17,10 @@ function initial_condition_simple(x, t, eq::BloodFlowEquations2D; R0=2.0)
     A0 = T(R0^2 / 2)
     QRθ = T(0.0)
     Qs = T(0.0)
-    E = T(1e7) 
+    E = T(1e7)
     return SVector(zero(T), QRθ, Qs, E, A0)
 end
 
-
 @doc raw"""
     curvature(x)
 
@@ -35,7 +34,6 @@ Curvature as a scalar.
 """
 curvature(x) = typeof(x)(1.0)
 
-
 @doc raw"""
     source_term_simple(u, x, t, eq::BloodFlowEquations2D)
 
@@ -57,16 +55,11 @@ function source_term_simple(u, x, t, eq::BloodFlowEquations2D)
     s1 = zero(T)
     k = friction(u, x, eq)
     R = radius(u, eq)
-    s2 = T(
-        2 * R / 3 * curvature(x[2]) * sin(x[1]) * Qs^2 / A + 3 * R * k * QRθ / A
-    )
-    s3 = T(
-        -2 * R / 3 * curvature(x[2]) * sin(x[1]) * Qs * QRθ / A + R * k * Qs / A
-    )
+    s2 = T(2 * R / 3 * curvature(x[2]) * sin(x[1]) * Qs^2 / A + 3 * R * k * QRθ / A)
+    s3 = T(-2 * R / 3 * curvature(x[2]) * sin(x[1]) * Qs * QRθ / A + R * k * Qs / A)
     return SVector(s1, s2, s3, 0, 0)
 end
 
-
 @doc raw"""
     boundary_condition_pressure_in(u_inner, orientation_or_normal, direction, x, t, surface_flux_function, eq::BloodFlowEquations2D)
 
@@ -84,18 +77,20 @@ Applies an inflow boundary condition with a prescribed pressure for the 2D blood
 ### Returns
 Boundary flux as an `SVector`.
 """
-function boundary_condition_pressure_in(u_inner, orientation_or_normal, direction, x, t, surface_flux_function, eq::BloodFlowEquations2D)
+function boundary_condition_pressure_in(
+    u_inner,
+    orientation_or_normal,
+    direction,
+    x,
+    t,
+    surface_flux_function,
+    eq::BloodFlowEquations2D,
+)
     Pin = ifelse(t < 0.125, 2e4 * sinpi(t / 0.125)^2, 0.0)
     Ain = inv_pressure(Pin, u_inner, eq)
     A0in = u_inner[5]
     ain = Ain - A0in
-    u_boundary = SVector(
-        ain,
-        u_inner[2],
-        u_inner[3],
-        u_inner[4],
-        u_inner[5]
-    )
+    u_boundary = SVector(ain, u_inner[2], u_inner[3], u_inner[4], u_inner[5])
     # Calculate the boundary flux
     if iseven(direction) # u_inner is "left" of boundary, u_boundary is "right" of boundary
         flux1 = surface_flux_function[1](u_inner, u_boundary, orientation_or_normal, eq)
@@ -104,10 +99,9 @@ function boundary_condition_pressure_in(u_inner, orientation_or_normal, directio
         flux1 = surface_flux_function[1](u_boundary, u_inner, orientation_or_normal, eq)
         flux2 = surface_flux_function[2](u_boundary, u_inner, orientation_or_normal, eq)
     end
-    return flux1,flux2
+    return flux1, flux2
 end
 
-
 @doc raw"""
     boundary_condition_pressure_in(u_inner, normal, x, t, surface_flux_function, eq::BloodFlowEquations2D)
 
@@ -124,19 +118,15 @@ Applies an inflow boundary condition with a prescribed pressure for the 2D blood
 ### Returns
 Boundary flux as an `SVector`.
 """
-function boundary_condition_pressure_in(u_inner, normal, x, t, surface_flux_function, eq::BloodFlowEquations2D)
+function boundary_condition_pressure_in(
+    u_inner, normal, x, t, surface_flux_function, eq::BloodFlowEquations2D
+)
     Pin = ifelse(t < 0.125, 2e4 * sinpi(t / 0.125)^2, 0.0)
     A0in = u_inner[5]
     Ain = inv_pressure(Pin, u_inner, eq)
     ain = Ain - A0in
-    u_boundary = SVector(
-        ain,
-        u_inner[2],
-        u_inner[3],
-        u_inner[4],
-        u_inner[5]
-    )
+    u_boundary = SVector(ain, u_inner[2], u_inner[3], u_inner[4], u_inner[5])
     flux1 = surface_flux_function[1](u_inner, u_boundary, normal, eq)
     flux2 = surface_flux_function[2](u_inner, u_boundary, normal, eq)
-    return flux1,flux2
+    return flux1, flux2
 end
diff --git a/src/2DModel/bc2d.jl b/src/2DModel/bc2d.jl
index 1c71a99..6bf1f9e 100644
--- a/src/2DModel/bc2d.jl
+++ b/src/2DModel/bc2d.jl
@@ -15,14 +15,21 @@ Applies an outflow boundary condition for the 2D blood flow model without reflec
 ### Returns
 Boundary flux as an `SVector`.
 """
-function boundary_condition_outflow(u_inner, orientation_or_normal, direction, x, t, surface_flux_function, eq::BloodFlowEquations2D)
+function boundary_condition_outflow(
+    u_inner,
+    orientation_or_normal,
+    direction,
+    x,
+    t,
+    surface_flux_function,
+    eq::BloodFlowEquations2D,
+)
     # Calculate the boundary flux without reflection
     flux1 = surface_flux_function[1](u_inner, u_inner, orientation_or_normal, eq)
     flux2 = surface_flux_function[2](u_inner, u_inner, orientation_or_normal, eq)
-    return flux1,flux2
+    return flux1, flux2
 end
 
-
 @doc raw"""
     boundary_condition_outflow(u_inner, orientation_or_normal, x, t, surface_flux_function, eq::BloodFlowEquations2D)
 
@@ -39,14 +46,15 @@ Applies an outflow boundary condition for the 2D blood flow model without reflec
 ### Returns
 Boundary flux as an `SVector`.
 """
-function boundary_condition_outflow(u_inner, orientation_or_normal, x, t, surface_flux_function, eq::BloodFlowEquations2D)
+function boundary_condition_outflow(
+    u_inner, orientation_or_normal, x, t, surface_flux_function, eq::BloodFlowEquations2D
+)
     # Calculate the boundary flux without reflection
     flux1 = surface_flux_function[1](u_inner, u_inner, orientation_or_normal, eq)
     flux2 = surface_flux_function[2](u_inner, u_inner, orientation_or_normal, eq)
-    return flux1,flux2
+    return flux1, flux2
 end
 
-
 @doc raw"""
     boundary_condition_slip_wall(u_inner, orientation_or_normal, direction, x, t, surface_flux_function, eq::BloodFlowEquations2D)
 
@@ -64,7 +72,15 @@ Applies a slip-wall boundary condition for the 2D blood flow model by reflecting
 ### Returns
 Boundary flux as an `SVector`.
 """
-function boundary_condition_slip_wall(u_inner, orientation_or_normal, direction, x, t, surface_flux_function, eq::BloodFlowEquations2D)
+function boundary_condition_slip_wall(
+    u_inner,
+    orientation_or_normal,
+    direction,
+    x,
+    t,
+    surface_flux_function,
+    eq::BloodFlowEquations2D,
+)
     # Create the external boundary solution state with reflected normal velocity
     u_boundary = SVector(u_inner[1], -u_inner[2], u_inner[3], u_inner[4])
 
@@ -76,5 +92,5 @@ function boundary_condition_slip_wall(u_inner, orientation_or_normal, direction,
         flux1 = surface_flux_function[1](u_boundary, u_inner, orientation_or_normal, eq)
         flux2 = surface_flux_function[2](u_boundary, u_inner, orientation_or_normal, eq)
     end
-    return flux1,flux2
+    return flux1, flux2
 end
diff --git a/src/2DModel/variables.jl b/src/2DModel/variables.jl
index ca64eb9..a9b2faf 100644
--- a/src/2DModel/variables.jl
+++ b/src/2DModel/variables.jl
@@ -12,7 +12,6 @@ Tuple containing the names of the conservative variables.
 """
 Trixi.varnames(::typeof(cons2cons), ::BloodFlowEquations2D) = ("a", "QRθ", "Qs", "E", "A0")
 
-
 @doc raw"""
     Trixi.varnames(::typeof(cons2prim), ::BloodFlowEquations2D)
 
@@ -27,7 +26,6 @@ Tuple containing the names of the primitive variables.
 """
 Trixi.varnames(::typeof(cons2prim), ::BloodFlowEquations2D) = ("A", "wθ", "ws", "P", "A0")
 
-
 @doc raw"""
     Trixi.varnames(::typeof(cons2entropy), ::BloodFlowEquations2D)
 
@@ -40,8 +38,9 @@ Returns the variable names in entropy form for the 2D blood flow model.
 ### Returns
 Tuple containing the names of the entropy variables.
 """
-Trixi.varnames(::typeof(cons2entropy), ::BloodFlowEquations2D) = ("A", "wθ", "ws", "En", "A0")
-
+Trixi.varnames(::typeof(cons2entropy), ::BloodFlowEquations2D) = (
+    "A", "wθ", "ws", "En", "A0"
+)
 
 @doc raw"""
     Trixi.prim2cons(u, eq::BloodFlowEquations2D)
@@ -64,7 +63,6 @@ function Trixi.prim2cons(u, eq::BloodFlowEquations2D)
     return SVector(a, QRθ, Qs, E, A0)
 end
 
-
 @doc raw"""
     Trixi.cons2prim(u, eq::BloodFlowEquations2D)
 
@@ -89,7 +87,6 @@ function Trixi.cons2prim(u, eq::BloodFlowEquations2D)
     return SVector(A, wθ, ws, P, A0)
 end
 
-
 @doc raw"""
     Trixi.cons2entropy(u, eq::BloodFlowEquations2D)
 
@@ -115,7 +112,6 @@ function Trixi.cons2entropy(u, eq::BloodFlowEquations2D)
     return SVector(A, wθ, ws, En, A0)
 end
 
-
 @doc raw"""
     friction(u, x, eq::BloodFlowEquations2D)
 
@@ -134,7 +130,6 @@ function friction(u, x, eq::BloodFlowEquations2D)
     return eltype(u)(-11 * eq.nu / R) # Return friction term based on viscosity and radius
 end
 
-
 @doc raw"""
     pressure(u, eq::BloodFlowEquations2D)
 
@@ -160,7 +155,6 @@ function pressure(u, eq::BloodFlowEquations2D)
     return T(b * (R - R0) / R0^2)
 end
 
-
 @doc raw"""
     radius(u, eq::BloodFlowEquations2D)
 
@@ -201,7 +195,6 @@ function inv_pressure(p, u, eq::BloodFlowEquations2D)
     return T((R0^2 * p / b + R0)^2 / 2)
 end
 
-
 @doc raw"""
     pressure_der(u, eq::BloodFlowEquations2D)
 
@@ -225,7 +218,6 @@ function pressure_der(u, eq::BloodFlowEquations2D)
     return T((b / sqrt(2)) * 0.5 / (sqrt(A) * A0))
 end
 
-
 @doc raw"""
     Trixi.entropy(u, eq::BloodFlowEquations2D)
 
diff --git a/src/2DModel/viz.jl b/src/2DModel/viz.jl
index e74bc38..d63f468 100644
--- a/src/2DModel/viz.jl
+++ b/src/2DModel/viz.jl
@@ -25,38 +25,54 @@ Named tuple containing:
 - `wtheta`: Flow angular velocity at each point.
 - `ws`: Flow axial velocity at each point.
 """
-function get3DData(eq::BloodFlowEquations2D,curve::F1,er::F2,semi,sol,time_index ::Int = 1;vtk ::Bool=false,out ::T="./datas") where {T<:AbstractString,F1<:Function,F2<:Function}
-    thetaval = semi.cache.elements.node_coordinates[1,:,:,:]
-    sval = semi.cache.elements.node_coordinates[2,:,:,:]
+function get3DData(
+    eq::BloodFlowEquations2D,
+    curve::F1,
+    er::F2,
+    semi,
+    sol,
+    time_index::Int=1;
+    vtk::Bool=false,
+    out::T="./datas",
+) where {T<:AbstractString,F1<:Function,F2<:Function}
+    thetaval = semi.cache.elements.node_coordinates[1, :, :, :]
+    sval = semi.cache.elements.node_coordinates[2, :, :, :]
     # Get unique values
     soltime = sol[time_index]
-    aval =@view(soltime[1:5:end]) 
-    Qthval =@view(soltime[2:5:end])
-    Qsval =@view(soltime[3:5:end])
-    Eval =@view(soltime[4:5:end]) 
+    aval = @view(soltime[1:5:end])
+    Qthval = @view(soltime[2:5:end])
+    Qsval = @view(soltime[3:5:end])
+    Eval = @view(soltime[4:5:end])
     A0val = @view(soltime[5:5:end])
-    Pval = map((a,Qth,Qs,E,A0)->BloodFlowTrixi.pressure(SA[a,Qth,Qs,E,A0],eq),aval,Qthval,Qsval,Eval,A0val)
-    M(theta,s,R) = curve(s) .+ R.*er(theta,s)
-    Typ = eltype(er(sval[1],thetaval[1]))
+    Pval = map(
+        (a, Qth, Qs, E, A0)->BloodFlowTrixi.pressure(SA[a, Qth, Qs, E, A0], eq),
+        aval,
+        Qthval,
+        Qsval,
+        Eval,
+        A0val,
+    )
+    M(theta, s, R) = curve(s) .+ R .* er(theta, s)
+    Typ = eltype(er(sval[1], thetaval[1]))
     s = length(thetaval)
-    x = zeros(Typ,s)
-    y = zeros(Typ,s)
-    z = zeros(Typ,s)
-    A = zeros(Typ,s)
-    wtheta = zeros(Typ,s)
-    ws = zeros(Typ,s)
-    P = zeros(Typ,s)
+    x = zeros(Typ, s)
+    y = zeros(Typ, s)
+    z = zeros(Typ, s)
+    A = zeros(Typ, s)
+    wtheta = zeros(Typ, s)
+    ws = zeros(Typ, s)
+    P = zeros(Typ, s)
     c=1
-    for i in eachindex(thetaval,sval,aval,Qthval,Qsval,A0val,Pval)
+    for i in eachindex(thetaval, sval, aval, Qthval, Qsval, A0val, Pval)
         thvali = thetaval[i]
         svali = sval[i]
         Avali = aval[i] + A0val[i]
         Rvali = sqrt(2*Avali)
         wthetavali = Typ(4/3*(Qthval[i]/Rvali)/Avali)
-        wthetavali =Qthval[i]/Avali 
-        wsvali = Qsval[i]/Avali 
+        wthetavali = Qthval[i]/Avali
+        wsvali = Qsval[i]/Avali
         Pvali = Pval[i]
-        xi,yi,zi = M(thvali,svali,Rvali)
+        xi, yi, zi = M(thvali, svali, Rvali)
         x[c] = xi
         y[c] = yi
         z[c] = zi
@@ -68,18 +84,17 @@ function get3DData(eq::BloodFlowEquations2D,curve::F1,er::F2,semi,sol,time_index
     end
     if vtk
         npoints = length(x)
-        cells = [MeshCell(VTKCellTypes.VTK_VERTEX, (i, )) for i = 1:npoints]
-        vtk_grid(joinpath(out,"./points$time_index"), x, y, z, cells) do vtk
-            vtk["Area", VTKPointData()] = A 
+        cells = [MeshCell(VTKCellTypes.VTK_VERTEX, (i,)) for i in 1:npoints]
+        vtk_grid(joinpath(out, "./points$time_index"), x, y, z, cells) do vtk
+            vtk["Area", VTKPointData()] = A
             vtk["Angular_Speed", VTKPointData()] = wtheta
             vtk["Axial", VTKPointData()] = ws
             vtk["Pressure", VTKPointData()] = P
-        end     
+        end
     end
-    return (;x=x,y=y,z=z,A=A,wtheta=wtheta,ws=ws)
+    return (; x=x, y=y, z=z, A=A, wtheta=wtheta, ws=ws)
 end
 
-
 @doc raw"""
     get3DData(eq::BloodFlowEquations2D,curve::F1,tanj::F2,nor::F3,semi,sol,time_index ::Int = 1;vtk ::Bool=false,out ::T="./datas") where {T<:AbstractString,F1<:Function,F2<:Function,F3<:Function}
 
@@ -109,10 +124,22 @@ Named tuple containing:
 - `ws`: Flow axial velocity at each point.
 
 """
-function get3DData(eq::BloodFlowEquations2D,curve::F1,tanj::F2,nor::F3,semi,sol,time_index ::Int = 1;vtk ::Bool=false,out ::T="./datas") where {T<:AbstractString,F1<:Function,F2<:Function,F3<:Function}
-    ∧(v,w) = SA[v[2]*w[3]-v[3]*w[2],v[3]*w[1]-v[1]*w[3],v[1]*w[2]-v[2]*w[1]]
-    er(theta,s) = cos(theta).*nor(s) .+ sin(theta).*∧(tanj(s),nor(s))
-    return get3DData(eq,curve,er,semi,sol,time_index;vtk=vtk,out=out)
+function get3DData(
+    eq::BloodFlowEquations2D,
+    curve::F1,
+    tanj::F2,
+    nor::F3,
+    semi,
+    sol,
+    time_index::Int=1;
+    vtk::Bool=false,
+    out::T="./datas",
+) where {T<:AbstractString,F1<:Function,F2<:Function,F3<:Function}
+    ∧(v, w) = SA[
+        v[2] * w[3] - v[3] * w[2], v[3] * w[1] - v[1] * w[3], v[1] * w[2] - v[2] * w[1]
+    ]
+    er(theta, s) = cos(theta) .* nor(s) .+ sin(theta) .* ∧(tanj(s), nor(s))
+    return get3DData(eq, curve, er, semi, sol, time_index; vtk=vtk, out=out)
 end
 
 @doc raw"""
@@ -141,15 +168,24 @@ Named tuple containing:
 - `ws`: Flow axial velocity at each point.
 
 """
-function get3DData(eq::BloodFlowEquations2D,semi,sol,time_index ::Int = 1;vtk ::Bool=false,out ::T="./datas") where {T<:AbstractString}
-    curve(s) = SA[s,0.0,0.0]
-    tanj(s) = SA[1.0,0.0,0.0]
-    nor(s) = SA[0.0,1.0,0.0]
-    ∧(v,w) = SA[v[2]*w[3]-v[3]*w[2],v[3]*w[1]-v[1]*w[3],v[1]*w[2]-v[2]*w[1]]
-    er(theta,s) = cos(theta).*nor(s) .+ sin(theta).*∧(tanj(s),nor(s))
-    return get3DData(eq,curve,er,semi,sol,time_index;vtk=vtk,out=out)
+function get3DData(
+    eq::BloodFlowEquations2D,
+    semi,
+    sol,
+    time_index::Int=1;
+    vtk::Bool=false,
+    out::T="./datas",
+) where {T<:AbstractString}
+    curve(s) = SA[s, 0.0, 0.0]
+    tanj(s) = SA[1.0, 0.0, 0.0]
+    nor(s) = SA[0.0, 1.0, 0.0]
+    ∧(v, w) = SA[
+        v[2] * w[3] - v[3] * w[2], v[3] * w[1] - v[1] * w[3], v[1] * w[2] - v[2] * w[1]
+    ]
+    er(theta, s) = cos(theta) .* nor(s) .+ sin(theta) .* ∧(tanj(s), nor(s))
+    return get3DData(eq, curve, er, semi, sol, time_index; vtk=vtk, out=out)
 end
 
 function interpolate_curve(xyz_data)
     @error "add DataInterpolations.jl to use this function"
-end
\ No newline at end of file
+end
diff --git a/src/BloodFlowTrixi.jl b/src/BloodFlowTrixi.jl
index 68d6e47..409859f 100644
--- a/src/BloodFlowTrixi.jl
+++ b/src/BloodFlowTrixi.jl
@@ -12,11 +12,26 @@ Docs under https://yolhan83.github.io/BloodFlowTrixi.jl
 $(isnothing(get(ENV, "CI", nothing)) ? ("\n" * "Package local path: " * pathof(BloodFlowTrixi)) : "") 
 """
 module BloodFlowTrixi
-    using Trixi,WriteVTK,StaticArrays
-    # Write your package code here.
-    abstract type AbstractBloodFlowEquations{NDIMS, NVARS} <:Trixi.AbstractEquations{NDIMS, NVARS} end
-    include("1DModel/1dmodel.jl")
-    include("2DModel/2dmodel.jl")
+using Trixi, WriteVTK, StaticArrays
+# Write your package code here.
+abstract type AbstractBloodFlowEquations{NDIMS,NVARS} <:
+              Trixi.AbstractEquations{NDIMS,NVARS} end
+include("1DModel/1dmodel.jl")
+include("2DModel/2dmodel.jl")
 
-    export BloodFlowEquations1D,BloodFlowEquations1DOrd2, BloodFlowEquations2D,flux_nonconservative,source_term_simple,source_term_simple_ord2,boundary_condition_pressure_in,initial_condition_simple,friction,pressure,radius,boundary_condition_outflow,boundary_condition_slip_wall,get3DData,interpolate_curve
+export BloodFlowEquations1D,
+    BloodFlowEquations1DOrd2,
+    BloodFlowEquations2D,
+    flux_nonconservative,
+    source_term_simple,
+    source_term_simple_ord2,
+    boundary_condition_pressure_in,
+    initial_condition_simple,
+    friction,
+    pressure,
+    radius,
+    boundary_condition_outflow,
+    boundary_condition_slip_wall,
+    get3DData,
+    interpolate_curve
 end
diff --git a/test/Aqua/aquatest.jl b/test/Aqua/aquatest.jl
index 85bde45..0fba377 100644
--- a/test/Aqua/aquatest.jl
+++ b/test/Aqua/aquatest.jl
@@ -1,5 +1,5 @@
 using Aqua
 
 @testset "Code quality (Aqua.jl)" begin
-    Aqua.test_all(BloodFlowTrixi; ambiguities = false,)
+    Aqua.test_all(BloodFlowTrixi; ambiguities=false)
 end
diff --git a/test/Extensions/DataInterpolationsTest.jl b/test/Extensions/DataInterpolationsTest.jl
index 2ee01e9..3fec52e 100644
--- a/test/Extensions/DataInterpolationsTest.jl
+++ b/test/Extensions/DataInterpolationsTest.jl
@@ -1,6 +1,5 @@
 using DataInterpolations
 
-
 @testset "2D Blood Flow Model with interpolation" begin
     include("../../exemples/Model2D/diexemple.jl")
-end
\ No newline at end of file
+end
diff --git a/test/runtests.jl b/test/runtests.jl
index b993357..36a63a7 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -1,7 +1,7 @@
 using BloodFlowTrixi
 using Test
 
-@testset "BloodFlowTrixi.jl"  begin
+@testset "BloodFlowTrixi.jl" begin
     include("./Aqua/aquatest.jl")
     include("./Extensions/DataInterpolationsTest.jl")
 
@@ -11,4 +11,4 @@ using Test
     @testset "2D Blood Flow Model" begin
         include("../exemples/Model2D/exemple.jl")
     end
-end
\ No newline at end of file
+end