max/TugJulia
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

style: applied linting recommendations

Browse Source
This commit is contained in:
nebmit 2023-12-04 19:21:23 +01:00
parent 1e57d8b5b5
commit 39eacff904
No known key found for this signature in database
10 changed files with 922 additions and 178 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

45
julia/TUG/src/AbstractSimulation.jl
View File

@ -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)

64
julia/TUG/src/Boundary.jl
View File

@ -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

243
julia/TUG/src/Core/BTCS.jl
View File

@ -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

290
julia/TUG/src/Core/FTCS.jl
View File

@ -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)

61
julia/TUG/src/DynamicSimulation.jl
View File

@ -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

131
julia/TUG/src/Grid.jl
View File

@ -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

85
julia/TUG/src/Simulation.jl
View File

@ -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

48
julia/TUG/src/TUG.jl
View File

@ -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")

78
julia/TUG/test/TestDynamicSimulation.jl
View File

@ -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
end

55
julia/TUG/test/TestSimulation.jl
View File

@ -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