docs: added docstrings for public functions

[skip ci]
This commit is contained in:
nebmit 2023-12-04 19:16:36 +01:00
parent 01b2247344
commit 1e57d8b5b5
No known key found for this signature in database
9 changed files with 610 additions and 111 deletions

View File

@ -1,7 +1,33 @@
using Printf using Printf
"""
@enum APPROACH BTCS FTCS
Enum representing the numerical approach for the simulation.
- `BTCS`: Backward Time, Central Space. A stable, implicit method.
- `FTCS`: Forward Time, Central Space. A conditionally stable, explicit method.
"""
@enum APPROACH BTCS FTCS @enum APPROACH BTCS FTCS
"""
@enum CONSOLE_OUTPUT CONSOLE_OUTPUT_OFF CONSOLE_OUTPUT_ON CONSOLE_OUTPUT_VERBOSE
Enum representing the level of console output for the simulation.
- `CONSOLE_OUTPUT_OFF`: No output to the console.
- `CONSOLE_OUTPUT_ON`: Standard output level.
- `CONSOLE_OUTPUT_VERBOSE`: Detailed output, including intermediate steps.
"""
@enum CONSOLE_OUTPUT CONSOLE_OUTPUT_OFF CONSOLE_OUTPUT_ON CONSOLE_OUTPUT_VERBOSE @enum CONSOLE_OUTPUT CONSOLE_OUTPUT_OFF CONSOLE_OUTPUT_ON CONSOLE_OUTPUT_VERBOSE
"""
@enum CSV_OUTPUT CSV_OUTPUT_OFF CSV_OUTPUT_ON CSV_OUTPUT_VERBOSE CSV_OUTPUT_XTREME
Enum representing the level of CSV file output for the simulation.
- `CSV_OUTPUT_OFF`: No output to CSV.
- `CSV_OUTPUT_ON`: Standard output level.
- `CSV_OUTPUT_VERBOSE`: Detailed output.
- `CSV_OUTPUT_XTREME`: Extremely detailed output, possibly including debugging information.
"""
@enum CSV_OUTPUT CSV_OUTPUT_OFF CSV_OUTPUT_ON CSV_OUTPUT_VERBOSE CSV_OUTPUT_XTREME @enum CSV_OUTPUT CSV_OUTPUT_OFF CSV_OUTPUT_ON CSV_OUTPUT_VERBOSE CSV_OUTPUT_XTREME
abstract type AbstractSimulation{T} end abstract type AbstractSimulation{T} end
@ -38,7 +64,7 @@ function adjustTimestep(grid::Grid{T}, approach::APPROACH, timestep::T, iteratio
return timestep, iterations return timestep, iterations
end end
function createCSVfile(simulation::AbstractSimulation{T}, grid::Union{Grid{T},Nothing}=nothing)::IOStream where {T} function createCSVfile(simulation::AbstractSimulation{T}, grid::Union{Grid{T},Nothing}=nothing) where {T}
grid = grid === nothing ? simulation.grid : grid grid = grid === nothing ? simulation.grid : grid
approachString = string(simulation.approach) approachString = string(simulation.approach)
rows = getRows(grid) rows = getRows(grid)
@ -68,8 +94,7 @@ function createCSVfile(simulation::AbstractSimulation{T}, grid::Union{Grid{T},No
end end
end end
file = open(filename, "a") open(filename, "a")
return file
end end
function writeBoundarySideValuesToCSV(file::IOStream, bc::Boundary{T}, side) where {T} function writeBoundarySideValuesToCSV(file::IOStream, bc::Boundary{T}, side) where {T}

View File

@ -3,18 +3,49 @@
# condition at the edges of the diffusion grid. # condition at the edges of the diffusion grid.
# Translated from C++'s Boundary.hpp. # Translated from C++'s Boundary.hpp.
"""
@enum TYPE CLOSED CONSTANT
Enum representing the type of boundary.
- `CLOSED`: Closed boundary condition.
- `CONSTANT`: Constant boundary condition.
"""
@enum TYPE CLOSED CONSTANT @enum TYPE CLOSED CONSTANT
"""
@enum SIDE LEFT RIGHT TOP BOTTOM
Enum representing the side of the boundary.
- `LEFT`: Left side.
- `RIGHT`: Right side.
- `TOP`: Top side.
- `BOTTOM`: Bottom side.
"""
@enum SIDE LEFT = 1 RIGHT = 2 TOP = 3 BOTTOM = 4 @enum SIDE LEFT = 1 RIGHT = 2 TOP = 3 BOTTOM = 4
"""
struct BoundaryElement{T}
Represents an element of a boundary condition in a diffusion grid.
It holds the type of boundary condition and its value if applicable.
# Fields
- `type::TYPE`: The type of boundary condition (CLOSED, CONSTANT).
- `value::T`: The value of the boundary condition, used if type is CONSTANT.
# Constructors
- `BoundaryElement{T}()` creates a CLOSED boundary element.
- `BoundaryElement{T}(value::T)` creates a CONSTANT boundary element with specified value.
"""
struct BoundaryElement{T} struct BoundaryElement{T}
type::TYPE type::TYPE
value::T value::T
# Generic constructor for closed case function BoundaryElement{T}()::BoundaryElement{T} where {T}
BoundaryElement{T}() where {T} = new{T}(CLOSED, convert(T, -1)) new{T}(CLOSED, convert(T, -1))
end
# Constructor for constant case function BoundaryElement{T}(value::T)::BoundaryElement{T} where {T}
function BoundaryElement{T}(value::T) where {T}
if value < 0 if value < 0
throw(ArgumentError("No negative concentration allowed.")) throw(ArgumentError("No negative concentration allowed."))
end end
@ -22,26 +53,72 @@ struct BoundaryElement{T}
end end
end end
"""
getType(be::BoundaryElement{T})::TYPE where {T}
Retrieves the type of the boundary element.
# Arguments
- `be::BoundaryElement{T}`: The boundary element.
# Returns
The type of the boundary element (`TYPE`).
"""
function getType(be::BoundaryElement{T})::TYPE where {T} function getType(be::BoundaryElement{T})::TYPE where {T}
be.type be.type
end end
"""
getValue(be::BoundaryElement{T})::T where {T}
Retrieves the value of the boundary element. Applicable if the boundary type is CONSTANT.
# Arguments
- `be::BoundaryElement{T}`: The boundary element.
# Returns
The value of the boundary element.
"""
function getValue(be::BoundaryElement{T})::T where {T} function getValue(be::BoundaryElement{T})::T where {T}
be.value be.value
end end
"""
getValue(be::Vector{BoundaryElement{T}})::Vector{T} where {T}
Retrieves the values of the boundary elements. Applicable if the boundary type is CONSTANT.
# Arguments
- `be::Vector{BoundaryElement{T}}`: The boundary elements.
# Returns
The values of the boundary elements.
"""
function getValue(be::Vector{BoundaryElement{T}})::Vector{T} where {T} function getValue(be::Vector{BoundaryElement{T}})::Vector{T} where {T}
[b.value for b in be] [b.value for b in be]
end end
"""
struct Boundary{T}
Represents the boundary conditions of a diffusion grid.
It includes the dimension of the grid and boundary elements for each side.
# Fields
- `dim::UInt8`: Dimension of the grid (1 or 2).
- `cols::UInt32`, `rows::UInt32`: Number of columns and rows in the grid.
- `boundaries::Vector{Vector{BoundaryElement{T}}}`: A vector of boundary elements for each side.
# Constructor
- `Boundary(grid::Grid{T})` constructs the boundary based on the given grid's dimensions.
"""
struct Boundary{T} struct Boundary{T}
dim::UInt8 dim::UInt8
cols::UInt32 cols::UInt32
rows::UInt32 rows::UInt32
boundaries::Vector{Vector{BoundaryElement{T}}} boundaries::Vector{Vector{BoundaryElement{T}}}
# Constructor function Boundary(grid::Grid{T})::Boundary{T} where {T}
function Boundary(grid::Grid{T}) where {T}
dim = grid.dim dim = grid.dim
cols = grid.cols cols = grid.cols
rows = grid.rows rows = grid.rows
@ -63,6 +140,19 @@ struct Boundary{T}
end end
end end
"""
getBoundaryElementType(boundary::Boundary{T}, side::SIDE, index::Int)::TYPE where {T}
Retrieves the type of the boundary element at the specified side and index.
# Arguments
- `boundary::Boundary{T}`: The boundary.
- `side::SIDE`: The side of the boundary.
- `index::Int`: The index of the boundary element.
# 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) 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."))
@ -74,6 +164,20 @@ function getBoundaryElementType(boundary::Boundary{T}, side::SIDE, index::Int)::
getType(boundary.boundaries[Int(side)][index]) getType(boundary.boundaries[Int(side)][index])
end end
"""
getBoundaryElementValue(boundary::Boundary{T}, side::SIDE, index::Int)::T where {T}
Retrieves the value of the boundary element at the specified side and index.
Applicable if the boundary type is CONSTANT.
# Arguments
- `boundary::Boundary{T}`: The boundary.
- `side::SIDE`: The side of the boundary.
- `index::Int`: The index of the boundary element.
# Returns
The value of the boundary element.
"""
function getBoundaryElementValue(boundary::Boundary{T}, side::SIDE, index::Int)::T where {T} function getBoundaryElementValue(boundary::Boundary{T}, side::SIDE, index::Int)::T where {T}
if boundary.dim == 1 && (side == BOTTOM || side == TOP) 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."))
@ -85,6 +189,18 @@ function getBoundaryElementValue(boundary::Boundary{T}, side::SIDE, index::Int):
getValue(boundary.boundaries[Int(side)][index]) getValue(boundary.boundaries[Int(side)][index])
end end
"""
getBoundarySide(boundary::Boundary{T}, side::SIDE)::Vector{BoundaryElement{T}} where {T}
Retrieves the boundary elements at the specified side.
# Arguments
- `boundary::Boundary{T}`: The boundary.
- `side::SIDE`: The side of the boundary.
# 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) 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."))
@ -93,7 +209,19 @@ function getBoundarySide(boundary::Boundary{T}, side::SIDE)::Vector{BoundaryElem
boundary.boundaries[Int(side)] boundary.boundaries[Int(side)]
end end
function setBoundarySideClosed!(boundary::Boundary{T}, side::SIDE) where {T} """
setBoundarySideClosed!(boundary::Boundary{T}, side::SIDE)::Vector{BoundaryElement{T}} where {T}
Sets the boundary elements at the specified side to CLOSED.
# Arguments
- `boundary::Boundary{T}`: The boundary.
- `side::SIDE`: The side of the boundary.
# Returns
The boundary elements at the specified side.
"""
function setBoundarySideClosed!(boundary::Boundary{T}, side::SIDE)::Vector{BoundaryElement{T}} where {T}
if boundary.dim == 1 && (side == BOTTOM || side == TOP) 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 end
@ -104,7 +232,20 @@ function setBoundarySideClosed!(boundary::Boundary{T}, side::SIDE) where {T}
boundary.boundaries[Int(side)] = [BoundaryElement{T}() for _ in 1:n] boundary.boundaries[Int(side)] = [BoundaryElement{T}() for _ in 1:n]
end end
function setBoundarySideConstant!(boundary::Boundary{T}, side::SIDE, value::T) where {T} """
setBoundarySideConstant!(boundary::Boundary{T}, side::SIDE, values::Vector{T})::Vector{BoundaryElement{T}} where {T}
Sets the boundary elements at the specified side to CONSTANT with the specified values.
# Arguments
- `boundary::Boundary{T}`: The boundary.
- `side::SIDE`: The side of the boundary.
- `value::T`: The value of the boundary elements.
# Returns
The boundary elements at the specified side.
"""
function setBoundarySideConstant!(boundary::Boundary{T}, side::SIDE, value::T)::Vector{BoundaryElement{T}} where {T}
if boundary.dim == 1 && (side == BOTTOM || side == TOP) 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 end

View File

@ -8,78 +8,51 @@ using Base.Threads
using LinearAlgebra using LinearAlgebra
using SparseArrays using SparseArrays
function calcBoundaryCoeffClosed(alpha_center::T, alpha_side::T, sx::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) alpha = calcAlphaIntercell(alpha_center, alpha_side)
centerCoeff = 1 + sx * alpha
sideCoeff = -sx * alpha sideCoeff = -sx * alpha
return (centerCoeff, sideCoeff)
end
function calcBoundaryCoeffConstant(alpha_center::T, alpha_side::T, sx::T) where {T} if bcType == CONSTANT
alpha = calcAlphaIntercell(alpha_center, alpha_side) centerCoeff = 1 + sx * (alpha + 2 * alpha_center)
centerCoeff = 1 + sx * (alpha + 2 * alpha_center) elseif bcType == CLOSED
sideCoeff = -sx * alpha centerCoeff = 1 + sx * alpha
return (centerCoeff, sideCoeff)
end
# creates coefficient matrix for next time step from alphas in x-direction
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}
# Precompute boundary condition type check for efficiency
bcTypeLeft = getType(bcLeft[rowIndex])
# Determine left side boundary coefficients based on boundary condition
centerCoeffTop, rightCoeffTop = if bcTypeLeft == CONSTANT
calcBoundaryCoeffConstant(alpha[rowIndex, 1], alpha[rowIndex, 2], sx)
elseif bcTypeLeft == CLOSED
calcBoundaryCoeffClosed(alpha[rowIndex, 1], alpha[rowIndex, 2], sx)
else else
error("Undefined Boundary Condition Type on Left!") error("Undefined Boundary Condition Type!")
end end
# Precompute boundary condition type check for efficiency return (centerCoeff, sideCoeff)
bcTypeRight = getType(bcRight[rowIndex]) 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}
# Determine left side boundary coefficients based on boundary condition
centerCoeffTop, rightCoeffTop = calcBoundaryCoeff(alpha[rowIndex, 1], alpha[rowIndex, 2], sx, getType(bcLeft[rowIndex]))
# Determine right side boundary coefficients based on boundary condition # Determine right side boundary coefficients based on boundary condition
centerCoeffBottom, leftCoeffBottom = if bcTypeRight == CONSTANT centerCoeffBottom, leftCoeffBottom = calcBoundaryCoeff(alpha[rowIndex, numCols], alpha[rowIndex, numCols-1], sx, getType(bcRight[rowIndex]))
calcBoundaryCoeffConstant(alpha[rowIndex, numCols], alpha[rowIndex, numCols-1], sx)
elseif bcTypeRight == CLOSED
calcBoundaryCoeffClosed(alpha[rowIndex, numCols], alpha[rowIndex, numCols-1], sx)
else
error("Undefined Boundary Condition Type on Right!")
end
dl = [-sx .* alpha_left; leftCoeffBottom] dl = [-sx .* alpha_left; leftCoeffBottom]
d = [centerCoeffTop; 1 .+ sx .* (alpha_right + alpha_left); centerCoeffBottom] d = [centerCoeffTop; 1 .+ sx .* (alpha_right + alpha_left); centerCoeffBottom]
du = [rightCoeffTop; -sx .* alpha_right] du = [rightCoeffTop; -sx .* alpha_right]
alpha_diagonal = Tridiagonal(dl, d, du) return Tridiagonal(dl, d, du)
return alpha_diagonal
end end
function calcExplicitConcentrationsBoundaryClosed(conc_center::T, alpha_center::T, alpha_neighbor::T, sy::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) alpha = calcAlphaIntercell(alpha_center, alpha_neighbor)
return (sy * alpha + (1 - sy * alpha)) * conc_center
(sy * alpha + (1 - sy * alpha)) * conc_center
end end
function calcExplicitConcentrationsBoundaryConstant(conc_center::T, conc_bc::T, alpha_center::T, alpha_neighbor::T, sy::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_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)
sy * alpha_center_neighbor * conc_center + return sy * alpha_center_neighbor * conc_center +
(1 - sy * (alpha_center_center + 2 * alpha_center)) * conc_center + (1 - sy * (alpha_center_center + 2 * alpha_center)) * conc_center +
sy * alpha_center * conc_bc sy * alpha_center * conc_bc
end end
function calcExplicitConcentrationsBoundaryClosed(conc_center::Vector{T}, alpha::Vector{T}, sy::T) where {T} function calcExplicitConcentrationsBoundaryClosed(conc_center::Vector{T}, alpha::Vector{T}, sy::T)::T where {T}
(sy .* alpha .+ (1 .- sy .* alpha)) .* conc_center return (sy .* alpha .+ (1 .- sy .* alpha)) .* conc_center
end
function calcExplicitConcentrationsBoundaryConstant(conc_center::Vector{T}, conc_bc::Vector{T}, alpha_center::Vector{T}, alpha_neighbor::Vector{T}, sy::T) where {T}
sy .* alpha_neighbor .* conc_center[rowIndex, :] +
(1 .- sy .* (alpha_center .+ 2 .* alphaY[rowIndex, :])) .* conc_center[rowIndex, :] +
sy .* alphaY[rowIndex, :] .* conc_bc
end end
# creates a solution vector for next time step from the current state of concentrations # creates a solution vector for next time step from the current state of concentrations
@ -87,17 +60,6 @@ function writeSolutionVector!(sv::Vector{T}, concentrations::Matrix{T}, alphaX::
numRows = size(concentrations, 1) numRows = size(concentrations, 1)
length = size(sv, 1) length = size(sv, 1)
# 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]
end
end
# First row # First row
if rowIndex == 1 if rowIndex == 1
@inbounds for i = 1:length @inbounds for i = 1:length
@ -111,6 +73,17 @@ function writeSolutionVector!(sv::Vector{T}, concentrations::Matrix{T}, alphaX::
end end
end end
# 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]
end
end
# Last row # Last row
if rowIndex == numRows if rowIndex == numRows
@inbounds for i = 1:length @inbounds for i = 1:length
@ -135,9 +108,6 @@ function writeSolutionVector!(sv::Vector{T}, concentrations::Matrix{T}, alphaX::
end end
end end
# BTCS solution for 1D grid
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)) sx = timestep / (getDeltaCol(grid) * getDeltaCol(grid))
@ -161,7 +131,6 @@ function BTCS_1D(grid::Grid{T}, bc::Boundary{T}, alpha_left::Matrix{T}, alpha_ri
concentrations[1, :] = A \ b concentrations[1, :] = A \ b
end end
# BTCS solution for 2D grid
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) rows = getRows(grid)
cols = getCols(grid) cols = getCols(grid)
@ -183,7 +152,6 @@ function BTCS_2D(grid::Grid{T}, bc::Boundary{T}, alphaX_left::Matrix{T}, alphaX_
bcTop = getBoundarySide(bc, TOP) bcTop = getBoundarySide(bc, TOP)
bcBottom = getBoundarySide(bc, BOTTOM) bcBottom = getBoundarySide(bc, BOTTOM)
localBs = [zeros(T, cols) for _ in 1:Threads.nthreads()] localBs = [zeros(T, cols) for _ in 1:Threads.nthreads()]
Threads.@threads for i = 1:rows Threads.@threads for i = 1:rows
localB = localBs[Threads.threadid()] localB = localBs[Threads.threadid()]
@ -193,11 +161,9 @@ function BTCS_2D(grid::Grid{T}, bc::Boundary{T}, alphaX_left::Matrix{T}, alphaX_
concentrations_intermediate[i, :] = A \ localB concentrations_intermediate[i, :] = A \ localB
end end
concentrations_intermediate = copy(transpose(concentrations_intermediate)) concentrations_intermediate = copy(transpose(concentrations_intermediate))
concentrations_t = fetch(concentrations_t_task) concentrations_t = fetch(concentrations_t_task)
# Swap alphas, boundary conditions and sx/sy for column-wise calculation # 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 _ in 1:Threads.nthreads()]
Threads.@threads for i = 1:cols Threads.@threads for i = 1:cols
@ -253,5 +219,4 @@ function runBTCS(grid::Grid{T}, bc::Boundary{T}, timestep::T, iterations::Int, s
else else
error("BTCS only implemented for 1D and 2D grids") error("BTCS only implemented for 1D and 2D grids")
end end
end end

View File

@ -89,8 +89,6 @@ function calcVerticalChangeBottomBoundary(boundaryType::TYPE, alphaY_current::T,
end end
end end
function FTCS_1D(grid::Grid{T}, bc::Boundary{T}, timestep::T) where {T} function FTCS_1D(grid::Grid{T}, bc::Boundary{T}, timestep::T) where {T}
colMax = getCols(grid) colMax = getCols(grid)
sx = timestep / (getDeltaCol(grid)^2) sx = timestep / (getDeltaCol(grid)^2)
@ -122,7 +120,6 @@ function FTCS_1D(grid::Grid{T}, bc::Boundary{T}, timestep::T) where {T}
setConcentrations!(grid, concentrations_t1) setConcentrations!(grid, concentrations_t1)
end end
function FTCS_2D(grid::Grid{T}, bc::Boundary{T}, timestep::T) where {T} function FTCS_2D(grid::Grid{T}, bc::Boundary{T}, timestep::T) where {T}
rowMax = getRows(grid) rowMax = getRows(grid)
colMax = getCols(grid) colMax = getCols(grid)

View File

@ -2,10 +2,6 @@ function calcAlphaIntercell(alpha1::T, alpha2::T) where {T}
2 / ((1 / alpha1) + (1 / alpha2)) 2 / ((1 / alpha1) + (1 / alpha2))
end end
function calcAlphaIntercell(alpha1::Vector{T}, alpha2::Vector{T}) where {T}
2 ./ ((1 ./ alpha1) .+ (1 ./ alpha2))
end
function calcAlphaIntercell(alpha1::Matrix{T}, alpha2::Matrix{T}) where {T} function calcAlphaIntercell(alpha1::Matrix{T}, alpha2::Matrix{T}) where {T}
2 ./ ((1 ./ alpha1) .+ (1 ./ alpha2)) 2 ./ ((1 ./ alpha1) .+ (1 ./ alpha2))
end end

View File

@ -1,5 +1,22 @@
using Distributed using Distributed
"""
struct DynamicSimulation{T} <: AbstractSimulation{T}
Represents a dynamic simulation environment with multiple grid states.
Allows the manipulation and running of simulations on different grid states in parallel.
# Fields
- `grid::Grid{T}`: The initial grid state.
- `grids::Vector{Grid{T}}`: A vector of grid states for the simulation.
- `bc::Boundary{T}`: The boundary conditions for the simulation.
- `approach::APPROACH`: The numerical approach (BTCS or FTCS).
- `iterations::Int`: Number of iterations to run each simulation.
- `timestep::T`: The timestep for each iteration.
# Constructor
- `DynamicSimulation(grid, bc, approach, timestep)` creates a new dynamic simulation with specified parameters.
"""
struct DynamicSimulation{T} <: AbstractSimulation{T} struct DynamicSimulation{T} <: AbstractSimulation{T}
grid::Grid{T} grid::Grid{T}
grids::Vector{Grid{T}} grids::Vector{Grid{T}}
@ -9,34 +26,79 @@ struct DynamicSimulation{T} <: AbstractSimulation{T}
iterations::Int iterations::Int
timestep::T timestep::T
# Constructor 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) where {T}
timestep, iterations = adjustTimestep(grid, approach, timestep, 1, false) timestep, iterations = adjustTimestep(grid, approach, timestep, 1, false)
new{T}(grid, Vector{Grid{T}}(), bc, approach, iterations, timestep) new{T}(grid, Vector{Grid{T}}(), bc, approach, iterations, timestep)
end end
end end
"""
createGrid(simulation::DynamicSimulation{T})::Int where {T}
Creates a new grid state based on the initial state and adds it to the simulation's grid vector.
# Arguments
- `simulation::DynamicSimulation{T}`: The dynamic simulation object.
# Returns
The index of the newly created grid in the simulation's grid vector.
"""
function createGrid(simulation::DynamicSimulation{T})::Int where {T} function createGrid(simulation::DynamicSimulation{T})::Int where {T}
new_grid = clone(simulation.grid) new_grid = clone(simulation.grid)
push!(simulation.grids, new_grid) push!(simulation.grids, new_grid)
return length(simulation.grids) return length(simulation.grids)
end end
function next(simulation::DynamicSimulation{T}) where {T} """
next(simulation::DynamicSimulation{T})::Nothing where {T}
Runs the simulation for all grid states in parallel.
# Arguments
- `simulation::DynamicSimulation{T}`: The dynamic simulation object.
# Returns
Nothing, but updates all grids in the simulation according to the defined approach and iterations.
"""
function next(simulation::DynamicSimulation{T})::Nothing where {T}
pmap(grid -> runSimulationForGrid(simulation, grid), simulation.grids) pmap(grid -> runSimulationForGrid(simulation, grid), simulation.grids)
end end
function printConcentrations(simulation::DynamicSimulation{T}, gridIndex::Int) where {T} """
printConcentrations(simulation::DynamicSimulation{T}, gridIndex::Int)::Nothing where {T}
Prints the concentrations of a specific grid state to the console.
# Arguments
- `simulation::DynamicSimulation{T}`: The dynamic simulation object.
- `gridIndex::Int`: Index of the grid in the simulation's grid vector.
# Returns
Nothing, but outputs the grid's concentration matrix to the console.
"""
function printConcentrations(simulation::DynamicSimulation{T}, gridIndex::Int)::Nothing where {T}
writeConcentrationsToCLI(simulation, simulation.grids[gridIndex]) writeConcentrationsToCLI(simulation, simulation.grids[gridIndex])
end end
function printConcentrationsCSV(simulation::DynamicSimulation{T}, gridIndex::Int) where {T} """
printConcentrationsCSV(simulation::DynamicSimulation{T}, gridIndex::Int)::Nothing where {T}
Writes the concentrations of a specific grid state to a CSV file.
# Arguments
- `simulation::DynamicSimulation{T}`: The dynamic simulation object.
- `gridIndex::Int`: Index of the grid in the simulation's grid vector.
# Returns
Nothing, but writes the grid's concentration matrix to a CSV file.
"""
function printConcentrationsCSV(simulation::DynamicSimulation{T}, gridIndex::Int)::Nothing where {T}
file = createCSVfile(simulation, simulation.grids[gridIndex]) file = createCSVfile(simulation, simulation.grids[gridIndex])
writeConcentrationsToCSV(file, simulation, simulation.grids[gridIndex]) writeConcentrationsToCSV(file, simulation, simulation.grids[gridIndex])
close(file) close(file)
end end
function runSimulationForGrid(simulation::DynamicSimulation{T}, grid::Grid{T}) where {T} function runSimulationForGrid(simulation::DynamicSimulation{T}, grid::Grid{T})::Nothing where {T}
if simulation.approach == BTCS 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 elseif simulation.approach == FTCS
@ -46,18 +108,69 @@ function runSimulationForGrid(simulation::DynamicSimulation{T}, grid::Grid{T}) w
end end
end end
"""
getConcentrations(simulation::DynamicSimulation{T}, gridIndex::Int)::Matrix{T} where {T}
Retrieves the concentrations of a specific grid state.
# Arguments
- `simulation::DynamicSimulation{T}`: The dynamic simulation object.
- `gridIndex::Int`: Index of the grid in the simulation's grid vector.
# 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]) getConcentrations(simulation.grids[gridIndex])
end end
function setConcentrations!(simulation::DynamicSimulation{T}, gridIndex::Int, concentrations::Matrix{T}) where {T} """
setConcentrations!(simulation::DynamicSimulation{T}, gridIndex::Int, concentrations::Matrix{T})::Nothing where {T}
Sets the concentrations of a specific grid state.
# Arguments
- `simulation::DynamicSimulation{T}`: The dynamic simulation object.
- `gridIndex::Int`: Index of the grid in the simulation's grid vector.
- `concentrations::Matrix{T}`: The new concentration matrix.
# Returns
Nothing, but updates the concentration matrix of the specified grid.
"""
function setConcentrations!(simulation::DynamicSimulation{T}, gridIndex::Int, concentrations::Matrix{T})::Nothing where {T}
setConcentrations!(simulation.grids[gridIndex], concentrations) setConcentrations!(simulation.grids[gridIndex], concentrations)
end end
function setAlphaX!(simulation::DynamicSimulation{T}, gridIndex::Int, alphaX::Matrix{T}) where {T} """
setAlphaX!(simulation::DynamicSimulation{T}, gridIndex::Int, alphaX::Matrix{T})::Nothing where {T}
Sets the alphaX matrix of a specific grid state.
# Arguments
- `simulation::DynamicSimulation{T}`: The dynamic simulation object.
- `gridIndex::Int`: Index of the grid in the simulation's grid vector.
- `alphaX::Matrix{T}`: The new alphaX matrix.
# Returns
Nothing, but updates the alphaX matrix of the specified grid.
"""
function setAlphaX!(simulation::DynamicSimulation{T}, gridIndex::Int, alphaX::Matrix{T})::Nothing where {T}
setAlphaX!(simulation.grids[gridIndex], alphaX) setAlphaX!(simulation.grids[gridIndex], alphaX)
end end
function setAlphaY!(simulation::DynamicSimulation{T}, gridIndex::Int, alphaY::Matrix{T}) where {T} """
setAlphaY!(simulation::DynamicSimulation{T}, gridIndex::Int, alphaY::Matrix{T})::Nothing where {T}
Sets the alphaY matrix of a specific grid state.
# Arguments
- `simulation::DynamicSimulation{T}`: The dynamic simulation object.
- `gridIndex::Int`: Index of the grid in the simulation's grid vector.
- `alphaY::Matrix{T}`: The new alphaY matrix.
# Returns
Nothing, but updates the alphaY matrix of the specified grid.
"""
function setAlphaY!(simulation::DynamicSimulation{T}, gridIndex::Int, alphaY::Matrix{T})::Nothing where {T}
setAlphaY!(simulation.grids[gridIndex], alphaY) setAlphaY!(simulation.grids[gridIndex], alphaY)
end end

View File

@ -5,6 +5,28 @@
using LinearAlgebra using LinearAlgebra
"""
struct Grid{T}
Represents a computational grid with concentration values and alpha coefficients.
It can be either 1D or 2D and includes dimensions, domain sizes, delta values, and matrices for concentrations and alpha coefficients.
# Fields
- `cols::Int`: Number of columns in the grid.
- `rows::Int`: Number of rows in the grid.
- `dim::Int`: Dimension of the grid (1 for 1D, 2 for 2D).
- `domainCol::T`, `domainRow::T`: Size of the grid's domain in column and row direction.
- `deltaCol::T`, `deltaRow::T`: Delta values for columns and rows.
- `concentrations::Ref{Matrix{T}}`: Reference to the matrix holding concentration values.
- `alphaX::Ref{Matrix{T}}`: Reference to the matrix of alpha coefficients in the X direction.
- `alphaY::Union{Ref{Matrix{T}},Nothing}`: Reference to the matrix of alpha coefficients in the Y direction (if applicable).
- `alphaX_t`, `alphaY_t`: Transposed matrices of alpha coefficients.
# Constructors
- `Grid(length, alphaX)` creates a 1D grid with the given length and alphaX matrix.
- `Grid(rows, cols, alphaX, alphaY)` creates a 2D grid with the given rows, columns, and alphaX and alphaY matrices.
- `Grid(rows, cols, dim, domainCol, domainRow, deltaCol, deltaRow, concentrations, alphaX, alphaY, alphaX_t, alphaY_t)` creates a grid with the given parameters.
"""
struct Grid{T} struct Grid{T}
cols::Int cols::Int
rows::Int rows::Int
@ -19,8 +41,7 @@ struct Grid{T}
alphaX_t::Union{Ref{Matrix{T}},Nothing} alphaX_t::Union{Ref{Matrix{T}},Nothing}
alphaY_t::Union{Ref{Matrix{T}},Nothing} alphaY_t::Union{Ref{Matrix{T}},Nothing}
# Constructor for 1D-Grid function Grid{T}(length::Int, alphaX::Matrix{T})::Grid{T} where {T}
function Grid{T}(length::Int, alphaX::Matrix{T}) where {T}
if length <= 3 if length <= 3
throw(ArgumentError("Given grid length too small. Must be greater than 3.")) throw(ArgumentError("Given grid length too small. Must be greater than 3."))
end end
@ -33,8 +54,7 @@ struct Grid{T}
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 end
# Constructor for 2D-Grid 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}) where {T}
if rows <= 3 || cols <= 3 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 end
@ -42,18 +62,29 @@ struct Grid{T}
throw(ArgumentError("Given matrices of alpha coefficients mismatch with Grid dimensions!")) throw(ArgumentError("Given matrices of alpha coefficients mismatch with Grid dimensions!"))
end end
# Precompute alphaX_t and alphaY_t
alphaX_t = alphaX' alphaX_t = alphaX'
alphaY_t = alphaY' 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 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}) where {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) new{T}(cols, rows, dim, domainCol, domainRow, deltaCol, deltaRow, concentrations, alphaX, alphaY, alphaX_t, alphaY_t)
end end
end end
"""
clone(grid::Grid{T})::Grid{T} where {T}
Creates a deep copy of the specified grid. The new grid will have identical properties
but will be a distinct object in memory.
# Arguments
- `grid::Grid{T}`: The grid to be cloned.
# Returns
A new `Grid{T}` object that is a deep copy of the input grid.
"""
function clone(grid::Grid{T})::Grid{T} where {T} function clone(grid::Grid{T})::Grid{T} where {T}
if grid.dim == 1 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)
@ -61,10 +92,33 @@ function clone(grid::Grid{T})::Grid{T} where {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[]))) 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 end
"""
getAlphaX(grid::Grid{T})::Matrix{T} where {T}
Retrieves the alpha coefficients matrix in the X direction from the specified grid.
# Arguments
- `grid::Grid{T}`: The grid from which to retrieve the alphaX matrix.
# Returns
The `alphaX` matrix of the grid.
"""
function getAlphaX(grid::Grid{T})::Matrix{T} where {T} function getAlphaX(grid::Grid{T})::Matrix{T} where {T}
grid.alphaX[] grid.alphaX[]
end end
"""
getAlphaY(grid::Grid{T})::Matrix{T} where {T}
Retrieves the alpha coefficients matrix in the Y direction from the specified grid.
Raises an error if the grid is 1D (no alphaY matrix).
# Arguments
- `grid::Grid{T}`: The grid from which to retrieve the alphaY matrix.
# Returns
The `alphaY` matrix of the grid, if applicable.
"""
function getAlphaY(grid::Grid{T})::Matrix{T} where {T} function getAlphaY(grid::Grid{T})::Matrix{T} where {T}
if grid.dim == 1 if grid.dim == 1
error("Grid is 1D, so there is no alphaY matrix!") error("Grid is 1D, so there is no alphaY matrix!")
@ -73,10 +127,33 @@ function getAlphaY(grid::Grid{T})::Matrix{T} where {T}
grid.alphaY[] grid.alphaY[]
end end
"""
getAlphaX_t(grid::Grid{T})::Matrix{T} where {T}
Retrieves the transposed alpha coefficients matrix in the X direction from the specified grid.
# Arguments
- `grid::Grid{T}`: The grid from which to retrieve the alphaX_t matrix.
# Returns
The transposed `alphaX_t` matrix of the grid.
"""
function getAlphaX_t(grid::Grid{T})::Matrix{T} where {T} function getAlphaX_t(grid::Grid{T})::Matrix{T} where {T}
grid.alphaX_t[] grid.alphaX_t[]
end end
"""
getAlphaY_t(grid::Grid{T})::Matrix{T} where {T}
Retrieves the transposed alpha coefficients matrix in the Y direction from the specified grid.
Raises an error if the grid is 1D (no alphaY_t matrix).
# Arguments
- `grid::Grid{T}`: The grid from which to retrieve the alphaY_t matrix.
# Returns
The transposed `alphaY_t` matrix of the grid, if applicable.
"""
function getAlphaY_t(grid::Grid{T})::Matrix{T} where {T} function getAlphaY_t(grid::Grid{T})::Matrix{T} where {T}
if grid.dim == 1 if grid.dim == 1
error("Grid is 1D, so there is no alphaY_t matrix!") error("Grid is 1D, so there is no alphaY_t matrix!")
@ -85,40 +162,133 @@ function getAlphaY_t(grid::Grid{T})::Matrix{T} where {T}
grid.alphaY_t[] grid.alphaY_t[]
end end
"""
getCols(grid::Grid{T})::Int where {T}
Retrieves the number of columns in the specified grid.
# Arguments
- `grid::Grid{T}`: The grid from which to retrieve the number of columns.
# Returns
The number of columns in the grid.
"""
function getCols(grid::Grid{T})::Int where {T} function getCols(grid::Grid{T})::Int where {T}
grid.cols grid.cols
end end
"""
getConcentrations(grid::Grid{T})::Matrix{T} where {T}
Retrieves the concentration matrix from the specified grid.
# Arguments
- `grid::Grid{T}`: The grid from which to retrieve the concentration matrix.
# Returns
The concentration matrix of the grid.
"""
function getConcentrations(grid::Grid{T})::Matrix{T} where {T} function getConcentrations(grid::Grid{T})::Matrix{T} where {T}
grid.concentrations[] grid.concentrations[]
end end
"""
getDeltaCol(grid::Grid{T})::T where {T}
Retrieves the delta value for columns in the specified grid.
# Arguments
- `grid::Grid{T}`: The grid from which to retrieve the delta column value.
# Returns
The delta value for columns in the grid.
"""
function getDeltaCol(grid::Grid{T})::T where {T} function getDeltaCol(grid::Grid{T})::T where {T}
grid.deltaCol grid.deltaCol
end end
"""
getDeltaRow(grid::Grid{T})::T where {T}
Retrieves the delta value for rows in the specified grid.
# Arguments
- `grid::Grid{T}`: The grid from which to retrieve the delta row value.
# Returns
The delta value for rows in the grid.
"""
function getDeltaRow(grid::Grid{T})::T where {T} function getDeltaRow(grid::Grid{T})::T where {T}
grid.deltaRow grid.deltaRow
end end
"""
getDim(grid::Grid{T})::Int where {T}
Retrieves the dimension (1D or 2D) of the specified grid.
# Arguments
- `grid::Grid{T}`: The grid from which to retrieve the dimension.
# Returns
The dimension of the grid (1 for 1D, 2 for 2D).
"""
function getDim(grid::Grid{T})::Int where {T} function getDim(grid::Grid{T})::Int where {T}
grid.dim grid.dim
end end
"""
getRows(grid::Grid{T})::Int where {T}
Retrieves the number of rows in the specified grid.
# Arguments
- `grid::Grid{T}`: The grid from which to retrieve the number of rows.
# Returns
The number of rows in the grid.
"""
function getRows(grid::Grid{T})::Int where {T} function getRows(grid::Grid{T})::Int where {T}
grid.rows grid.rows
end end
function setAlphaX!(grid::Grid{T}, new_alphaX::Matrix{T}) where {T} """
setAlphaX!(grid::Grid{T}, new_alphaX::Matrix{T})::Nothing where {T}
Sets the alpha coefficients matrix in the X direction for the specified grid.
Throws an error if the dimensions of the new matrix don't match the grid's dimensions.
# Arguments
- `grid::Grid{T}`: The grid for which to set the new alphaX matrix.
- `new_alphaX::Matrix{T}`: The new alpha coefficients matrix in the X direction.
# Returns
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[]) 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 end
grid.alphaX[] = new_alphaX grid.alphaX[] = new_alphaX
grid.alphaX_t[] = new_alphaX' grid.alphaX_t[] = new_alphaX'
return
end end
function setAlphaY!(grid::Grid{T}, new_alphaY::Matrix{T}) where {T} """
setAlphaY!(grid::Grid{T}, new_alphaY::Matrix{T})::Nothing where {T}
Sets the alpha coefficients matrix in the Y direction for the specified grid.
Throws an error if the grid is 1D or if the dimensions of the new matrix don't match the grid's dimensions.
# Arguments
- `grid::Grid{T}`: The grid for which to set the new alphaY matrix.
- `new_alphaY::Matrix{T}`: The new alpha coefficients matrix in the Y direction.
# Returns
Nothing, but modifies the alphaY matrix of the grid.
"""
function setAlphaY!(grid::Grid{T}, new_alphaY::Matrix{T})::Nothing where {T}
if grid.dim == 1 if grid.dim == 1
error("Grid is 1D, so there is no alphaY matrix!") error("Grid is 1D, so there is no alphaY matrix!")
end end
@ -128,12 +298,27 @@ function setAlphaY!(grid::Grid{T}, new_alphaY::Matrix{T}) where {T}
grid.alphaY[] = new_alphaY grid.alphaY[] = new_alphaY
grid.alphaY_t[] = new_alphaY' grid.alphaY_t[] = new_alphaY'
return
end end
function setConcentrations!(grid::Grid{T}, new_concentrations::Matrix{T}) where {T} """
setConcentrations!(grid::Grid{T}, new_concentrations::Matrix{T})::Nothing where {T}
Sets the concentration matrix for the specified grid.
Throws an error if the dimensions of the new matrix don't match the grid's dimensions.
# Arguments
- `grid::Grid{T}`: The grid for which to set the new concentrations matrix.
- `new_concentrations::Matrix{T}`: The new concentrations matrix.
# Returns
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[]) 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 end
grid.concentrations[] = new_concentrations grid.concentrations[] = new_concentrations
return
end end

View File

@ -4,6 +4,25 @@
# options. Simulation object also holds a predefined Grid and Boundary object. # options. Simulation object also holds a predefined Grid and Boundary object.
# Translated from C++'s Simulation.hpp. # Translated from C++'s Simulation.hpp.
"""
struct Simulation{T} <: AbstractSimulation{T}
Represents a simulation environment with a specified grid and boundary conditions.
Holds all necessary information for a simulation run including timestep, number
of iterations, output options, and simulation approach.
# Fields
- `grid::Grid{T}`: Grid object defining the simulation space.
- `bc::Boundary{T}`: Boundary conditions for the simulation.
- `approach::APPROACH`: Numerical approach for the simulation (e.g., BTCS, FTCS).
- `iterations::Int`: Number of iterations the simulation will run.
- `timestep::T`: Time step for each simulation iteration.
- `consoleOutput::CONSOLE_OUTPUT`: Option for console output level.
- `csvOutput::CSV_OUTPUT`: Option for CSV file output level.
# Constructor
- `Simulation(grid, bc, approach, iterations, timestep, consoleOutput, csvOutput)` creates a new simulation with specified parameters.
"""
struct Simulation{T} <: AbstractSimulation{T} struct Simulation{T} <: AbstractSimulation{T}
grid::Grid{T} grid::Grid{T}
bc::Boundary{T} bc::Boundary{T}
@ -15,14 +34,25 @@ struct Simulation{T} <: AbstractSimulation{T}
consoleOutput::CONSOLE_OUTPUT consoleOutput::CONSOLE_OUTPUT
csvOutput::CSV_OUTPUT csvOutput::CSV_OUTPUT
# Constructor
function Simulation(grid::Grid{T}, bc::Boundary{T}, approach::APPROACH=BTCS, iterations::Int=1, timestep::T=0.1, 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) where {T} 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) new{T}(grid, bc, approach, iterations, timestep, consoleOutput, csvOutput)
end end
end end
function run(simulation::Simulation{T}) where {T} """
run(simulation::Simulation{T})::Nothing where {T}
Executes the simulation based on the defined parameters within the `simulation` object.
Depending on the output settings, it may write data to console and/or a CSV file.
# Arguments
- `simulation::Simulation{T}`: The simulation object to be run.
# Returns
Nothing, but executes the simulation and handles outputs.
"""
function run(simulation::Simulation{T})::Nothing where {T}
file = nothing file = nothing
try try
if simulation.csvOutput > CSV_OUTPUT_OFF if simulation.csvOutput > CSV_OUTPUT_OFF
@ -54,7 +84,6 @@ function run(simulation::Simulation{T}) where {T}
if simulation.csvOutput >= CSV_OUTPUT_ON if simulation.csvOutput >= CSV_OUTPUT_ON
writeConcentrationsToCSV(file, simulation) writeConcentrationsToCSV(file, simulation)
end end
finally finally
if file !== nothing if file !== nothing
close(file) close(file)
@ -62,18 +91,66 @@ function run(simulation::Simulation{T}) where {T}
end end
end end
"""
setIterations(simulation::Simulation{T}, iterations::Int)::Simulation{T} where {T}
Sets the number of iterations for the given simulation.
# Arguments
- `simulation::Simulation{T}`: The simulation object.
- `iterations::Int`: The new number of iterations to be set.
# Returns
A new `Simulation` object with updated iterations.
"""
function setIterations(simulation::Simulation{T}, iterations::Int)::Simulation{T} where {T} 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 end
"""
setOutputConsole(simulation::Simulation{T}, consoleOutput::CONSOLE_OUTPUT)::Simulation{T} where {T}
Sets the console output level for the simulation.
# Arguments
- `simulation::Simulation{T}`: The simulation object.
- `consoleOutput::CONSOLE_OUTPUT`: The new console output level.
# Returns
A new `Simulation` object with updated console output setting.
"""
function setOutputConsole(simulation::Simulation{T}, consoleOutput::CONSOLE_OUTPUT)::Simulation{T} where {T} 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) return Simulation(simulation.grid, simulation.bc, simulation.approach, simulation.iterations, simulation.timestep, consoleOutput, simulation.csvOutput)
end end
"""
setOutputCSV(simulation::Simulation{T}, csvOutput::CSV_OUTPUT)::Simulation{T} where {T}
Sets the CSV output level for the simulation.
# Arguments
- `simulation::Simulation{T}`: The simulation object.
- `csvOutput::CSV_OUTPUT`: The new CSV output level.
# Returns
A new `Simulation` object with updated CSV output setting.
"""
function setOutputCSV(simulation::Simulation{T}, csvOutput::CSV_OUTPUT)::Simulation{T} where {T} 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) return Simulation(simulation.grid, simulation.bc, simulation.approach, simulation.iterations, simulation.timestep, simulation.consoleOutput, csvOutput)
end end
"""
setTimestep(simulation::Simulation{T}, timestep::T)::Simulation{T} where {T}
Sets the timestep for the simulation.
# Arguments
- `simulation::Simulation{T}`: The simulation object.
- `timestep::T`: The new timestep to be set.
# Returns
A new `Simulation` object with updated timestep.
"""
function setTimestep(simulation::Simulation{T}, timestep::T)::Simulation{T} where {T} 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 end

View File

@ -3,7 +3,7 @@ module TUG
include("Grid.jl") include("Grid.jl")
export Grid export Grid
export getAlphaX, getAlphaY, 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") include("Boundary.jl")
@ -12,7 +12,7 @@ export getType, getValue, setBoundarySideClosed!, setBoundarySideConstant!, getB
include("AbstractSimulation.jl") include("AbstractSimulation.jl")
export AbstractSimulation, 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") include("Simulation.jl")
@ -22,7 +22,7 @@ export run, setTimestep, setIterations, setOutputConsole, setOutputCSV
include("DynamicSimulation.jl") include("DynamicSimulation.jl")
export DynamicSimulation export DynamicSimulation
export createGrid, getConcentrations, setConcentrations!, setAlphaX!, setAlphaY!, next, printConcentrationsCSV, printConcentrations export createGrid, next, printConcentrations, printConcentrationsCSV, getConcentrations, setConcentrations!, setAlphaX!, setAlphaY!
include("Core/Utils.jl") include("Core/Utils.jl")
include("Core/BTCS.jl") include("Core/BTCS.jl")