diff --git a/julia/TUG/src/AbstractSimulation.jl b/julia/TUG/src/AbstractSimulation.jl
index 05b8e6b..b3f8c66 100644
--- a/julia/TUG/src/AbstractSimulation.jl
+++ b/julia/TUG/src/AbstractSimulation.jl
@@ -32,7 +32,13 @@ Enum representing the level of CSV file output for the simulation.
 
 abstract type AbstractSimulation{T} end
 
-function adjustTimestep(grid::Grid{T}, approach::APPROACH, timestep::T, iterations::Int, verbose::Bool)::Tuple{T,Int} where {T}
+function adjustTimestep(
+    grid::Grid{T},
+    approach::APPROACH,
+    timestep::T,
+    iterations::Int,
+    verbose::Bool,
+)::Tuple{T,Int} where {T}
     if approach == FTCS
         if getDim(grid) == 1
             deltaSquare = getDeltaCol(grid)
@@ -55,7 +61,13 @@ function adjustTimestep(grid::Grid{T}, approach::APPROACH, timestep::T, iteratio
             iterations = iterations * innerIterations
             timestep = timestep / innerIterations
             if verbose
-                println("Warning: Timestep is too large for FTCS approach. Adjusting timestep to ", timestep, " and iterations to ", iterations, ".")
+                println(
+                    "Warning: Timestep is too large for FTCS approach. Adjusting timestep to ",
+                    timestep,
+                    " and iterations to ",
+                    iterations,
+                    ".",
+                )
             end
         else
             timestep = timestep
@@ -64,7 +76,10 @@ function adjustTimestep(grid::Grid{T}, approach::APPROACH, timestep::T, iteratio
     return timestep, iterations
 end
 
-function createCSVfile(simulation::AbstractSimulation{T}, grid::Union{Grid{T},Nothing}=nothing) where {T}
+function createCSVfile(
+    simulation::AbstractSimulation{T},
+    grid::Union{Grid{T},Nothing} = nothing,
+) where {T}
     grid = grid === nothing ? simulation.grid : grid
     approachString = string(simulation.approach)
     rows = getRows(grid)
@@ -76,7 +91,18 @@ function createCSVfile(simulation::AbstractSimulation{T}, grid::Union{Grid{T},No
     appendIdent = 0
     while isfile(filename)
         appendIdent += 1
-        filename = string(approachString, "_", rows, "_", cols, "_", numIterations, "-", appendIdent, ".csv")
+        filename = string(
+            approachString,
+            "_",
+            rows,
+            "_",
+            cols,
+            "_",
+            numIterations,
+            "-",
+            appendIdent,
+            ".csv",
+        )
     end
 
     # Write boundary conditions if required
@@ -103,12 +129,19 @@ function writeBoundarySideValuesToCSV(file::IOStream, bc::Boundary{T}, side) whe
     write(file, formatted_values, "\n")
 end
 
-function writeConcentrationsToCLI(simulation::AbstractSimulation{T}, grid::Union{Grid{T},Nothing}=nothing) where {T}
+function writeConcentrationsToCLI(
+    simulation::AbstractSimulation{T},
+    grid::Union{Grid{T},Nothing} = nothing,
+) where {T}
     grid = grid === nothing ? simulation.grid : grid
     println(getConcentrations(grid))
 end
 
-function writeConcentrationsToCSV(file::IOStream, simulation::AbstractSimulation{T}, grid::Union{Grid{T},Nothing}=nothing) where {T}
+function writeConcentrationsToCSV(
+    file::IOStream,
+    simulation::AbstractSimulation{T},
+    grid::Union{Grid{T},Nothing} = nothing,
+) where {T}
     grid = grid === nothing ? simulation.grid : grid
     concentrations = getConcentrations(grid)
 
diff --git a/julia/TUG/src/Boundary.jl b/julia/TUG/src/Boundary.jl
index 716784e..b3c3909 100644
--- a/julia/TUG/src/Boundary.jl
+++ b/julia/TUG/src/Boundary.jl
@@ -128,10 +128,10 @@ struct Boundary{T}
             boundaries[Int(LEFT)] = [BoundaryElement{T}()]
             boundaries[Int(RIGHT)] = [BoundaryElement{T}()]
         elseif dim == 2
-            boundaries[Int(LEFT)] = [BoundaryElement{T}() for _ in 1:rows]
-            boundaries[Int(RIGHT)] = [BoundaryElement{T}() for _ in 1:rows]
-            boundaries[Int(TOP)] = [BoundaryElement{T}() for _ in 1:cols]
-            boundaries[Int(BOTTOM)] = [BoundaryElement{T}() for _ in 1:cols]
+            boundaries[Int(LEFT)] = [BoundaryElement{T}() for _ = 1:rows]
+            boundaries[Int(RIGHT)] = [BoundaryElement{T}() for _ = 1:rows]
+            boundaries[Int(TOP)] = [BoundaryElement{T}() for _ = 1:cols]
+            boundaries[Int(BOTTOM)] = [BoundaryElement{T}() for _ = 1:cols]
         else
             throw(ArgumentError("Only 1- and 2-dimensional grids are defined!"))
         end
@@ -153,9 +153,17 @@ Retrieves the type of the boundary element at the specified side and index.
 # Returns
 The type of the boundary element (`TYPE`).
 """
-function getBoundaryElementType(boundary::Boundary{T}, side::SIDE, index::Int)::TYPE where {T}
+function getBoundaryElementType(
+    boundary::Boundary{T},
+    side::SIDE,
+    index::Int,
+)::TYPE where {T}
     if boundary.dim == 1 && (side == BOTTOM || side == TOP)
-        throw(ArgumentError("For the one-dimensional case, only the left and right borders exist."))
+        throw(
+            ArgumentError(
+                "For the one-dimensional case, only the left and right borders exist.",
+            ),
+        )
     end
     if index < 1 || index > (side in (LEFT, RIGHT) ? boundary.rows : boundary.cols)
         throw(ArgumentError("Index out of bounds!"))
@@ -180,7 +188,11 @@ The value of the boundary element.
 """
 function getBoundaryElementValue(boundary::Boundary{T}, side::SIDE, index::Int)::T where {T}
     if boundary.dim == 1 && (side == BOTTOM || side == TOP)
-        throw(ArgumentError("For the one-dimensional case, only the left and right borders exist."))
+        throw(
+            ArgumentError(
+                "For the one-dimensional case, only the left and right borders exist.",
+            ),
+        )
     end
     if index < 1 || index > (side in (LEFT, RIGHT) ? boundary.rows : boundary.cols)
         throw(ArgumentError("Index out of bounds!"))
@@ -201,9 +213,16 @@ Retrieves the boundary elements at the specified side.
 # Returns
 The boundary elements at the specified side.
 """
-function getBoundarySide(boundary::Boundary{T}, side::SIDE)::Vector{BoundaryElement{T}} where {T}
+function getBoundarySide(
+    boundary::Boundary{T},
+    side::SIDE,
+)::Vector{BoundaryElement{T}} where {T}
     if boundary.dim == 1 && (side == BOTTOM || side == TOP)
-        throw(ArgumentError("For the one-dimensional case, only the left and right borders exist."))
+        throw(
+            ArgumentError(
+                "For the one-dimensional case, only the left and right borders exist.",
+            ),
+        )
     end
 
     boundary.boundaries[Int(side)]
@@ -221,15 +240,22 @@ Sets the boundary elements at the specified side to CLOSED.
 # Returns
 The boundary elements at the specified side.
 """
-function setBoundarySideClosed!(boundary::Boundary{T}, side::SIDE)::Vector{BoundaryElement{T}} where {T}
+function setBoundarySideClosed!(
+    boundary::Boundary{T},
+    side::SIDE,
+)::Vector{BoundaryElement{T}} where {T}
     if boundary.dim == 1 && (side == BOTTOM || side == TOP)
-        throw(ArgumentError("For the one-dimensional case, only the left and right borders exist."))
+        throw(
+            ArgumentError(
+                "For the one-dimensional case, only the left and right borders exist.",
+            ),
+        )
     end
 
     is_vertical = side in (LEFT, RIGHT)
     n = is_vertical ? boundary.rows : boundary.cols
 
-    boundary.boundaries[Int(side)] = [BoundaryElement{T}() for _ in 1:n]
+    boundary.boundaries[Int(side)] = [BoundaryElement{T}() for _ = 1:n]
 end
 
 """
@@ -245,13 +271,21 @@ Sets the boundary elements at the specified side to CONSTANT with the specified
 # Returns
 The boundary elements at the specified side.
 """
-function setBoundarySideConstant!(boundary::Boundary{T}, side::SIDE, value::T)::Vector{BoundaryElement{T}} where {T}
+function setBoundarySideConstant!(
+    boundary::Boundary{T},
+    side::SIDE,
+    value::T,
+)::Vector{BoundaryElement{T}} where {T}
     if boundary.dim == 1 && (side == BOTTOM || side == TOP)
-        throw(ArgumentError("For the one-dimensional case, only the left and right borders exist."))
+        throw(
+            ArgumentError(
+                "For the one-dimensional case, only the left and right borders exist.",
+            ),
+        )
     end
 
     is_vertical = side in (LEFT, RIGHT)
     n = is_vertical ? boundary.rows : boundary.cols
 
-    boundary.boundaries[Int(side)] = [BoundaryElement{T}(value) for _ in 1:n]
+    boundary.boundaries[Int(side)] = [BoundaryElement{T}(value) for _ = 1:n]
 end
diff --git a/julia/TUG/src/Core/BTCS.jl b/julia/TUG/src/Core/BTCS.jl
index 26a058a..70e31c3 100644
--- a/julia/TUG/src/Core/BTCS.jl
+++ b/julia/TUG/src/Core/BTCS.jl
@@ -8,7 +8,12 @@ using Base.Threads
 using LinearAlgebra
 using SparseArrays
 
-function calcBoundaryCoeff(alpha_center::T, alpha_side::T, sx::T, bcType::TYPE)::Tuple{T,T} where {T}
+function calcBoundaryCoeff(
+    alpha_center::T,
+    alpha_side::T,
+    sx::T,
+    bcType::TYPE,
+)::Tuple{T,T} where {T}
     alpha = calcAlphaIntercell(alpha_center, alpha_side)
     sideCoeff = -sx * alpha
 
@@ -23,12 +28,31 @@ function calcBoundaryCoeff(alpha_center::T, alpha_side::T, sx::T, bcType::TYPE):
     return (centerCoeff, sideCoeff)
 end
 
-function createCoeffMatrix(alpha::Matrix{T}, alpha_left::Vector{T}, alpha_right::Vector{T}, bcLeft::Vector{BoundaryElement{T}}, bcRight::Vector{BoundaryElement{T}}, numCols::Int, rowIndex::Int, sx::T)::Tridiagonal{T} where {T}
+function createCoeffMatrix(
+    alpha::Matrix{T},
+    alpha_left::Vector{T},
+    alpha_right::Vector{T},
+    bcLeft::Vector{BoundaryElement{T}},
+    bcRight::Vector{BoundaryElement{T}},
+    numCols::Int,
+    rowIndex::Int,
+    sx::T,
+)::Tridiagonal{T} where {T}
     # Determine left side boundary coefficients based on boundary condition
-    centerCoeffTop, rightCoeffTop = calcBoundaryCoeff(alpha[rowIndex, 1], alpha[rowIndex, 2], sx, getType(bcLeft[rowIndex]))
+    centerCoeffTop, rightCoeffTop = calcBoundaryCoeff(
+        alpha[rowIndex, 1],
+        alpha[rowIndex, 2],
+        sx,
+        getType(bcLeft[rowIndex]),
+    )
 
     # Determine right side boundary coefficients based on boundary condition
-    centerCoeffBottom, leftCoeffBottom = calcBoundaryCoeff(alpha[rowIndex, numCols], alpha[rowIndex, numCols-1], sx, getType(bcRight[rowIndex]))
+    centerCoeffBottom, leftCoeffBottom = calcBoundaryCoeff(
+        alpha[rowIndex, numCols],
+        alpha[rowIndex, numCols-1],
+        sx,
+        getType(bcRight[rowIndex]),
+    )
 
     dl = [-sx .* alpha_left; leftCoeffBottom]
     d = [centerCoeffTop; 1 .+ sx .* (alpha_right + alpha_left); centerCoeffBottom]
@@ -37,26 +61,55 @@ function createCoeffMatrix(alpha::Matrix{T}, alpha_left::Vector{T}, alpha_right:
 end
 
 
-function calcExplicitConcentrationsBoundaryClosed(conc_center::T, alpha_center::T, alpha_neighbor::T, sy::T)::T where {T}
+function calcExplicitConcentrationsBoundaryClosed(
+    conc_center::T,
+    alpha_center::T,
+    alpha_neighbor::T,
+    sy::T,
+)::T where {T}
     alpha = calcAlphaIntercell(alpha_center, alpha_neighbor)
     return (sy * alpha + (1 - sy * alpha)) * conc_center
 end
 
-function calcExplicitConcentrationsBoundaryConstant(conc_center::T, conc_bc::T, alpha_center::T, alpha_neighbor::T, sy::T)::T where {T}
+function calcExplicitConcentrationsBoundaryConstant(
+    conc_center::T,
+    conc_bc::T,
+    alpha_center::T,
+    alpha_neighbor::T,
+    sy::T,
+)::T where {T}
     alpha_center_neighbor = calcAlphaIntercell(alpha_center, alpha_neighbor)
-    alpha_center_center = alpha_center == alpha_neighbor ? alpha_center_neighbor : calcAlphaIntercell(alpha_center, alpha_center)
+    alpha_center_center =
+        alpha_center == alpha_neighbor ? alpha_center_neighbor :
+        calcAlphaIntercell(alpha_center, alpha_center)
 
     return sy * alpha_center_neighbor * conc_center +
            (1 - sy * (alpha_center_center + 2 * alpha_center)) * conc_center +
            sy * alpha_center * conc_bc
 end
 
-function calcExplicitConcentrationsBoundaryClosed(conc_center::Vector{T}, alpha::Vector{T}, sy::T)::T where {T}
+function calcExplicitConcentrationsBoundaryClosed(
+    conc_center::Vector{T},
+    alpha::Vector{T},
+    sy::T,
+)::T where {T}
     return (sy .* alpha .+ (1 .- sy .* alpha)) .* conc_center
 end
 
 # creates a solution vector for next time step from the current state of concentrations
-function writeSolutionVector!(sv::Vector{T}, concentrations::Matrix{T}, alphaX::Matrix{T}, alphaY::Matrix{T}, bcLeft::Vector{BoundaryElement{T}}, bcRight::Vector{BoundaryElement{T}}, bcTop::Vector{BoundaryElement{T}}, bcBottom::Vector{BoundaryElement{T}}, rowIndex::Int, sx::T, sy::T) where {T}
+function writeSolutionVector!(
+    sv::Vector{T},
+    concentrations::Matrix{T},
+    alphaX::Matrix{T},
+    alphaY::Matrix{T},
+    bcLeft::Vector{BoundaryElement{T}},
+    bcRight::Vector{BoundaryElement{T}},
+    bcTop::Vector{BoundaryElement{T}},
+    bcBottom::Vector{BoundaryElement{T}},
+    rowIndex::Int,
+    sx::T,
+    sy::T,
+) where {T}
     numRows = size(concentrations, 1)
     length = size(sv, 1)
 
@@ -64,9 +117,20 @@ function writeSolutionVector!(sv::Vector{T}, concentrations::Matrix{T}, alphaX::
     if rowIndex == 1
         @inbounds for i = 1:length
             if getType(bcTop[i]) == CONSTANT
-                sv[i] = calcExplicitConcentrationsBoundaryConstant(concentrations[rowIndex, i], getValue(bcTop[i]), alphaY[rowIndex, i], alphaY[rowIndex+1, i], sy)
+                sv[i] = calcExplicitConcentrationsBoundaryConstant(
+                    concentrations[rowIndex, i],
+                    getValue(bcTop[i]),
+                    alphaY[rowIndex, i],
+                    alphaY[rowIndex+1, i],
+                    sy,
+                )
             elseif getType(bcTop[i]) == CLOSED
-                sv[i] = calcExplicitConcentrationsBoundaryClosed(concentrations[rowIndex, i], alphaY[rowIndex, i], alphaY[rowIndex+1, i], sy)
+                sv[i] = calcExplicitConcentrationsBoundaryClosed(
+                    concentrations[rowIndex, i],
+                    alphaY[rowIndex, i],
+                    alphaY[rowIndex+1, i],
+                    sy,
+                )
             else
                 error("Undefined Boundary Condition Type somewhere on Left or Top!")
             end
@@ -76,11 +140,16 @@ function writeSolutionVector!(sv::Vector{T}, concentrations::Matrix{T}, alphaX::
     # Inner rows
     if rowIndex > 1 && rowIndex < numRows
         @inbounds for i = 1:length
-            alpha_here_below = calcAlphaIntercell(alphaY[rowIndex, i], alphaY[rowIndex+1, i])
-            alpha_here_above = alphaY[rowIndex+1, i] == alphaY[rowIndex-1, i] ? alpha_here_below : calcAlphaIntercell(alphaY[rowIndex-1, i], alphaY[rowIndex, i]) # calcAlphaIntercell is symmetric, so we can use it for both directions
-            sv[i] = sy * alpha_here_below * concentrations[rowIndex+1, i] +
-                    (1 - sy * (alpha_here_below + alpha_here_above)) * concentrations[rowIndex, i] +
-                    sy * alpha_here_above * concentrations[rowIndex-1, i]
+            alpha_here_below =
+                calcAlphaIntercell(alphaY[rowIndex, i], alphaY[rowIndex+1, i])
+            alpha_here_above =
+                alphaY[rowIndex+1, i] == alphaY[rowIndex-1, i] ? alpha_here_below :
+                calcAlphaIntercell(alphaY[rowIndex-1, i], alphaY[rowIndex, i]) # calcAlphaIntercell is symmetric, so we can use it for both directions
+            sv[i] =
+                sy * alpha_here_below * concentrations[rowIndex+1, i] +
+                (1 - sy * (alpha_here_below + alpha_here_above)) *
+                concentrations[rowIndex, i] +
+                sy * alpha_here_above * concentrations[rowIndex-1, i]
         end
     end
 
@@ -88,9 +157,20 @@ function writeSolutionVector!(sv::Vector{T}, concentrations::Matrix{T}, alphaX::
     if rowIndex == numRows
         @inbounds for i = 1:length
             if getType(bcBottom[i]) == CONSTANT
-                sv[i] = calcExplicitConcentrationsBoundaryConstant(concentrations[rowIndex, i], getValue(bcBottom[i]), alphaY[rowIndex, i], alphaY[rowIndex-1, i], sy)
+                sv[i] = calcExplicitConcentrationsBoundaryConstant(
+                    concentrations[rowIndex, i],
+                    getValue(bcBottom[i]),
+                    alphaY[rowIndex, i],
+                    alphaY[rowIndex-1, i],
+                    sy,
+                )
             elseif getType(bcBottom[i]) == CLOSED
-                sv[i] = calcExplicitConcentrationsBoundaryClosed(concentrations[rowIndex, i], alphaY[rowIndex, i], alphaY[rowIndex-1, i], sy)
+                sv[i] = calcExplicitConcentrationsBoundaryClosed(
+                    concentrations[rowIndex, i],
+                    alphaY[rowIndex, i],
+                    alphaY[rowIndex-1, i],
+                    sy,
+                )
             else
                 error("Undefined Boundary Condition Type somewhere on Right or Bottom!")
             end
@@ -108,7 +188,13 @@ function writeSolutionVector!(sv::Vector{T}, concentrations::Matrix{T}, alphaX::
     end
 end
 
-function BTCS_1D(grid::Grid{T}, bc::Boundary{T}, alpha_left::Matrix{T}, alpha_right::Matrix{T}, timestep::T) where {T}
+function BTCS_1D(
+    grid::Grid{T},
+    bc::Boundary{T},
+    alpha_left::Matrix{T},
+    alpha_right::Matrix{T},
+    timestep::T,
+) where {T}
     sx = timestep / (getDeltaCol(grid) * getDeltaCol(grid))
 
     alpha = getAlphaX(grid)
@@ -118,7 +204,16 @@ function BTCS_1D(grid::Grid{T}, bc::Boundary{T}, alpha_left::Matrix{T}, alpha_ri
 
     concentrations::Matrix{T} = getConcentrations(grid)
 
-    A::Tridiagonal{T} = createCoeffMatrix(alpha, alpha_left[1, :], alpha_right[1, :], bcLeft, bcRight, length, 1, sx)
+    A::Tridiagonal{T} = createCoeffMatrix(
+        alpha,
+        alpha_left[1, :],
+        alpha_right[1, :],
+        bcLeft,
+        bcRight,
+        length,
+        1,
+        sx,
+    )
     b = concentrations[1, :]
 
     if getType(bcLeft[1]) == CONSTANT
@@ -131,7 +226,15 @@ function BTCS_1D(grid::Grid{T}, bc::Boundary{T}, alpha_left::Matrix{T}, alpha_ri
     concentrations[1, :] = A \ b
 end
 
-function BTCS_2D(grid::Grid{T}, bc::Boundary{T}, alphaX_left::Matrix{T}, alphaX_right::Matrix{T}, alphaY_t_left::Matrix{T}, alphaY_t_right::Matrix{T}, timestep::T) where {T}
+function BTCS_2D(
+    grid::Grid{T},
+    bc::Boundary{T},
+    alphaX_left::Matrix{T},
+    alphaX_right::Matrix{T},
+    alphaY_t_left::Matrix{T},
+    alphaY_t_right::Matrix{T},
+    timestep::T,
+) where {T}
     rows = getRows(grid)
     cols = getCols(grid)
     sx = timestep / (2 * getDeltaCol(grid) * getDeltaCol(grid))
@@ -152,11 +255,32 @@ function BTCS_2D(grid::Grid{T}, bc::Boundary{T}, alphaX_left::Matrix{T}, alphaX_
     bcTop = getBoundarySide(bc, TOP)
     bcBottom = getBoundarySide(bc, BOTTOM)
 
-    localBs = [zeros(T, cols) for _ in 1:Threads.nthreads()]
+    localBs = [zeros(T, cols) for _ = 1:Threads.nthreads()]
     Threads.@threads for i = 1:rows
         localB = localBs[Threads.threadid()]
-        A::Tridiagonal{T} = createCoeffMatrix(alphaX, alphaX_left[i, :], alphaX_right[i, :], bcLeft, bcRight, cols, i, sx)
-        writeSolutionVector!(localB, concentrations, alphaX, alphaY, bcLeft, bcRight, bcTop, bcBottom, i, sx, sy)
+        A::Tridiagonal{T} = createCoeffMatrix(
+            alphaX,
+            alphaX_left[i, :],
+            alphaX_right[i, :],
+            bcLeft,
+            bcRight,
+            cols,
+            i,
+            sx,
+        )
+        writeSolutionVector!(
+            localB,
+            concentrations,
+            alphaX,
+            alphaY,
+            bcLeft,
+            bcRight,
+            bcTop,
+            bcBottom,
+            i,
+            sx,
+            sy,
+        )
 
         concentrations_intermediate[i, :] = A \ localB
     end
@@ -165,11 +289,32 @@ function BTCS_2D(grid::Grid{T}, bc::Boundary{T}, alphaX_left::Matrix{T}, alphaX_
     concentrations_t = fetch(concentrations_t_task)
 
     # Swap alphas, boundary conditions and sx/sy for column-wise calculation
-    localBs = [zeros(T, rows) for _ in 1:Threads.nthreads()]
+    localBs = [zeros(T, rows) for _ = 1:Threads.nthreads()]
     Threads.@threads for i = 1:cols
         localB = localBs[Threads.threadid()]
-        A::Tridiagonal{T} = createCoeffMatrix(alphaY_t, alphaY_t_left[i, :], alphaY_t_right[i, :], bcTop, bcBottom, rows, i, sy)
-        writeSolutionVector!(localB, concentrations_intermediate, alphaY_t, alphaX_t, bcTop, bcBottom, bcLeft, bcRight, i, sy, sx)
+        A::Tridiagonal{T} = createCoeffMatrix(
+            alphaY_t,
+            alphaY_t_left[i, :],
+            alphaY_t_right[i, :],
+            bcTop,
+            bcBottom,
+            rows,
+            i,
+            sy,
+        )
+        writeSolutionVector!(
+            localB,
+            concentrations_intermediate,
+            alphaY_t,
+            alphaX_t,
+            bcTop,
+            bcBottom,
+            bcLeft,
+            bcRight,
+            i,
+            sy,
+            sx,
+        )
 
         concentrations_t[i, :] = A \ localB
     end
@@ -179,18 +324,28 @@ function BTCS_2D(grid::Grid{T}, bc::Boundary{T}, alphaX_left::Matrix{T}, alphaX_
     setConcentrations!(grid, concentrations)
 end
 
-function runBTCS(grid::Grid{T}, bc::Boundary{T}, timestep::T, iterations::Int, stepCallback::Function) where {T}
+function runBTCS(
+    grid::Grid{T},
+    bc::Boundary{T},
+    timestep::T,
+    iterations::Int,
+    stepCallback::Function,
+) where {T}
     if getDim(grid) == 1
         length = getCols(grid)
         alpha = getAlphaX(grid)
 
-        alpha_left_task = Threads.@spawn calcAlphaIntercell(alpha[:, 1:(length-2)], alpha[:, 2:(length-1)])
-        alpha_right_task = Threads.@spawn calcAlphaIntercell(alpha[:, 2:(length-1)], alpha[:, 3:length])
+        alpha_left_task = Threads.@spawn calcAlphaIntercell(
+            alpha[:, 1:(length-2)],
+            alpha[:, 2:(length-1)],
+        )
+        alpha_right_task =
+            Threads.@spawn calcAlphaIntercell(alpha[:, 2:(length-1)], alpha[:, 3:length])
 
         alpha_left = fetch(alpha_left_task)
         alpha_right = fetch(alpha_right_task)
 
-        for _ in 1:(iterations)
+        for _ = 1:(iterations)
             BTCS_1D(grid, bc, alpha_left, alpha_right, timestep)
 
             stepCallback()
@@ -201,18 +356,32 @@ function runBTCS(grid::Grid{T}, bc::Boundary{T}, timestep::T, iterations::Int, s
         alphaX = getAlphaX(grid)
         alphaY_t = getAlphaY_t(grid)
 
-        alphaX_left_task = Threads.@spawn calcAlphaIntercell(alphaX[:, 1:(cols-2)], alphaX[:, 2:(cols-1)])
-        alphaX_right_task = Threads.@spawn calcAlphaIntercell(alphaX[:, 2:(cols-1)], alphaX[:, 3:cols])
-        alphaY_t_left_task = Threads.@spawn calcAlphaIntercell(alphaY_t[:, 1:(rows-2)], alphaY_t[:, 2:(rows-1)])
-        alphaY_t_right_task = Threads.@spawn calcAlphaIntercell(alphaY_t[:, 2:(rows-1)], alphaY_t[:, 3:rows])
+        alphaX_left_task =
+            Threads.@spawn calcAlphaIntercell(alphaX[:, 1:(cols-2)], alphaX[:, 2:(cols-1)])
+        alphaX_right_task =
+            Threads.@spawn calcAlphaIntercell(alphaX[:, 2:(cols-1)], alphaX[:, 3:cols])
+        alphaY_t_left_task = Threads.@spawn calcAlphaIntercell(
+            alphaY_t[:, 1:(rows-2)],
+            alphaY_t[:, 2:(rows-1)],
+        )
+        alphaY_t_right_task =
+            Threads.@spawn calcAlphaIntercell(alphaY_t[:, 2:(rows-1)], alphaY_t[:, 3:rows])
 
         alphaX_left = fetch(alphaX_left_task)
         alphaX_right = fetch(alphaX_right_task)
         alphaY_t_left = fetch(alphaY_t_left_task)
         alphaY_t_right = fetch(alphaY_t_right_task)
 
-        for _ in 1:(iterations)
-            BTCS_2D(grid, bc, alphaX_left, alphaX_right, alphaY_t_left, alphaY_t_right, timestep)
+        for _ = 1:(iterations)
+            BTCS_2D(
+                grid,
+                bc,
+                alphaX_left,
+                alphaX_right,
+                alphaY_t_left,
+                alphaY_t_right,
+                timestep,
+            )
 
             stepCallback()
         end
diff --git a/julia/TUG/src/Core/FTCS.jl b/julia/TUG/src/Core/FTCS.jl
index 26aede4..c60fbc4 100644
--- a/julia/TUG/src/Core/FTCS.jl
+++ b/julia/TUG/src/Core/FTCS.jl
@@ -7,8 +7,14 @@ using Base.Threads
 using LinearAlgebra
 using SparseArrays
 
-function calcHorizontalChange(alphaX_next::T, alphaX_prev::T, alphaX_current::T,
-    concentration_next::T, concentration_prev::T, concentration_current::T) where {T}
+function calcHorizontalChange(
+    alphaX_next::T,
+    alphaX_prev::T,
+    alphaX_current::T,
+    concentration_next::T,
+    concentration_prev::T,
+    concentration_current::T,
+) where {T}
 
     intercellAlpha_next = calcAlphaIntercell(alphaX_next, alphaX_current)
     intercellAlpha_prev = calcAlphaIntercell(alphaX_prev, alphaX_current)
@@ -18,8 +24,14 @@ function calcHorizontalChange(alphaX_next::T, alphaX_prev::T, alphaX_current::T,
            intercellAlpha_prev * concentration_prev
 end
 
-function calcHorizontalChangeLeftBoundary(boundaryType::TYPE, alphaX_next::T, alphaX_current::T,
-    concentration_next::T, concentration_current::T, boundaryValue::T) where {T}
+function calcHorizontalChangeLeftBoundary(
+    boundaryType::TYPE,
+    alphaX_next::T,
+    alphaX_current::T,
+    concentration_next::T,
+    concentration_current::T,
+    boundaryValue::T,
+) where {T}
     intercellAlpha = calcAlphaIntercell(alphaX_next, alphaX_current)
 
     if boundaryType == CONSTANT
@@ -33,8 +45,14 @@ function calcHorizontalChangeLeftBoundary(boundaryType::TYPE, alphaX_next::T, al
     end
 end
 
-function calcHorizontalChangeRightBoundary(boundaryType::TYPE, alphaX_prev::T, alphaX_current::T,
-    concentration_prev::T, concentration_current::T, boundaryValue::T) where {T}
+function calcHorizontalChangeRightBoundary(
+    boundaryType::TYPE,
+    alphaX_prev::T,
+    alphaX_current::T,
+    concentration_prev::T,
+    concentration_current::T,
+    boundaryValue::T,
+) where {T}
     intercellAlpha = calcAlphaIntercell(alphaX_prev, alphaX_current)
 
     if boundaryType == CONSTANT
@@ -48,8 +66,14 @@ function calcHorizontalChangeRightBoundary(boundaryType::TYPE, alphaX_prev::T, a
     end
 end
 
-function calcVerticalChange(alphaY_above::T, alphaY_below::T, alphaY_current::T,
-    concentration_above::T, concentration_below::T, concentration_current::T) where {T}
+function calcVerticalChange(
+    alphaY_above::T,
+    alphaY_below::T,
+    alphaY_current::T,
+    concentration_above::T,
+    concentration_below::T,
+    concentration_current::T,
+) where {T}
     intercellAlpha_above = calcAlphaIntercell(alphaY_above, alphaY_current)
     intercellAlpha_below = calcAlphaIntercell(alphaY_below, alphaY_current)
 
@@ -58,8 +82,14 @@ function calcVerticalChange(alphaY_above::T, alphaY_below::T, alphaY_current::T,
            intercellAlpha_below * concentration_below
 end
 
-function calcVerticalChangeTopBoundary(boundaryType::TYPE, alphaY_above::T, alphaY_current::T,
-    concentration_above::T, concentration_current::T, boundaryValue::T) where {T}
+function calcVerticalChangeTopBoundary(
+    boundaryType::TYPE,
+    alphaY_above::T,
+    alphaY_current::T,
+    concentration_above::T,
+    concentration_current::T,
+    boundaryValue::T,
+) where {T}
     intercellAlpha = calcAlphaIntercell(alphaY_above, alphaY_current)
 
     if boundaryType == CONSTANT
@@ -74,8 +104,14 @@ function calcVerticalChangeTopBoundary(boundaryType::TYPE, alphaY_above::T, alph
 end
 
 
-function calcVerticalChangeBottomBoundary(boundaryType::TYPE, alphaY_current::T, alphaY_below::T,
-    concentration_current::T, concentration_below::T, boundaryValue::T) where {T}
+function calcVerticalChangeBottomBoundary(
+    boundaryType::TYPE,
+    alphaY_current::T,
+    alphaY_below::T,
+    concentration_current::T,
+    concentration_below::T,
+    boundaryValue::T,
+) where {T}
     intercellAlpha = calcAlphaIntercell(alphaY_current, alphaY_below)
 
     if boundaryType == CONSTANT
@@ -101,21 +137,42 @@ function FTCS_1D(grid::Grid{T}, bc::Boundary{T}, timestep::T) where {T}
 
     # inner cells
     for col = 2:colMax-1
-        concentrations_t1[row, col] += sx * calcHorizontalChange(alphaX[row, col+1], alphaX[row, col-1], alphaX[row, col],
-            concentrations[row, col+1], concentrations[row, col-1], concentrations[row, col])
+        concentrations_t1[row, col] +=
+            sx * calcHorizontalChange(
+                alphaX[row, col+1],
+                alphaX[row, col-1],
+                alphaX[row, col],
+                concentrations[row, col+1],
+                concentrations[row, col-1],
+                concentrations[row, col],
+            )
     end
 
     # left boundary
     leftBoundaryType = getBoundaryElementType(bc, LEFT, row)
     leftBoundaryValue = getBoundaryElementValue(bc, LEFT, row)
-    concentrations_t1[row, 1] += sx * calcHorizontalChangeLeftBoundary(leftBoundaryType, alphaX[row, 2], alphaX[row, 1],
-        concentrations[row, 2], concentrations[row, 1], leftBoundaryValue)
+    concentrations_t1[row, 1] +=
+        sx * calcHorizontalChangeLeftBoundary(
+            leftBoundaryType,
+            alphaX[row, 2],
+            alphaX[row, 1],
+            concentrations[row, 2],
+            concentrations[row, 1],
+            leftBoundaryValue,
+        )
 
     # right boundary
     rightBoundaryType = getBoundaryElementType(bc, RIGHT, row)
     rightBoundaryValue = getBoundaryElementValue(bc, RIGHT, row)
-    concentrations_t1[row, colMax] += sx * calcHorizontalChangeRightBoundary(rightBoundaryType, alphaX[row, colMax-1], alphaX[row, colMax],
-        concentrations[row, colMax-1], concentrations[row, colMax], rightBoundaryValue)
+    concentrations_t1[row, colMax] +=
+        sx * calcHorizontalChangeRightBoundary(
+            rightBoundaryType,
+            alphaX[row, colMax-1],
+            alphaX[row, colMax],
+            concentrations[row, colMax-1],
+            concentrations[row, colMax],
+            rightBoundaryValue,
+        )
 
     setConcentrations!(grid, concentrations_t1)
 end
@@ -139,78 +196,201 @@ function FTCS_2D(grid::Grid{T}, bc::Boundary{T}, timestep::T) where {T}
 
     Threads.@threads for row = 2:rowMax-1
         for col = 2:colMax-1
-            concentrations_t1[row, col] += sy * calcVerticalChange(alphaY[row+1, col], alphaY[row-1, col], alphaY[row, col],
-                concentrations[row+1, col], concentrations[row-1, col], concentrations[row, col]) +
-                                           sx * calcHorizontalChange(alphaX[row, col+1], alphaX[row, col-1], alphaX[row, col],
-                concentrations[row, col+1], concentrations[row, col-1], concentrations[row, col])
+            concentrations_t1[row, col] +=
+                sy * calcVerticalChange(
+                    alphaY[row+1, col],
+                    alphaY[row-1, col],
+                    alphaY[row, col],
+                    concentrations[row+1, col],
+                    concentrations[row-1, col],
+                    concentrations[row, col],
+                ) +
+                sx * calcHorizontalChange(
+                    alphaX[row, col+1],
+                    alphaX[row, col-1],
+                    alphaX[row, col],
+                    concentrations[row, col+1],
+                    concentrations[row, col-1],
+                    concentrations[row, col],
+                )
         end
 
         # Boundary conditions for each row
         # Left boundary without corners
         leftBoundaryValue = getBoundaryElementValue(bc, LEFT, row)
-        concentrations_t1[row, 1] += sx * calcHorizontalChangeLeftBoundary(leftBoundaryType, alphaX[row, 2], alphaX[row, 1],
-            concentrations[row, 2], concentrations[row, 1], leftBoundaryValue) +
-                                     sy * calcVerticalChange(alphaY[row+1, 1], alphaY[row-1, 1], alphaY[row, 1],
-            concentrations[row+1, 1], concentrations[row-1, 1], concentrations[row, 1])
+        concentrations_t1[row, 1] +=
+            sx * calcHorizontalChangeLeftBoundary(
+                leftBoundaryType,
+                alphaX[row, 2],
+                alphaX[row, 1],
+                concentrations[row, 2],
+                concentrations[row, 1],
+                leftBoundaryValue,
+            ) +
+            sy * calcVerticalChange(
+                alphaY[row+1, 1],
+                alphaY[row-1, 1],
+                alphaY[row, 1],
+                concentrations[row+1, 1],
+                concentrations[row-1, 1],
+                concentrations[row, 1],
+            )
 
         # Right boundary without corners
         rightBoundaryValue = getBoundaryElementValue(bc, RIGHT, row)
-        concentrations_t1[row, colMax] += sx * calcHorizontalChangeRightBoundary(rightBoundaryType, alphaX[row, colMax-1], alphaX[row, colMax],
-            concentrations[row, colMax-1], concentrations[row, colMax], rightBoundaryValue) +
-                                          sy * calcVerticalChange(alphaY[row+1, colMax], alphaY[row-1, colMax], alphaY[row, colMax],
-            concentrations[row+1, colMax], concentrations[row-1, colMax], concentrations[row, colMax])
+        concentrations_t1[row, colMax] +=
+            sx * calcHorizontalChangeRightBoundary(
+                rightBoundaryType,
+                alphaX[row, colMax-1],
+                alphaX[row, colMax],
+                concentrations[row, colMax-1],
+                concentrations[row, colMax],
+                rightBoundaryValue,
+            ) +
+            sy * calcVerticalChange(
+                alphaY[row+1, colMax],
+                alphaY[row-1, colMax],
+                alphaY[row, colMax],
+                concentrations[row+1, colMax],
+                concentrations[row-1, colMax],
+                concentrations[row, colMax],
+            )
     end
 
     # Handle top/bottom boundaries
     Threads.@threads for col = 2:colMax-1
         # Top boundary
         topBoundaryValue = getBoundaryElementValue(bc, TOP, col)
-        concentrations_t1[1, col] += sy * calcVerticalChangeTopBoundary(topBoundaryType, alphaY[2, col], alphaY[1, col],
-            concentrations[2, col], concentrations[1, col], topBoundaryValue) +
-                                     sx * calcHorizontalChange(alphaX[1, col+1], alphaX[1, col-1], alphaX[1, col],
-            concentrations[1, col+1], concentrations[1, col-1], concentrations[1, col])
+        concentrations_t1[1, col] +=
+            sy * calcVerticalChangeTopBoundary(
+                topBoundaryType,
+                alphaY[2, col],
+                alphaY[1, col],
+                concentrations[2, col],
+                concentrations[1, col],
+                topBoundaryValue,
+            ) +
+            sx * calcHorizontalChange(
+                alphaX[1, col+1],
+                alphaX[1, col-1],
+                alphaX[1, col],
+                concentrations[1, col+1],
+                concentrations[1, col-1],
+                concentrations[1, col],
+            )
 
         # Bottom boundary
         bottomBoundaryValue = getBoundaryElementValue(bc, BOTTOM, col)
-        concentrations_t1[rowMax, col] += sy * calcVerticalChangeBottomBoundary(bottomBoundaryType, alphaY[rowMax, col], alphaY[rowMax-1, col],
-            concentrations[rowMax, col], concentrations[rowMax-1, col], bottomBoundaryValue) +
-                                          sx * calcHorizontalChange(alphaX[rowMax, col+1], alphaX[rowMax, col-1], alphaX[rowMax, col],
-            concentrations[rowMax, col+1], concentrations[rowMax, col-1], concentrations[rowMax, col])
+        concentrations_t1[rowMax, col] +=
+            sy * calcVerticalChangeBottomBoundary(
+                bottomBoundaryType,
+                alphaY[rowMax, col],
+                alphaY[rowMax-1, col],
+                concentrations[rowMax, col],
+                concentrations[rowMax-1, col],
+                bottomBoundaryValue,
+            ) +
+            sx * calcHorizontalChange(
+                alphaX[rowMax, col+1],
+                alphaX[rowMax, col-1],
+                alphaX[rowMax, col],
+                concentrations[rowMax, col+1],
+                concentrations[rowMax, col-1],
+                concentrations[rowMax, col],
+            )
     end
 
     # Handle corners
     # Top left corner
     topBoundaryValue = getBoundaryElementValue(bc, TOP, 1)
-    concentrations_t1[1, 1] += sx * calcHorizontalChange(alphaX[1, 2], alphaX[1, 1], alphaX[1, 1],
-        concentrations[1, 2], concentrations[1, 1], concentrations[1, 1]) +
-                               sy * calcVerticalChangeTopBoundary(topBoundaryType, alphaY[2, 1], alphaY[1, 1],
-        concentrations[2, 1], concentrations[1, 1], topBoundaryValue)
+    concentrations_t1[1, 1] +=
+        sx * calcHorizontalChange(
+            alphaX[1, 2],
+            alphaX[1, 1],
+            alphaX[1, 1],
+            concentrations[1, 2],
+            concentrations[1, 1],
+            concentrations[1, 1],
+        ) +
+        sy * calcVerticalChangeTopBoundary(
+            topBoundaryType,
+            alphaY[2, 1],
+            alphaY[1, 1],
+            concentrations[2, 1],
+            concentrations[1, 1],
+            topBoundaryValue,
+        )
 
     # Top right corner
     topBoundaryValue = getBoundaryElementValue(bc, TOP, colMax)
-    concentrations_t1[1, colMax] += sx * calcHorizontalChange(alphaX[1, colMax-1], alphaX[1, colMax], alphaX[1, colMax],
-        concentrations[1, colMax-1], concentrations[1, colMax], concentrations[1, colMax]) +
-                                    sy * calcVerticalChangeTopBoundary(topBoundaryType, alphaY[2, colMax], alphaY[1, colMax],
-        concentrations[2, colMax], concentrations[1, colMax], topBoundaryValue)
+    concentrations_t1[1, colMax] +=
+        sx * calcHorizontalChange(
+            alphaX[1, colMax-1],
+            alphaX[1, colMax],
+            alphaX[1, colMax],
+            concentrations[1, colMax-1],
+            concentrations[1, colMax],
+            concentrations[1, colMax],
+        ) +
+        sy * calcVerticalChangeTopBoundary(
+            topBoundaryType,
+            alphaY[2, colMax],
+            alphaY[1, colMax],
+            concentrations[2, colMax],
+            concentrations[1, colMax],
+            topBoundaryValue,
+        )
 
     # Bottom left corner
     bottomBoundaryValue = getBoundaryElementValue(bc, BOTTOM, 1)
-    concentrations_t1[rowMax, 1] += sx * calcHorizontalChange(alphaX[rowMax, 2], alphaX[rowMax, 1], alphaX[rowMax, 1],
-        concentrations[rowMax, 2], concentrations[rowMax, 1], concentrations[rowMax, 1]) +
-                                    sy * calcVerticalChangeBottomBoundary(bottomBoundaryType, alphaY[rowMax, 1], alphaY[rowMax-1, 1],
-        concentrations[rowMax, 1], concentrations[rowMax-1, 1], bottomBoundaryValue)
+    concentrations_t1[rowMax, 1] +=
+        sx * calcHorizontalChange(
+            alphaX[rowMax, 2],
+            alphaX[rowMax, 1],
+            alphaX[rowMax, 1],
+            concentrations[rowMax, 2],
+            concentrations[rowMax, 1],
+            concentrations[rowMax, 1],
+        ) +
+        sy * calcVerticalChangeBottomBoundary(
+            bottomBoundaryType,
+            alphaY[rowMax, 1],
+            alphaY[rowMax-1, 1],
+            concentrations[rowMax, 1],
+            concentrations[rowMax-1, 1],
+            bottomBoundaryValue,
+        )
 
     # Bottom right corner
     bottomBoundaryValue = getBoundaryElementValue(bc, BOTTOM, colMax)
-    concentrations_t1[rowMax, colMax] += sx * calcHorizontalChange(alphaX[rowMax, colMax-1], alphaX[rowMax, colMax], alphaX[rowMax, colMax],
-        concentrations[rowMax, colMax-1], concentrations[rowMax, colMax], concentrations[rowMax, colMax]) +
-                                         sy * calcVerticalChangeBottomBoundary(bottomBoundaryType, alphaY[rowMax, colMax], alphaY[rowMax-1, colMax],
-        concentrations[rowMax, colMax], concentrations[rowMax-1, colMax], bottomBoundaryValue)
+    concentrations_t1[rowMax, colMax] +=
+        sx * calcHorizontalChange(
+            alphaX[rowMax, colMax-1],
+            alphaX[rowMax, colMax],
+            alphaX[rowMax, colMax],
+            concentrations[rowMax, colMax-1],
+            concentrations[rowMax, colMax],
+            concentrations[rowMax, colMax],
+        ) +
+        sy * calcVerticalChangeBottomBoundary(
+            bottomBoundaryType,
+            alphaY[rowMax, colMax],
+            alphaY[rowMax-1, colMax],
+            concentrations[rowMax, colMax],
+            concentrations[rowMax-1, colMax],
+            bottomBoundaryValue,
+        )
 
     setConcentrations!(grid, concentrations_t1)
 end
 
-function runFTCS(grid::Grid{T}, bc::Boundary{T}, timestep::T, iterations::Int, stepCallback::Function) where {T}
+function runFTCS(
+    grid::Grid{T},
+    bc::Boundary{T},
+    timestep::T,
+    iterations::Int,
+    stepCallback::Function,
+) where {T}
     if getDim(grid) == 1
         for _ = 1:iterations
             FTCS_1D(grid, bc, timestep)
diff --git a/julia/TUG/src/DynamicSimulation.jl b/julia/TUG/src/DynamicSimulation.jl
index 756566d..f9cb461 100644
--- a/julia/TUG/src/DynamicSimulation.jl
+++ b/julia/TUG/src/DynamicSimulation.jl
@@ -26,7 +26,12 @@ struct DynamicSimulation{T} <: AbstractSimulation{T}
     iterations::Int
     timestep::T
 
-    function DynamicSimulation(grid::Grid{T}, bc::Boundary{T}, approach::APPROACH, timestep::T)::DynamicSimulation{T} where {T}
+    function DynamicSimulation(
+        grid::Grid{T},
+        bc::Boundary{T},
+        approach::APPROACH,
+        timestep::T,
+    )::DynamicSimulation{T} where {T}
         timestep, iterations = adjustTimestep(grid, approach, timestep, 1, false)
         new{T}(grid, Vector{Grid{T}}(), bc, approach, iterations, timestep)
     end
@@ -76,7 +81,10 @@ Prints the concentrations of a specific grid state to the console.
 # Returns
 Nothing, but outputs the grid's concentration matrix to the console.
 """
-function printConcentrations(simulation::DynamicSimulation{T}, gridIndex::Int)::Nothing where {T}
+function printConcentrations(
+    simulation::DynamicSimulation{T},
+    gridIndex::Int,
+)::Nothing where {T}
     writeConcentrationsToCLI(simulation, simulation.grids[gridIndex])
 end
 
@@ -92,17 +100,35 @@ Writes the concentrations of a specific grid state to a CSV file.
 # Returns
 Nothing, but writes the grid's concentration matrix to a CSV file.
 """
-function printConcentrationsCSV(simulation::DynamicSimulation{T}, gridIndex::Int)::Nothing where {T}
+function printConcentrationsCSV(
+    simulation::DynamicSimulation{T},
+    gridIndex::Int,
+)::Nothing where {T}
     file = createCSVfile(simulation, simulation.grids[gridIndex])
     writeConcentrationsToCSV(file, simulation, simulation.grids[gridIndex])
     close(file)
 end
 
-function runSimulationForGrid(simulation::DynamicSimulation{T}, grid::Grid{T})::Nothing where {T}
+function runSimulationForGrid(
+    simulation::DynamicSimulation{T},
+    grid::Grid{T},
+)::Nothing where {T}
     if simulation.approach == BTCS
-        runBTCS(grid, simulation.bc, simulation.timestep, simulation.iterations, () -> nothing)
+        runBTCS(
+            grid,
+            simulation.bc,
+            simulation.timestep,
+            simulation.iterations,
+            () -> nothing,
+        )
     elseif simulation.approach == FTCS
-        runFTCS(grid, simulation.bc, simulation.timestep, simulation.iterations, () -> nothing)
+        runFTCS(
+            grid,
+            simulation.bc,
+            simulation.timestep,
+            simulation.iterations,
+            () -> nothing,
+        )
     else
         error("Undefined approach!")
     end
@@ -120,7 +146,10 @@ Retrieves the concentrations of a specific grid state.
 # Returns
 The concentration matrix of the specified grid.
 """
-function getConcentrations(simulation::DynamicSimulation{T}, gridIndex::Int)::Matrix{T} where {T}
+function getConcentrations(
+    simulation::DynamicSimulation{T},
+    gridIndex::Int,
+)::Matrix{T} where {T}
     getConcentrations(simulation.grids[gridIndex])
 end
 
@@ -137,7 +166,11 @@ Sets the concentrations of a specific grid state.
 # Returns
 Nothing, but updates the concentration matrix of the specified grid.
 """
-function setConcentrations!(simulation::DynamicSimulation{T}, gridIndex::Int, concentrations::Matrix{T})::Nothing where {T}
+function setConcentrations!(
+    simulation::DynamicSimulation{T},
+    gridIndex::Int,
+    concentrations::Matrix{T},
+)::Nothing where {T}
     setConcentrations!(simulation.grids[gridIndex], concentrations)
 end
 
@@ -154,7 +187,11 @@ Sets the alphaX matrix of a specific grid state.
 # Returns
 Nothing, but updates the alphaX matrix of the specified grid.
 """
-function setAlphaX!(simulation::DynamicSimulation{T}, gridIndex::Int, alphaX::Matrix{T})::Nothing where {T}
+function setAlphaX!(
+    simulation::DynamicSimulation{T},
+    gridIndex::Int,
+    alphaX::Matrix{T},
+)::Nothing where {T}
     setAlphaX!(simulation.grids[gridIndex], alphaX)
 end
 
@@ -171,6 +208,10 @@ Sets the alphaY matrix of a specific grid state.
 # Returns
 Nothing, but updates the alphaY matrix of the specified grid.
 """
-function setAlphaY!(simulation::DynamicSimulation{T}, gridIndex::Int, alphaY::Matrix{T})::Nothing where {T}
+function setAlphaY!(
+    simulation::DynamicSimulation{T},
+    gridIndex::Int,
+    alphaY::Matrix{T},
+)::Nothing where {T}
     setAlphaY!(simulation.grids[gridIndex], alphaY)
 end
diff --git a/julia/TUG/src/Grid.jl b/julia/TUG/src/Grid.jl
index 321051f..f7074f7 100644
--- a/julia/TUG/src/Grid.jl
+++ b/julia/TUG/src/Grid.jl
@@ -46,30 +46,99 @@ struct Grid{T}
             throw(ArgumentError("Given grid length too small. Must be greater than 3."))
         end
         if size(alphaX, 1) != 1 || size(alphaX, 2) != length
-            throw(ArgumentError("Given matrix of alpha coefficients mismatch with Grid dimensions!"))
+            throw(
+                ArgumentError(
+                    "Given matrix of alpha coefficients mismatch with Grid dimensions!",
+                ),
+            )
         end
 
         alphaX_t = alphaX'
 
-        new{T}(length, 1, 1, T(length), 0, T(1), 0, Ref(fill(T(0), 1, length)), alphaX, nothing, alphaX_t, nothing)
+        new{T}(
+            length,
+            1,
+            1,
+            T(length),
+            0,
+            T(1),
+            0,
+            Ref(fill(T(0), 1, length)),
+            alphaX,
+            nothing,
+            alphaX_t,
+            nothing,
+        )
     end
 
-    function Grid{T}(rows::Int, cols::Int, alphaX::Matrix{T}, alphaY::Matrix{T})::Grid{T} where {T}
+    function Grid{T}(
+        rows::Int,
+        cols::Int,
+        alphaX::Matrix{T},
+        alphaY::Matrix{T},
+    )::Grid{T} where {T}
         if rows <= 3 || cols <= 3
-            throw(ArgumentError("Given grid dimensions too small. Must each be greater than 3."))
+            throw(
+                ArgumentError(
+                    "Given grid dimensions too small. Must each be greater than 3.",
+                ),
+            )
         end
         if size(alphaX) != (rows, cols) || size(alphaY) != (rows, cols)
-            throw(ArgumentError("Given matrices of alpha coefficients mismatch with Grid dimensions!"))
+            throw(
+                ArgumentError(
+                    "Given matrices of alpha coefficients mismatch with Grid dimensions!",
+                ),
+            )
         end
 
         alphaX_t = alphaX'
         alphaY_t = alphaY'
 
-        new{T}(cols, rows, 2, T(cols), T(rows), T(1), T(1), Ref(fill(T(0), rows, cols)), alphaX, alphaY, alphaX_t, alphaY_t)
+        new{T}(
+            cols,
+            rows,
+            2,
+            T(cols),
+            T(rows),
+            T(1),
+            T(1),
+            Ref(fill(T(0), rows, cols)),
+            alphaX,
+            alphaY,
+            alphaX_t,
+            alphaY_t,
+        )
     end
 
-    function Grid{T}(rows::Int, cols::Int, dim::Int, domainCol::T, domainRow::T, deltaCol::T, deltaRow::T, concentrations::Ref{Matrix{T}}, alphaX::Ref{Matrix{T}}, alphaY::Union{Ref{Matrix{T}},Nothing}, alphaX_t::Union{Ref{Matrix{T}},Nothing}, alphaY_t::Union{Ref{Matrix{T}},Nothing})::Grid{T} where {T}
-        new{T}(cols, rows, dim, domainCol, domainRow, deltaCol, deltaRow, concentrations, alphaX, alphaY, alphaX_t, alphaY_t)
+    function Grid{T}(
+        rows::Int,
+        cols::Int,
+        dim::Int,
+        domainCol::T,
+        domainRow::T,
+        deltaCol::T,
+        deltaRow::T,
+        concentrations::Ref{Matrix{T}},
+        alphaX::Ref{Matrix{T}},
+        alphaY::Union{Ref{Matrix{T}},Nothing},
+        alphaX_t::Union{Ref{Matrix{T}},Nothing},
+        alphaY_t::Union{Ref{Matrix{T}},Nothing},
+    )::Grid{T} where {T}
+        new{T}(
+            cols,
+            rows,
+            dim,
+            domainCol,
+            domainRow,
+            deltaCol,
+            deltaRow,
+            concentrations,
+            alphaX,
+            alphaY,
+            alphaX_t,
+            alphaY_t,
+        )
     end
 end
 
@@ -87,9 +156,35 @@ A new `Grid{T}` object that is a deep copy of the input grid.
 """
 function clone(grid::Grid{T})::Grid{T} where {T}
     if grid.dim == 1
-        return Grid{T}(1, grid.cols, grid.dim, grid.domainCol, grid.domainRow, grid.deltaCol, grid.deltaRow, Ref(copy(grid.concentrations[])), Ref(copy(grid.alphaX[])), nothing, Ref(copy(grid.alphaX_t[])), nothing)
+        return Grid{T}(
+            1,
+            grid.cols,
+            grid.dim,
+            grid.domainCol,
+            grid.domainRow,
+            grid.deltaCol,
+            grid.deltaRow,
+            Ref(copy(grid.concentrations[])),
+            Ref(copy(grid.alphaX[])),
+            nothing,
+            Ref(copy(grid.alphaX_t[])),
+            nothing,
+        )
     end
-    Grid{T}(grid.rows, grid.cols, grid.dim, grid.domainCol, grid.domainRow, grid.deltaCol, grid.deltaRow, Ref(copy(grid.concentrations[])), Ref(copy(grid.alphaX[])), Ref(copy(grid.alphaY[])), Ref(copy(grid.alphaX_t[])), Ref(copy(grid.alphaY_t[])))
+    Grid{T}(
+        grid.rows,
+        grid.cols,
+        grid.dim,
+        grid.domainCol,
+        grid.domainRow,
+        grid.deltaCol,
+        grid.deltaRow,
+        Ref(copy(grid.concentrations[])),
+        Ref(copy(grid.alphaX[])),
+        Ref(copy(grid.alphaY[])),
+        Ref(copy(grid.alphaX_t[])),
+        Ref(copy(grid.alphaY_t[])),
+    )
 end
 
 """
@@ -267,7 +362,11 @@ Nothing, but modifies the alphaX matrix of the grid.
 """
 function setAlphaX!(grid::Grid{T}, new_alphaX::Matrix{T})::Nothing where {T}
     if size(new_alphaX) != size(grid.alphaX[])
-        throw(ArgumentError("Given matrix of alpha coefficients mismatch with Grid dimensions!"))
+        throw(
+            ArgumentError(
+                "Given matrix of alpha coefficients mismatch with Grid dimensions!",
+            ),
+        )
     end
 
     grid.alphaX[] = new_alphaX
@@ -293,7 +392,11 @@ function setAlphaY!(grid::Grid{T}, new_alphaY::Matrix{T})::Nothing where {T}
         error("Grid is 1D, so there is no alphaY matrix!")
     end
     if size(new_alphaY) != size(grid.alphaY[])
-        throw(ArgumentError("Given matrix of alpha coefficients mismatch with Grid dimensions!"))
+        throw(
+            ArgumentError(
+                "Given matrix of alpha coefficients mismatch with Grid dimensions!",
+            ),
+        )
     end
 
     grid.alphaY[] = new_alphaY
@@ -316,7 +419,9 @@ Nothing, but modifies the concentration matrix of the grid.
 """
 function setConcentrations!(grid::Grid{T}, new_concentrations::Matrix{T})::Nothing where {T}
     if size(new_concentrations) != size(grid.concentrations[])
-        throw(ArgumentError("Given matrix of concentrations mismatch with Grid dimensions!"))
+        throw(
+            ArgumentError("Given matrix of concentrations mismatch with Grid dimensions!"),
+        )
     end
 
     grid.concentrations[] = new_concentrations
diff --git a/julia/TUG/src/Simulation.jl b/julia/TUG/src/Simulation.jl
index 6e27542..ee1c660 100644
--- a/julia/TUG/src/Simulation.jl
+++ b/julia/TUG/src/Simulation.jl
@@ -34,8 +34,15 @@ struct Simulation{T} <: AbstractSimulation{T}
     consoleOutput::CONSOLE_OUTPUT
     csvOutput::CSV_OUTPUT
 
-    function Simulation(grid::Grid{T}, bc::Boundary{T}, approach::APPROACH=BTCS, iterations::Int=1, timestep::T=0.1,
-        consoleOutput::CONSOLE_OUTPUT=CONSOLE_OUTPUT_OFF, csvOutput::CSV_OUTPUT=CSV_OUTPUT_OFF)::Simulation{T} where {T}
+    function Simulation(
+        grid::Grid{T},
+        bc::Boundary{T},
+        approach::APPROACH = BTCS,
+        iterations::Int = 1,
+        timestep::T = 0.1,
+        consoleOutput::CONSOLE_OUTPUT = CONSOLE_OUTPUT_OFF,
+        csvOutput::CSV_OUTPUT = CSV_OUTPUT_OFF,
+    )::Simulation{T} where {T}
         new{T}(grid, bc, approach, iterations, timestep, consoleOutput, csvOutput)
     end
 end
@@ -69,10 +76,28 @@ function run(simulation::Simulation{T})::Nothing where {T}
         end
 
         if simulation.approach == BTCS
-            runBTCS(simulation.grid, simulation.bc, simulation.timestep, simulation.iterations, simulationStepCallback)
+            runBTCS(
+                simulation.grid,
+                simulation.bc,
+                simulation.timestep,
+                simulation.iterations,
+                simulationStepCallback,
+            )
         elseif simulation.approach == FTCS
-            timestep, iterations = adjustTimestep(simulation.grid, simulation.approach, simulation.timestep, simulation.iterations, simulation.consoleOutput >= CONSOLE_OUTPUT_VERBOSE)
-            runFTCS(simulation.grid, simulation.bc, timestep, iterations, simulationStepCallback)
+            timestep, iterations = adjustTimestep(
+                simulation.grid,
+                simulation.approach,
+                simulation.timestep,
+                simulation.iterations,
+                simulation.consoleOutput >= CONSOLE_OUTPUT_VERBOSE,
+            )
+            runFTCS(
+                simulation.grid,
+                simulation.bc,
+                timestep,
+                iterations,
+                simulationStepCallback,
+            )
         else
             error("Undefined approach!")
         end
@@ -104,7 +129,15 @@ Sets the number of iterations for the given simulation.
 A new `Simulation` object with updated iterations.
 """
 function setIterations(simulation::Simulation{T}, iterations::Int)::Simulation{T} where {T}
-    return Simulation(simulation.grid, simulation.bc, simulation.approach, iterations, simulation.timestep, simulation.consoleOutput, simulation.csvOutput)
+    return Simulation(
+        simulation.grid,
+        simulation.bc,
+        simulation.approach,
+        iterations,
+        simulation.timestep,
+        simulation.consoleOutput,
+        simulation.csvOutput,
+    )
 end
 
 """
@@ -119,8 +152,19 @@ Sets the console output level for the simulation.
 # Returns
 A new `Simulation` object with updated console output setting.
 """
-function setOutputConsole(simulation::Simulation{T}, consoleOutput::CONSOLE_OUTPUT)::Simulation{T} where {T}
-    return Simulation(simulation.grid, simulation.bc, simulation.approach, simulation.iterations, simulation.timestep, consoleOutput, simulation.csvOutput)
+function setOutputConsole(
+    simulation::Simulation{T},
+    consoleOutput::CONSOLE_OUTPUT,
+)::Simulation{T} where {T}
+    return Simulation(
+        simulation.grid,
+        simulation.bc,
+        simulation.approach,
+        simulation.iterations,
+        simulation.timestep,
+        consoleOutput,
+        simulation.csvOutput,
+    )
 end
 
 """
@@ -135,8 +179,19 @@ Sets the CSV output level for the simulation.
 # Returns
 A new `Simulation` object with updated CSV output setting.
 """
-function setOutputCSV(simulation::Simulation{T}, csvOutput::CSV_OUTPUT)::Simulation{T} where {T}
-    return Simulation(simulation.grid, simulation.bc, simulation.approach, simulation.iterations, simulation.timestep, simulation.consoleOutput, csvOutput)
+function setOutputCSV(
+    simulation::Simulation{T},
+    csvOutput::CSV_OUTPUT,
+)::Simulation{T} where {T}
+    return Simulation(
+        simulation.grid,
+        simulation.bc,
+        simulation.approach,
+        simulation.iterations,
+        simulation.timestep,
+        simulation.consoleOutput,
+        csvOutput,
+    )
 end
 
 """
@@ -152,5 +207,13 @@ Sets the timestep for the simulation.
 A new `Simulation` object with updated timestep.
 """
 function setTimestep(simulation::Simulation{T}, timestep::T)::Simulation{T} where {T}
-    return Simulation(simulation.grid, simulation.bc, simulation.approach, simulation.iterations, timestep, simulation.consoleOutput, simulation.csvOutput)
+    return Simulation(
+        simulation.grid,
+        simulation.bc,
+        simulation.approach,
+        simulation.iterations,
+        timestep,
+        simulation.consoleOutput,
+        simulation.csvOutput,
+    )
 end
diff --git a/julia/TUG/src/TUG.jl b/julia/TUG/src/TUG.jl
index fa2afd4..069e646 100644
--- a/julia/TUG/src/TUG.jl
+++ b/julia/TUG/src/TUG.jl
@@ -3,16 +3,47 @@ module TUG
 include("Grid.jl")
 
 export Grid
-export clone, getAlphaX, getAlphaY, getAlphaX_t, getAlphaY_t, getConcentrations, setConcentrations!, setAlphaX!, setAlphaY!, getDomainCol, getDomainRow, getDeltaCol, getDeltaRow, getDim
+export clone,
+    getAlphaX,
+    getAlphaY,
+    getAlphaX_t,
+    getAlphaY_t,
+    getConcentrations,
+    setConcentrations!,
+    setAlphaX!,
+    setAlphaY!,
+    getDomainCol,
+    getDomainRow,
+    getDeltaCol,
+    getDeltaRow,
+    getDim
 
 include("Boundary.jl")
 
-export Boundary, BoundaryElement, TYPE, SIDE, BC_TYPE_CLOSED, BC_TYPE_CONSTANT, LEFT, RIGHT, TOP, BOTTOM
-export getType, getValue, setBoundarySideClosed!, setBoundarySideConstant!, getBoundarySide, getBoundaryElementType, getBoundaryElementValue
+export Boundary,
+    BoundaryElement, TYPE, SIDE, BC_TYPE_CLOSED, BC_TYPE_CONSTANT, LEFT, RIGHT, TOP, BOTTOM
+export getType,
+    getValue,
+    setBoundarySideClosed!,
+    setBoundarySideConstant!,
+    getBoundarySide,
+    getBoundaryElementType,
+    getBoundaryElementValue
 
 include("AbstractSimulation.jl")
 
-export APPROACH, CONSOLE_OUTPUT, CSV_OUTPUT, BTCS, FTCS, CONSOLE_OUTPUT_OFF, CONSOLE_OUTPUT_ON, CONSOLE_OUTPUT_VERBOSE, CSV_OUTPUT_OFF, CSV_OUTPUT_ON, CSV_OUTPUT_VERBOSE, CSV_OUTPUT_XTREME
+export APPROACH,
+    CONSOLE_OUTPUT,
+    CSV_OUTPUT,
+    BTCS,
+    FTCS,
+    CONSOLE_OUTPUT_OFF,
+    CONSOLE_OUTPUT_ON,
+    CONSOLE_OUTPUT_VERBOSE,
+    CSV_OUTPUT_OFF,
+    CSV_OUTPUT_ON,
+    CSV_OUTPUT_VERBOSE,
+    CSV_OUTPUT_XTREME
 
 include("Simulation.jl")
 
@@ -22,7 +53,14 @@ export run, setTimestep, setIterations, setOutputConsole, setOutputCSV
 include("DynamicSimulation.jl")
 
 export DynamicSimulation
-export createGrid, next, printConcentrations, printConcentrationsCSV, getConcentrations, setConcentrations!, setAlphaX!, setAlphaY!
+export createGrid,
+    next,
+    printConcentrations,
+    printConcentrationsCSV,
+    getConcentrations,
+    setConcentrations!,
+    setAlphaX!,
+    setAlphaY!
 
 include("Core/Utils.jl")
 include("Core/BTCS.jl")
diff --git a/julia/TUG/test/TestDynamicSimulation.jl b/julia/TUG/test/TestDynamicSimulation.jl
index a07f4d7..6b369ff 100644
--- a/julia/TUG/test/TestDynamicSimulation.jl
+++ b/julia/TUG/test/TestDynamicSimulation.jl
@@ -33,11 +33,16 @@
         boundary = TUG.Boundary(grid)
         simulation = TUG.DynamicSimulation(grid, boundary, BTCS, 0.01)
         TUG.createGrid(simulation)
-        for _ in 1:20
+        for _ = 1:20
             TUG.next(simulation)
         end
-        expected_concentrations = [1.281106278320615 3.5643693033301567 14.309048836698485 3.5643693033301598 1.281106278320616]
-        @test isapprox(TUG.getConcentrations(simulation, 1), expected_concentrations, atol=1e-6)
+        expected_concentrations =
+            [1.281106278320615 3.5643693033301567 14.309048836698485 3.5643693033301598 1.281106278320616]
+        @test isapprox(
+            TUG.getConcentrations(simulation, 1),
+            expected_concentrations,
+            atol = 1e-6,
+        )
 
         grid = TUG.Grid{Float64}(5, ones(1, 5))
         TUG.setConcentrations!(grid, [1.0 1.0 20.0 1.0 1.0])
@@ -45,34 +50,77 @@
         TUG.setBoundarySideConstant!(boundary, LEFT, 5.0)
         simulation = TUG.DynamicSimulation(grid, boundary, BTCS, 0.01)
         TUG.createGrid(simulation)
-        for _ in 1:20
+        for _ = 1:20
             TUG.next(simulation)
         end
-        expected_concentrations = [2.4416160635284823 3.6810808789967466 14.317333805802393 3.5648326408458035 1.2811288426376255]
-        @test isapprox(TUG.getConcentrations(simulation, 1), expected_concentrations, atol=1e-6)
+        expected_concentrations =
+            [2.4416160635284823 3.6810808789967466 14.317333805802393 3.5648326408458035 1.2811288426376255]
+        @test isapprox(
+            TUG.getConcentrations(simulation, 1),
+            expected_concentrations,
+            atol = 1e-6,
+        )
     end
     @testset "2D-Run" begin
         grid = TUG.Grid{Float64}(5, 5, ones(5, 5), ones(5, 5))
-        TUG.setConcentrations!(grid, [1.0 1.0 20.0 1.0 1.0; 1.0 1.0 20.0 1.0 1.0; 1.0 1.0 20.0 1.0 1.0; 1.0 1.0 20.0 1.0 1.0; 1.0 1.0 20.0 1.0 1.0])
+        TUG.setConcentrations!(
+            grid,
+            [
+                1.0 1.0 20.0 1.0 1.0
+                1.0 1.0 20.0 1.0 1.0
+                1.0 1.0 20.0 1.0 1.0
+                1.0 1.0 20.0 1.0 1.0
+                1.0 1.0 20.0 1.0 1.0
+            ],
+        )
         boundary = TUG.Boundary(grid)
         simulation = TUG.DynamicSimulation(grid, boundary, BTCS, 0.01)
         TUG.createGrid(simulation)
-        for _ in 1:20
+        for _ = 1:20
             TUG.next(simulation)
         end
-        expected_concentrations = [1.141904802011076 3.591390417498421 14.249599956958917 3.5913904174984217 1.1419048020110782; 1.1419048020110738 3.5913904174984173 14.2495999569589 3.5913904174984177 1.1419048020110767; 1.1419048020110725 3.591390417498413 14.249599956958875 3.5913904174984137 1.1419048020110751; 1.1419048020110738 3.5913904174984164 14.249599956958901 3.5913904174984173 1.141904802011077; 1.1419048020110774 3.5913904174984297 14.24959995695894 3.5913904174984297 1.1419048020110796]
-        @test isapprox(TUG.getConcentrations(simulation, 1), expected_concentrations, atol=1e-6)
+        expected_concentrations = [
+            1.141904802011076 3.591390417498421 14.249599956958917 3.5913904174984217 1.1419048020110782
+            1.1419048020110738 3.5913904174984173 14.2495999569589 3.5913904174984177 1.1419048020110767
+            1.1419048020110725 3.591390417498413 14.249599956958875 3.5913904174984137 1.1419048020110751
+            1.1419048020110738 3.5913904174984164 14.249599956958901 3.5913904174984173 1.141904802011077
+            1.1419048020110774 3.5913904174984297 14.24959995695894 3.5913904174984297 1.1419048020110796
+        ]
+        @test isapprox(
+            TUG.getConcentrations(simulation, 1),
+            expected_concentrations,
+            atol = 1e-6,
+        )
 
         grid = TUG.Grid{Float64}(5, 5, ones(5, 5), ones(5, 5))
-        TUG.setConcentrations!(grid, [1.0 1.0 20.0 1.0 1.0; 1.0 1.0 20.0 1.0 1.0; 1.0 1.0 20.0 1.0 1.0; 1.0 1.0 20.0 1.0 1.0; 1.0 1.0 20.0 1.0 1.0])
+        TUG.setConcentrations!(
+            grid,
+            [
+                1.0 1.0 20.0 1.0 1.0
+                1.0 1.0 20.0 1.0 1.0
+                1.0 1.0 20.0 1.0 1.0
+                1.0 1.0 20.0 1.0 1.0
+                1.0 1.0 20.0 1.0 1.0
+            ],
+        )
         boundary = TUG.Boundary(grid)
         TUG.setBoundarySideConstant!(boundary, LEFT, 5.0)
         simulation = TUG.DynamicSimulation(grid, boundary, BTCS, 0.01)
         TUG.createGrid(simulation)
-        for _ in 1:20
+        for _ = 1:20
             TUG.next(simulation)
         end
-        expected_concentrations = [1.9866377371338924 3.67421468453773 14.255058363518529 3.5916629034159486 1.1419105589005596; 1.98663773713389 3.674214684537723 14.255058363518497 3.5916629034159406 1.1419105589005576; 1.9866377371338884 3.6742146845377186 14.255058363518481 3.591662903415937 1.1419105589005565; 1.9866377371338895 3.674214684537725 14.255058363518502 3.5916629034159424 1.1419105589005574; 1.9866377371338952 3.6742146845377377 14.255058363518547 3.591662903415955 1.141910558900562]
-        @test isapprox(TUG.getConcentrations(simulation, 1), expected_concentrations, atol=1e-6)
+        expected_concentrations = [
+            1.9866377371338924 3.67421468453773 14.255058363518529 3.5916629034159486 1.1419105589005596
+            1.98663773713389 3.674214684537723 14.255058363518497 3.5916629034159406 1.1419105589005576
+            1.9866377371338884 3.6742146845377186 14.255058363518481 3.591662903415937 1.1419105589005565
+            1.9866377371338895 3.674214684537725 14.255058363518502 3.5916629034159424 1.1419105589005574
+            1.9866377371338952 3.6742146845377377 14.255058363518547 3.591662903415955 1.141910558900562
+        ]
+        @test isapprox(
+            TUG.getConcentrations(simulation, 1),
+            expected_concentrations,
+            atol = 1e-6,
+        )
     end
-end
\ No newline at end of file
+end
diff --git a/julia/TUG/test/TestSimulation.jl b/julia/TUG/test/TestSimulation.jl
index edf26a9..3367454 100644
--- a/julia/TUG/test/TestSimulation.jl
+++ b/julia/TUG/test/TestSimulation.jl
@@ -13,7 +13,8 @@
 
         grid = TUG.Grid{Float64}(5, ones(1, 5))
         boundary = TUG.Boundary(grid)
-        simulation = TUG.Simulation(grid, boundary, FTCS, 2, 0.2, CONSOLE_OUTPUT_ON, CSV_OUTPUT_ON)
+        simulation =
+            TUG.Simulation(grid, boundary, FTCS, 2, 0.2, CONSOLE_OUTPUT_ON, CSV_OUTPUT_ON)
         @test simulation.grid == grid
         @test simulation.bc == boundary
         @test simulation.approach == FTCS
@@ -28,8 +29,9 @@
         boundary = TUG.Boundary(grid)
         simulation = TUG.Simulation(grid, boundary, BTCS, 20, 0.01)
         TUG.run(simulation)
-        expected_concentrations = [1.281106278320615 3.5643693033301567 14.309048836698485 3.5643693033301598 1.281106278320616]
-        @test isapprox(TUG.getConcentrations(grid), expected_concentrations, atol=1e-6)
+        expected_concentrations =
+            [1.281106278320615 3.5643693033301567 14.309048836698485 3.5643693033301598 1.281106278320616]
+        @test isapprox(TUG.getConcentrations(grid), expected_concentrations, atol = 1e-6)
 
         grid = TUG.Grid{Float64}(5, ones(1, 5))
         TUG.setConcentrations!(grid, [1.0 1.0 20.0 1.0 1.0])
@@ -37,25 +39,56 @@
         TUG.setBoundarySideConstant!(boundary, LEFT, 5.0)
         simulation = TUG.Simulation(grid, boundary, BTCS, 20, 0.01)
         TUG.run(simulation)
-        expected_concentrations = [2.4416160635284823 3.6810808789967466 14.317333805802393 3.5648326408458035 1.2811288426376255]
-        @test isapprox(TUG.getConcentrations(grid), expected_concentrations, atol=1e-6)
+        expected_concentrations =
+            [2.4416160635284823 3.6810808789967466 14.317333805802393 3.5648326408458035 1.2811288426376255]
+        @test isapprox(TUG.getConcentrations(grid), expected_concentrations, atol = 1e-6)
     end
     @testset "2D-Run" begin
         grid = TUG.Grid{Float64}(5, 5, ones(5, 5), ones(5, 5))
-        TUG.setConcentrations!(grid, [1.0 1.0 20.0 1.0 1.0; 1.0 1.0 20.0 1.0 1.0; 1.0 1.0 20.0 1.0 1.0; 1.0 1.0 20.0 1.0 1.0; 1.0 1.0 20.0 1.0 1.0])
+        TUG.setConcentrations!(
+            grid,
+            [
+                1.0 1.0 20.0 1.0 1.0
+                1.0 1.0 20.0 1.0 1.0
+                1.0 1.0 20.0 1.0 1.0
+                1.0 1.0 20.0 1.0 1.0
+                1.0 1.0 20.0 1.0 1.0
+            ],
+        )
         boundary = TUG.Boundary(grid)
         simulation = TUG.Simulation(grid, boundary, BTCS, 20, 0.01)
         TUG.run(simulation)
-        expected_concentrations = [1.141904802011076 3.591390417498421 14.249599956958917 3.5913904174984217 1.1419048020110782; 1.1419048020110738 3.5913904174984173 14.2495999569589 3.5913904174984177 1.1419048020110767; 1.1419048020110725 3.591390417498413 14.249599956958875 3.5913904174984137 1.1419048020110751; 1.1419048020110738 3.5913904174984164 14.249599956958901 3.5913904174984173 1.141904802011077; 1.1419048020110774 3.5913904174984297 14.24959995695894 3.5913904174984297 1.1419048020110796]
-        @test isapprox(TUG.getConcentrations(grid), expected_concentrations, atol=1e-6)
+        expected_concentrations = [
+            1.141904802011076 3.591390417498421 14.249599956958917 3.5913904174984217 1.1419048020110782
+            1.1419048020110738 3.5913904174984173 14.2495999569589 3.5913904174984177 1.1419048020110767
+            1.1419048020110725 3.591390417498413 14.249599956958875 3.5913904174984137 1.1419048020110751
+            1.1419048020110738 3.5913904174984164 14.249599956958901 3.5913904174984173 1.141904802011077
+            1.1419048020110774 3.5913904174984297 14.24959995695894 3.5913904174984297 1.1419048020110796
+        ]
+        @test isapprox(TUG.getConcentrations(grid), expected_concentrations, atol = 1e-6)
 
         grid = TUG.Grid{Float64}(5, 5, ones(5, 5), ones(5, 5))
-        TUG.setConcentrations!(grid, [1.0 1.0 20.0 1.0 1.0; 1.0 1.0 20.0 1.0 1.0; 1.0 1.0 20.0 1.0 1.0; 1.0 1.0 20.0 1.0 1.0; 1.0 1.0 20.0 1.0 1.0])
+        TUG.setConcentrations!(
+            grid,
+            [
+                1.0 1.0 20.0 1.0 1.0
+                1.0 1.0 20.0 1.0 1.0
+                1.0 1.0 20.0 1.0 1.0
+                1.0 1.0 20.0 1.0 1.0
+                1.0 1.0 20.0 1.0 1.0
+            ],
+        )
         boundary = TUG.Boundary(grid)
         TUG.setBoundarySideConstant!(boundary, LEFT, 5.0)
         simulation = TUG.Simulation(grid, boundary, BTCS, 20, 0.01)
         TUG.run(simulation)
-        expected_concentrations = [1.9866377371338924 3.67421468453773 14.255058363518529 3.5916629034159486 1.1419105589005596; 1.98663773713389 3.674214684537723 14.255058363518497 3.5916629034159406 1.1419105589005576; 1.9866377371338884 3.6742146845377186 14.255058363518481 3.591662903415937 1.1419105589005565; 1.9866377371338895 3.674214684537725 14.255058363518502 3.5916629034159424 1.1419105589005574; 1.9866377371338952 3.6742146845377377 14.255058363518547 3.591662903415955 1.141910558900562]
-        @test isapprox(TUG.getConcentrations(grid), expected_concentrations, atol=1e-6)
+        expected_concentrations = [
+            1.9866377371338924 3.67421468453773 14.255058363518529 3.5916629034159486 1.1419105589005596
+            1.98663773713389 3.674214684537723 14.255058363518497 3.5916629034159406 1.1419105589005576
+            1.9866377371338884 3.6742146845377186 14.255058363518481 3.591662903415937 1.1419105589005565
+            1.9866377371338895 3.674214684537725 14.255058363518502 3.5916629034159424 1.1419105589005574
+            1.9866377371338952 3.6742146845377377 14.255058363518547 3.591662903415955 1.141910558900562
+        ]
+        @test isapprox(TUG.getConcentrations(grid), expected_concentrations, atol = 1e-6)
     end
 end