diff --git a/julia/tug/Boundary.jl b/julia/tug/Boundary.jl index d4e87ca..88810e3 100644 --- a/julia/tug/Boundary.jl +++ b/julia/tug/Boundary.jl @@ -77,6 +77,28 @@ struct Boundary{T} end end +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.")) + end + if index < 1 || index > (side in (LEFT, RIGHT) ? boundary.rows : boundary.cols) + throw(ArgumentError("Index out of bounds!")) + end + + getType(boundary.boundaries[Int(side)][index]) +end + +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.")) + end + if index < 1 || index > (side in (LEFT, RIGHT) ? boundary.rows : boundary.cols) + throw(ArgumentError("Index out of bounds!")) + end + + getValue(boundary.boundaries[Int(side)][index]) +end + 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.")) diff --git a/julia/tug/Core/BTCS.jl b/julia/tug/Core/BTCS.jl index 38e63d6..16b9058 100644 --- a/julia/tug/Core/BTCS.jl +++ b/julia/tug/Core/BTCS.jl @@ -10,19 +10,7 @@ using Base.Threads include("../Boundary.jl") include("../Grid.jl") - -function calcAlphaIntercell(alpha1::T, alpha2::T) where {T} - 2 / ((1 / alpha1) + (1 / alpha2)) -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} - 2 ./ ((1 ./ alpha1) .+ (1 ./ alpha2)) -end - +include("Utils.jl") function calcBoundaryCoeffClosed(alpha_center::T, alpha_side::T, sx::T) where {T} alpha = calcAlphaIntercell(alpha_center, alpha_side) diff --git a/julia/tug/Core/FTCS.jl b/julia/tug/Core/FTCS.jl new file mode 100644 index 0000000..8db17be --- /dev/null +++ b/julia/tug/Core/FTCS.jl @@ -0,0 +1,258 @@ +# FTCS.jl +# Implementation of heterogenous FTCS (forward time-centered space) +# solution of diffusion equation in 1D and 2D space. +# Translated from C++'s FTCS.hpp. + +include("../Boundary.jl") +include("../Grid.jl") +include("Utils.jl") + +function calcHorizontalChange(grid::Grid{T}, row::Int, col::Int) where {T} + return calcAlphaIntercell(getAlphaX(grid)[row, col+1], + getAlphaX(grid)[row, col]) * + getConcentrations(grid)[row, col+1] - + (calcAlphaIntercell(getAlphaX(grid)[row, col+1], + getAlphaX(grid)[row, col]) + + calcAlphaIntercell(getAlphaX(grid)[row, col-1], + getAlphaX(grid)[row, col])) * + getConcentrations(grid)[row, col] + + calcAlphaIntercell(getAlphaX(grid)[row, col-1], + getAlphaX(grid)[row, col]) * + getConcentrations(grid)[row, col-1] +end + +function calcHorizontalChangeLeftBoundary(grid::Grid{T}, bc::Boundary{T}, row::Int, col::Int) where {T} + if getBoundaryElementType(bc, LEFT, row) == CONSTANT + return calcHorizontalChangeLeftBoundaryConstant(grid, bc, row, col) + elseif getBoundaryElementType(bc, LEFT, row) == CLOSED + return calcHorizontalChangeLeftBoundaryClosed(grid, row, col) + else + error("Undefined Boundary Condition Type!") + end +end + +function calcHorizontalChangeLeftBoundaryConstant(grid::Grid{T}, bc::Boundary{T}, row::Int, col::Int) where {T} + return calcAlphaIntercell(getAlphaX(grid)[row, col+1], + getAlphaX(grid)[row, col]) * + getConcentrations(grid)[row, col+1] - + (calcAlphaIntercell(getAlphaX(grid)[row, col+1], + getAlphaX(grid)[row, col]) + + 2 * getAlphaX(grid)[row, col]) * + getConcentrations(grid)[row, col] + + 2 * getAlphaX(grid)[row, col] * + getBoundaryElementValue(bc, LEFT, row) +end + +function calcHorizontalChangeLeftBoundaryClosed(grid::Grid{T}, row::Int, col::Int) where {T} + return calcAlphaIntercell(getAlphaX(grid)[row, col+1], + getAlphaX(grid)[row, col]) * + (getConcentrations(grid)[row, col+1] - + getConcentrations(grid)[row, col]) +end + + +function calcHorizontalChangeRightBoundary(grid::Grid{T}, bc::Boundary{T}, row::Int, col::Int) where {T} + if getBoundaryElementType(bc, RIGHT, row) == CONSTANT + return calcHorizontalChangeRightBoundaryConstant(grid, bc, row, col) + elseif getBoundaryElementType(bc, RIGHT, row) == CLOSED + return calcHorizontalChangeRightBoundaryClosed(grid, row, col) + else + error("Undefined Boundary Condition Type!") + end +end + +function calcHorizontalChangeRightBoundaryConstant(grid::Grid{T}, bc::Boundary{T}, row::Int, col::Int) where {T} + return 2 * getAlphaX(grid)[row, col] * + getBoundaryElementValue(bc, RIGHT, row) - + (calcAlphaIntercell(getAlphaX(grid)[row, col-1], + getAlphaX(grid)[row, col]) + + 2 * getAlphaX(grid)[row, col]) * + getConcentrations(grid)[row, col] + + calcAlphaIntercell(getAlphaX(grid)[row, col-1], + getAlphaX(grid)[row, col]) * + getConcentrations(grid)[row, col-1] +end + +function calcHorizontalChangeRightBoundaryClosed(grid::Grid{T}, row::Int, col::Int) where {T} + return -(calcAlphaIntercell(getAlphaX(grid)[row, col-1], + getAlphaX(grid)[row, col]) * + (getConcentrations(grid)[row, col] - + getConcentrations(grid)[row, col-1])) +end + + +function calcVerticalChange(grid::Grid{T}, row::Int, col::Int) where {T} + return calcAlphaIntercell(getAlphaY(grid)[row+1, col], + getAlphaY(grid)[row, col]) * + getConcentrations(grid)[row+1, col] - + (calcAlphaIntercell(getAlphaY(grid)[row+1, col], + getAlphaY(grid)[row, col]) + + calcAlphaIntercell(getAlphaY(grid)[row-1, col], + getAlphaY(grid)[row, col])) * + getConcentrations(grid)[row, col] + + calcAlphaIntercell(getAlphaY(grid)[row-1, col], + getAlphaY(grid)[row, col]) * + getConcentrations(grid)[row-1, col] +end + +function calcVerticalChangeTopBoundary(grid::Grid{T}, bc::Boundary{T}, row::Int, col::Int) where {T} + if getBoundaryElementType(bc, TOP, col) == CONSTANT + return calcVerticalChangeTopBoundaryConstant(grid, bc, row, col) + elseif getBoundaryElementType(bc, TOP, col) == CLOSED + return calcVerticalChangeTopBoundaryClosed(grid, row, col) + else + error("Undefined Boundary Condition Type!") + end +end + +function calcVerticalChangeTopBoundaryConstant(grid::Grid{T}, bc::Boundary{T}, row::Int, col::Int) where {T} + return calcAlphaIntercell(getAlphaY(grid)[row+1, col], + getAlphaY(grid)[row, col]) * + getConcentrations(grid)[row+1, col] - + (calcAlphaIntercell(getAlphaY(grid)[row+1, col], + getAlphaY(grid)[row, col]) + + 2 * getAlphaY(grid)[row, col]) * + getConcentrations(grid)[row, col] + + 2 * getAlphaY(grid)[row, col] * + getBoundaryElementValue(bc, TOP, col) +end + +function calcVerticalChangeTopBoundaryClosed(grid::Grid{T}, row::Int, col::Int) where {T} + return calcAlphaIntercell(getAlphaY(grid)[row+1, col], + getAlphaY(grid)[row, col]) * + (getConcentrations(grid)[row+1, col] - + getConcentrations(grid)[row, col]) +end + +function calcVerticalChangeBottomBoundary(grid::Grid{T}, bc::Boundary{T}, row::Int, col::Int) where {T} + if getBoundaryElementType(bc, BOTTOM, col) == CONSTANT + return calcVerticalChangeBottomBoundaryConstant(grid, bc, row, col) + elseif getBoundaryElementType(bc, BOTTOM, col) == CLOSED + return calcVerticalChangeBottomBoundaryClosed(grid, row, col) + else + error("Undefined Boundary Condition Type!") + end +end + +function calcVerticalChangeBottomBoundaryConstant(grid::Grid{T}, bc::Boundary{T}, row::Int, col::Int) where {T} + return 2 * getAlphaY(grid)[row, col] * + getBoundaryElementValue(bc, BOTTOM, col) - + (calcAlphaIntercell(getAlphaY(grid)[row, col], + getAlphaY(grid)[row-1, col]) + + 2 * getAlphaY(grid)[row, col]) * + getConcentrations(grid)[row, col] + + calcAlphaIntercell(getAlphaY(grid)[row, col], + getAlphaY(grid)[row-1, col]) * + getConcentrations(grid)[row-1, col] +end + +function calcVerticalChangeBottomBoundaryClosed(grid::Grid{T}, row::Int, col::Int) where {T} + return -(calcAlphaIntercell(getAlphaY(grid)[row, col], + getAlphaY(grid)[row-1, col]) * + (getConcentrations(grid)[row, col] - + getConcentrations(grid)[row-1, col])) +end + + + +function FTCS_1D(grid::Grid{T}, bc::Boundary{T}, timestep::T) where {T} + colMax = getCols(grid) + sx = timestep / getDeltaCol(grid) * getDeltaCol(grid) + + # matrix for concentrations at time t+1 + concentrations_t1 = zeros(1, colMax) + + # only one row in 1D case -> row constant at index 1 + row = 1 + + # inner cells + # independent of boundary condition type + for col = 2:colMax-1 + concentrations_t1[row, col] = getConcentrations(grid)[row, col] + + sx * calcHorizontalChange(grid, row, col) + end + + # left boundary; hold column constant at index 1 + concentrations_t1[row, 1] = getConcentrations(grid)[row, 1] + + sx * calcHorizontalChangeLeftBoundary(grid, bc, row, 1) + + # right boundary; hold column constant at max index + concentrations_t1[row, colMax] = getConcentrations(grid)[row, colMax] + + sx * calcHorizontalChangeRightBoundary(grid, bc, row, colMax) + + # overwrite obsolete concentrations + setConcentrations!(grid, concentrations_t1) +end + +function FTCS_2D(grid::Grid{T}, bc::Boundary{T}, timestep::T) where {T} + rowMax = getRows(grid) + colMax = getCols(grid) + sx = timestep / (getDeltaCol(grid)^2) + sy = timestep / (getDeltaRow(grid)^2) + + concentrations = getConcentrations(grid) + concentrations_t1 = copy(concentrations) + + # inner cells + # these are independent of the boundary condition type + for row = 2:rowMax-1 + # inner cells + for col = 2:colMax-1 + concentrations_t1[row, col] += sy * calcVerticalChange(grid, row, col) + + sx * calcHorizontalChange(grid, row, col) + end + + # left boundary without corners + concentrations_t1[row, 1] += sx * calcHorizontalChangeLeftBoundary(grid, bc, row, 1) + + sy * calcVerticalChange(grid, row, 1) + # right boundary without corners + concentrations_t1[row, colMax] += sx * calcHorizontalChangeRightBoundary(grid, bc, row, colMax) + + sy * calcVerticalChange(grid, row, colMax) + end + + # top / bottom without corners / looping over columns + for col = 2:colMax-1 + # top + concentrations_t1[1, col] += sy * calcVerticalChangeTopBoundary(grid, bc, 1, col) + + sx * calcHorizontalChange(grid, 1, col) + # bottom + concentrations_t1[rowMax, col] += sy * calcVerticalChangeBottomBoundary(grid, bc, rowMax, col) + + sx * calcHorizontalChange(grid, rowMax, col) + end + + # corner top left + concentrations_t1[1, 1] += sx * calcHorizontalChangeLeftBoundary(grid, bc, 1, 1) + + sy * calcVerticalChangeTopBoundary(grid, bc, 1, 1) + + # corner top right + concentrations_t1[1, colMax] += sx * calcHorizontalChangeRightBoundary(grid, bc, 1, colMax) + + sy * calcVerticalChangeTopBoundary(grid, bc, 1, colMax) + + # corner bottom left + concentrations_t1[rowMax, 1] += sx * calcHorizontalChangeLeftBoundary(grid, bc, rowMax, 1) + + sy * calcVerticalChangeBottomBoundary(grid, bc, rowMax, 1) + + # corner bottom right + concentrations_t1[rowMax, colMax] += sx * calcHorizontalChangeRightBoundary(grid, bc, rowMax, colMax) + + sy * calcVerticalChangeBottomBoundary(grid, bc, rowMax, colMax) + + setConcentrations!(grid, concentrations_t1) +end + +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) + + stepCallback() + end + elseif getDim(grid) == 2 + for _ = 1:iterations + FTCS_2D(grid, bc, timestep) + + stepCallback() + end + else + error("FTCS only implemented for 1D and 2D grids") + end +end diff --git a/julia/tug/Core/Utils.jl b/julia/tug/Core/Utils.jl new file mode 100644 index 0000000..03578b5 --- /dev/null +++ b/julia/tug/Core/Utils.jl @@ -0,0 +1,12 @@ + +function calcAlphaIntercell(alpha1::T, alpha2::T) where {T} + 2 / ((1 / alpha1) + (1 / alpha2)) +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} + 2 ./ ((1 ./ alpha1) .+ (1 ./ alpha2)) +end diff --git a/julia/tug/Simulation.jl b/julia/tug/Simulation.jl index aa8a03b..9f0785f 100644 --- a/julia/tug/Simulation.jl +++ b/julia/tug/Simulation.jl @@ -9,8 +9,9 @@ using Printf include("Grid.jl") include("Boundary.jl") include("Core/BTCS.jl") +include("Core/FTCS.jl") -@enum APPROACH BTCS +@enum APPROACH BTCS FTCS @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 @@ -105,6 +106,8 @@ function run(simulation::Simulation{T}) where {T} if simulation.approach == BTCS runBTCS(simulation.grid, simulation.bc, simulation.timestep, simulation.iterations, simulationStepCallback) + elseif simulation.approach == FTCS + runFTCS(simulation.grid, simulation.bc, simulation.timestep, simulation.iterations, simulationStepCallback) else error("Undefined approach!") end