From 28cb5416f5d217e523a95823a8f53c66519e88af Mon Sep 17 00:00:00 2001 From: nebmit <76664673+nebmit@users.noreply.github.com> Date: Tue, 5 Dec 2023 18:26:00 +0100 Subject: [PATCH] feat: added option to modify singular boundary elements --- julia/TUG/src/Boundary.jl | 34 +++++++++++++++++++++ julia/TUG/src/Core/BTCS.jl | 10 +------ julia/TUG/src/Core/FTCS.jl | 61 ++++++++++---------------------------- julia/TUG/src/TUG.jl | 7 +++-- 4 files changed, 55 insertions(+), 57 deletions(-) diff --git a/julia/TUG/src/Boundary.jl b/julia/TUG/src/Boundary.jl index b3c3909..9abd377 100644 --- a/julia/TUG/src/Boundary.jl +++ b/julia/TUG/src/Boundary.jl @@ -289,3 +289,37 @@ function setBoundarySideConstant!( boundary.boundaries[Int(side)] = [BoundaryElement{T}(value) for _ = 1:n] end + +""" + setBoundarySideElement(boundary::Boundary{T}, side::SIDE, index::Int, be::BoundaryElement{T})::BoundaryElement{T} where {T} + +Sets 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. +- `be::BoundaryElement{T}`: The boundary element. + +# Returns +The boundary element at the specified side and index. +""" +function setBoundarySideElement( + boundary::Boundary{T}, + side::SIDE, + index::Int, + be::BoundaryElement{T}, +)::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.", + ), + ) + end + if index < 1 || index > (side in (LEFT, RIGHT) ? boundary.rows : boundary.cols) + throw(ArgumentError("Index out of bounds!")) + end + + boundary.boundaries[Int(side)][index] = be +end diff --git a/julia/TUG/src/Core/BTCS.jl b/julia/TUG/src/Core/BTCS.jl index aee378f..6f99ab2 100644 --- a/julia/TUG/src/Core/BTCS.jl +++ b/julia/TUG/src/Core/BTCS.jl @@ -37,7 +37,6 @@ function createCoeffMatrix( 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], @@ -45,7 +44,6 @@ function createCoeffMatrix( getType(bcLeft[rowIndex]), ) - # Determine right side boundary coefficients based on boundary condition centerCoeffBottom, leftCoeffBottom = calcBoundaryCoeff( alpha[rowIndex, numCols], alpha[rowIndex, numCols-1], @@ -95,7 +93,6 @@ function calcExplicitConcentrationsBoundaryClosed( 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}, @@ -112,7 +109,6 @@ function writeSolutionVector!( numRows = size(concentrations, 1) length = size(sv, 1) - # First row if rowIndex == 1 @inbounds for i = 1:length if getType(bcTop[i]) == CONSTANT @@ -136,14 +132,13 @@ function writeSolutionVector!( 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 + calcAlphaIntercell(alphaY[rowIndex-1, i], alphaY[rowIndex, i]) sv[i] = sy * alpha_here_below * concentrations[rowIndex+1, i] + (1 - sy * (alpha_here_below + alpha_here_above)) * @@ -152,7 +147,6 @@ function writeSolutionVector!( end end - # Last row if rowIndex == numRows @inbounds for i = 1:length if getType(bcBottom[i]) == CONSTANT @@ -176,12 +170,10 @@ function writeSolutionVector!( end end - # First column - additional fixed concentration change from perpendicular dimension in constant BC case if getType(bcLeft[rowIndex]) == CONSTANT sv[1] += 2 * sx * alphaX[rowIndex, 1] * getValue(bcLeft[rowIndex]) end - # Last column - additional fixed concentration change from perpendicular dimension in constant BC case if getType(bcRight[rowIndex]) == CONSTANT sv[end] += 2 * sx * alphaX[rowIndex, end] * getValue(bcRight[rowIndex]) end diff --git a/julia/TUG/src/Core/FTCS.jl b/julia/TUG/src/Core/FTCS.jl index 15cf188..93cd77a 100644 --- a/julia/TUG/src/Core/FTCS.jl +++ b/julia/TUG/src/Core/FTCS.jl @@ -14,7 +14,6 @@ function calcHorizontalChange( concentration_prev::T, concentration_current::T, ) where {T} - intercellAlpha_next = calcAlphaIntercell(alphaX_next, alphaX_current) intercellAlpha_prev = calcAlphaIntercell(alphaX_prev, alphaX_current) @@ -134,7 +133,6 @@ function FTCS_1D(grid::Grid{T}, bc::Boundary{T}, timestep::T) where {T} row = 1 - # inner cells for col = 2:colMax-1 concentrations_t1[row, col] += sx * calcHorizontalChange( @@ -147,7 +145,6 @@ function FTCS_1D(grid::Grid{T}, bc::Boundary{T}, timestep::T) where {T} ) end - # left boundary leftBoundaryType = getBoundaryElementType(bc, LEFT, row) leftBoundaryValue = getBoundaryElementValue(bc, LEFT, row) concentrations_t1[row, 1] += @@ -160,7 +157,6 @@ function FTCS_1D(grid::Grid{T}, bc::Boundary{T}, timestep::T) where {T} leftBoundaryValue, ) - # right boundary rightBoundaryType = getBoundaryElementType(bc, RIGHT, row) rightBoundaryValue = getBoundaryElementValue(bc, RIGHT, row) concentrations_t1[row, colMax] += @@ -187,12 +183,6 @@ function FTCS_2D(grid::Grid{T}, bc::Boundary{T}, timestep::T) where {T} concentrations = getConcentrations(grid) concentrations_t1 = copy(concentrations) - # Currently boundary sides can only be set for the whole side, not for each cell individually - leftBoundaryType = getBoundaryElementType(bc, LEFT, 1) - rightBoundaryType = getBoundaryElementType(bc, RIGHT, 1) - topBoundaryType = getBoundaryElementType(bc, TOP, 1) - bottomBoundaryType = getBoundaryElementType(bc, BOTTOM, 1) - Threads.@threads for row = 2:rowMax-1 for col = 2:colMax-1 concentrations_t1[row, col] += @@ -214,17 +204,14 @@ function FTCS_2D(grid::Grid{T}, bc::Boundary{T}, timestep::T) where {T} ) end - # Boundary conditions for each row - # Left boundary without corners - leftBoundaryValue = getBoundaryElementValue(bc, LEFT, row) concentrations_t1[row, 1] += sx * calcHorizontalChangeLeftBoundary( - leftBoundaryType, + getBoundaryElementType(bc, LEFT, row), alphaX[row, 2], alphaX[row, 1], concentrations[row, 2], concentrations[row, 1], - leftBoundaryValue, + getBoundaryElementValue(bc, LEFT, row), ) + sy * calcVerticalChange( alphaY[row+1, 1], @@ -235,16 +222,14 @@ function FTCS_2D(grid::Grid{T}, bc::Boundary{T}, timestep::T) where {T} concentrations[row, 1], ) - # Right boundary without corners - rightBoundaryValue = getBoundaryElementValue(bc, RIGHT, row) concentrations_t1[row, colMax] += sx * calcHorizontalChangeRightBoundary( - rightBoundaryType, + getBoundaryElementType(bc, RIGHT, row), alphaX[row, colMax-1], alphaX[row, colMax], concentrations[row, colMax-1], concentrations[row, colMax], - rightBoundaryValue, + getBoundaryElementValue(bc, RIGHT, row), ) + sy * calcVerticalChange( alphaY[row+1, colMax], @@ -256,18 +241,15 @@ function FTCS_2D(grid::Grid{T}, bc::Boundary{T}, timestep::T) where {T} ) 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, + getBoundaryElementType(bc, TOP, col), alphaY[2, col], alphaY[1, col], concentrations[2, col], concentrations[1, col], - topBoundaryValue, + getBoundaryElementValue(bc, TOP, col), ) + sx * calcHorizontalChange( alphaX[1, col+1], @@ -278,16 +260,14 @@ function FTCS_2D(grid::Grid{T}, bc::Boundary{T}, timestep::T) where {T} concentrations[1, col], ) - # Bottom boundary - bottomBoundaryValue = getBoundaryElementValue(bc, BOTTOM, col) concentrations_t1[rowMax, col] += sy * calcVerticalChangeBottomBoundary( - bottomBoundaryType, + getBoundaryElementType(bc, BOTTOM, col), alphaY[rowMax, col], alphaY[rowMax-1, col], concentrations[rowMax, col], concentrations[rowMax-1, col], - bottomBoundaryValue, + getBoundaryElementValue(bc, BOTTOM, col), ) + sx * calcHorizontalChange( alphaX[rowMax, col+1], @@ -299,9 +279,6 @@ function FTCS_2D(grid::Grid{T}, bc::Boundary{T}, timestep::T) where {T} ) end - # Handle corners - # Top left corner - topBoundaryValue = getBoundaryElementValue(bc, TOP, 1) concentrations_t1[1, 1] += sx * calcHorizontalChange( alphaX[1, 2], @@ -312,16 +289,14 @@ function FTCS_2D(grid::Grid{T}, bc::Boundary{T}, timestep::T) where {T} concentrations[1, 1], ) + sy * calcVerticalChangeTopBoundary( - topBoundaryType, + getBoundaryElementType(bc, TOP, 1), alphaY[2, 1], alphaY[1, 1], concentrations[2, 1], concentrations[1, 1], - topBoundaryValue, + getBoundaryElementValue(bc, TOP, 1), ) - # Top right corner - topBoundaryValue = getBoundaryElementValue(bc, TOP, colMax) concentrations_t1[1, colMax] += sx * calcHorizontalChange( alphaX[1, colMax-1], @@ -332,16 +307,14 @@ function FTCS_2D(grid::Grid{T}, bc::Boundary{T}, timestep::T) where {T} concentrations[1, colMax], ) + sy * calcVerticalChangeTopBoundary( - topBoundaryType, + getBoundaryElementType(bc, TOP, colMax), alphaY[2, colMax], alphaY[1, colMax], concentrations[2, colMax], concentrations[1, colMax], - topBoundaryValue, + getBoundaryElementValue(bc, TOP, colMax), ) - # Bottom left corner - bottomBoundaryValue = getBoundaryElementValue(bc, BOTTOM, 1) concentrations_t1[rowMax, 1] += sx * calcHorizontalChange( alphaX[rowMax, 2], @@ -352,16 +325,14 @@ function FTCS_2D(grid::Grid{T}, bc::Boundary{T}, timestep::T) where {T} concentrations[rowMax, 1], ) + sy * calcVerticalChangeBottomBoundary( - bottomBoundaryType, + getBoundaryElementType(bc, BOTTOM, 1), alphaY[rowMax, 1], alphaY[rowMax-1, 1], concentrations[rowMax, 1], concentrations[rowMax-1, 1], - bottomBoundaryValue, + getBoundaryElementValue(bc, BOTTOM, 1), ) - # Bottom right corner - bottomBoundaryValue = getBoundaryElementValue(bc, BOTTOM, colMax) concentrations_t1[rowMax, colMax] += sx * calcHorizontalChange( alphaX[rowMax, colMax-1], @@ -372,12 +343,12 @@ function FTCS_2D(grid::Grid{T}, bc::Boundary{T}, timestep::T) where {T} concentrations[rowMax, colMax], ) + sy * calcVerticalChangeBottomBoundary( - bottomBoundaryType, + getBoundaryElementType(bc, BOTTOM, colMax), alphaY[rowMax, colMax], alphaY[rowMax-1, colMax], concentrations[rowMax, colMax], concentrations[rowMax-1, colMax], - bottomBoundaryValue, + getBoundaryElementValue(bc, BOTTOM, colMax), ) setConcentrations!(grid, concentrations_t1) diff --git a/julia/TUG/src/TUG.jl b/julia/TUG/src/TUG.jl index 069e646..094f6a3 100644 --- a/julia/TUG/src/TUG.jl +++ b/julia/TUG/src/TUG.jl @@ -24,11 +24,12 @@ export Boundary, BoundaryElement, TYPE, SIDE, BC_TYPE_CLOSED, BC_TYPE_CONSTANT, LEFT, RIGHT, TOP, BOTTOM export getType, getValue, - setBoundarySideClosed!, - setBoundarySideConstant!, getBoundarySide, getBoundaryElementType, - getBoundaryElementValue + getBoundaryElementValue, + setBoundarySideClosed!, + setBoundarySideConstant!, + setBoundarySideElement include("AbstractSimulation.jl")