diff --git a/julia/TUG/src/Core/BTCS.jl b/julia/TUG/src/Core/BTCS.jl index 6f99ab2..ecbca66 100644 --- a/julia/TUG/src/Core/BTCS.jl +++ b/julia/TUG/src/Core/BTCS.jl @@ -28,28 +28,20 @@ function calcBoundaryCoeff( end function createCoeffMatrix( - alpha::Matrix{T}, - alpha_left::Vector{T}, - alpha_right::Vector{T}, + alpha::SubArray{T}, + alpha_left::SubArray{T}, + alpha_right::SubArray{T}, bcLeft::Vector{BoundaryElement{T}}, bcRight::Vector{BoundaryElement{T}}, numCols::Int, rowIndex::Int, sx::T, )::Tridiagonal{T} where {T} - centerCoeffTop, rightCoeffTop = calcBoundaryCoeff( - alpha[rowIndex, 1], - alpha[rowIndex, 2], - sx, - getType(bcLeft[rowIndex]), - ) + centerCoeffTop, rightCoeffTop = + calcBoundaryCoeff(alpha[1], alpha[2], sx, getType(bcLeft[rowIndex])) - centerCoeffBottom, leftCoeffBottom = calcBoundaryCoeff( - alpha[rowIndex, numCols], - alpha[rowIndex, numCols-1], - sx, - getType(bcRight[rowIndex]), - ) + centerCoeffBottom, leftCoeffBottom = + calcBoundaryCoeff(alpha[numCols], alpha[numCols-1], sx, getType(bcRight[rowIndex])) dl = [-sx .* alpha_left; leftCoeffBottom] d = [centerCoeffTop; 1 .+ sx .* (alpha_right + alpha_left); centerCoeffBottom] @@ -95,25 +87,24 @@ end function writeSolutionVector!( sv::Vector{T}, - concentrations::Matrix{T}, + concentrations_t::Matrix{T}, alphaX::Matrix{T}, alphaY::Matrix{T}, - bcLeft::Vector{BoundaryElement{T}}, - bcRight::Vector{BoundaryElement{T}}, + bcLeft::BoundaryElement{T}, + bcRight::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) + numRows = size(concentrations_t, 2) if rowIndex == 1 - @inbounds for i = 1:length + @inbounds for i in eachindex(sv) if getType(bcTop[i]) == CONSTANT sv[i] = calcExplicitConcentrationsBoundaryConstant( - concentrations[rowIndex, i], + concentrations_t[i, rowIndex], getValue(bcTop[i]), alphaY[rowIndex, i], alphaY[rowIndex+1, i], @@ -121,7 +112,7 @@ function writeSolutionVector!( ) elseif getType(bcTop[i]) == CLOSED sv[i] = calcExplicitConcentrationsBoundaryClosed( - concentrations[rowIndex, i], + concentrations_t[i, rowIndex], alphaY[rowIndex, i], alphaY[rowIndex+1, i], sy, @@ -133,25 +124,26 @@ function writeSolutionVector!( end if rowIndex > 1 && rowIndex < numRows - @inbounds for i = 1:length + @inbounds for i in eachindex(sv) 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]) + sv[i] = - sy * alpha_here_below * concentrations[rowIndex+1, i] + + sy * alpha_here_below * concentrations_t[i, rowIndex+1] + (1 - sy * (alpha_here_below + alpha_here_above)) * - concentrations[rowIndex, i] + - sy * alpha_here_above * concentrations[rowIndex-1, i] + concentrations_t[i, rowIndex] + + sy * alpha_here_above * concentrations_t[i, rowIndex-1] end end if rowIndex == numRows - @inbounds for i = 1:length + @inbounds for i in eachindex(sv) if getType(bcBottom[i]) == CONSTANT sv[i] = calcExplicitConcentrationsBoundaryConstant( - concentrations[rowIndex, i], + concentrations_t[i, rowIndex], getValue(bcBottom[i]), alphaY[rowIndex, i], alphaY[rowIndex-1, i], @@ -159,7 +151,7 @@ function writeSolutionVector!( ) elseif getType(bcBottom[i]) == CLOSED sv[i] = calcExplicitConcentrationsBoundaryClosed( - concentrations[rowIndex, i], + concentrations_t[i, rowIndex], alphaY[rowIndex, i], alphaY[rowIndex-1, i], sy, @@ -170,12 +162,12 @@ function writeSolutionVector!( end end - if getType(bcLeft[rowIndex]) == CONSTANT - sv[1] += 2 * sx * alphaX[rowIndex, 1] * getValue(bcLeft[rowIndex]) + if getType(bcLeft) == CONSTANT + sv[1] += 2 * sx * alphaX[rowIndex, 1] * getValue(bcLeft) end - if getType(bcRight[rowIndex]) == CONSTANT - sv[end] += 2 * sx * alphaX[rowIndex, end] * getValue(bcRight[rowIndex]) + if getType(bcRight) == CONSTANT + sv[end] += 2 * sx * alphaX[rowIndex, end] * getValue(bcRight) end end @@ -196,9 +188,9 @@ function BTCS_1D( concentrations::Matrix{T} = getConcentrations(grid) A::Tridiagonal{T} = createCoeffMatrix( - alpha, - alpha_left[1, :], - alpha_right[1, :], + view(alpha, 1, :), + view(alpha_left, 1, :), + view(alpha_right, 1, :), bcLeft, bcRight, length, @@ -239,20 +231,19 @@ function BTCS_2D( concentrations = getConcentrations(grid) concentrations_intermediate = similar(concentrations) - concentrations_t_task = Threads.@spawn copy(transpose(concentrations)) + concentrations_t = copy(transpose(concentrations)) bcLeft = getBoundarySide(bc, LEFT) bcRight = getBoundarySide(bc, RIGHT) bcTop = getBoundarySide(bc, TOP) bcBottom = getBoundarySide(bc, BOTTOM) - localBs = [zeros(T, cols) for _ = 1:Threads.nthreads()] Threads.@threads for i = 1:rows - localB = localBs[Threads.threadid()] + localB = zeros(T, cols) A::Tridiagonal{T} = createCoeffMatrix( - alphaX, - alphaX_left[i, :], - alphaX_right[i, :], + view(alphaX, i, :), + view(alphaX_left, i, :), + view(alphaX_right, i, :), bcLeft, bcRight, cols, @@ -261,11 +252,11 @@ function BTCS_2D( ) writeSolutionVector!( localB, - concentrations, + concentrations_t, alphaX, alphaY, - bcLeft, - bcRight, + bcLeft[i], + bcRight[i], bcTop, bcBottom, i, @@ -276,17 +267,13 @@ function BTCS_2D( concentrations_intermediate[i, :] = A \ localB end - concentrations_intermediate = copy(transpose(concentrations_intermediate)) - concentrations_t = fetch(concentrations_t_task) - # Swap alphas, boundary conditions and sx/sy for column-wise calculation - localBs = [zeros(T, rows) for _ = 1:Threads.nthreads()] Threads.@threads for i = 1:cols - localB = localBs[Threads.threadid()] + localB = zeros(T, cols) A::Tridiagonal{T} = createCoeffMatrix( - alphaY_t, - alphaY_t_left[i, :], - alphaY_t_right[i, :], + view(alphaY_t, i, :), + view(alphaY_t_left, i, :), + view(alphaY_t_right, i, :), bcTop, bcBottom, rows, @@ -298,8 +285,8 @@ function BTCS_2D( concentrations_intermediate, alphaY_t, alphaX_t, - bcTop, - bcBottom, + bcTop[i], + bcBottom[i], bcLeft, bcRight, i, diff --git a/julia/TUG/src/Core/FTCS.jl b/julia/TUG/src/Core/FTCS.jl index 93cd77a..1c02eb4 100644 --- a/julia/TUG/src/Core/FTCS.jl +++ b/julia/TUG/src/Core/FTCS.jl @@ -127,8 +127,8 @@ function FTCS_1D(grid::Grid{T}, bc::Boundary{T}, timestep::T) where {T} colMax = getCols(grid) sx = timestep / (getDeltaCol(grid)^2) - concentrations = getConcentrations(grid) alphaX = getAlphaX(grid) + concentrations = getConcentrations(grid) concentrations_t1 = copy(concentrations) row = 1 @@ -184,26 +184,6 @@ function FTCS_2D(grid::Grid{T}, bc::Boundary{T}, timestep::T) where {T} concentrations_t1 = copy(concentrations) 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], - ) - end - concentrations_t1[row, 1] += sx * calcHorizontalChangeLeftBoundary( getBoundaryElementType(bc, LEFT, row), @@ -242,6 +222,26 @@ function FTCS_2D(grid::Grid{T}, bc::Boundary{T}, timestep::T) where {T} end Threads.@threads for col = 2:colMax-1 + for row = 2:rowMax-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], + ) + end + concentrations_t1[1, col] += sy * calcVerticalChangeTopBoundary( getBoundaryElementType(bc, TOP, col), diff --git a/julia/TUG/src/Grid.jl b/julia/TUG/src/Grid.jl index f7074f7..92347e8 100644 --- a/julia/TUG/src/Grid.jl +++ b/julia/TUG/src/Grid.jl @@ -35,11 +35,11 @@ struct Grid{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} + concentrations::Matrix{T} + alphaX::Matrix{T} + alphaY::Union{Matrix{T},Nothing} + alphaX_t::Matrix{T} + alphaY_t::Union{Matrix{T},Nothing} function Grid{T}(length::Int, alphaX::Matrix{T})::Grid{T} where {T} if length <= 3 @@ -63,7 +63,7 @@ struct Grid{T} 0, T(1), 0, - Ref(fill(T(0), 1, length)), + fill(T(0), 1, length), alphaX, nothing, alphaX_t, @@ -103,7 +103,7 @@ struct Grid{T} T(rows), T(1), T(1), - Ref(fill(T(0), rows, cols)), + fill(T(0), rows, cols), alphaX, alphaY, alphaX_t, @@ -119,11 +119,11 @@ struct Grid{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}, + concentrations::Matrix{T}, + alphaX::Matrix{T}, + alphaY::Union{Matrix{T},Nothing}, + alphaX_t::Matrix{T}, + alphaY_t::Union{Matrix{T},Nothing}, )::Grid{T} where {T} new{T}( cols, @@ -164,10 +164,10 @@ function clone(grid::Grid{T})::Grid{T} where {T} grid.domainRow, grid.deltaCol, grid.deltaRow, - Ref(copy(grid.concentrations[])), - Ref(copy(grid.alphaX[])), + copy(grid.concentrations), + copy(grid.alphaX), nothing, - Ref(copy(grid.alphaX_t[])), + copy(grid.alphaX_t), nothing, ) end @@ -179,11 +179,11 @@ function clone(grid::Grid{T})::Grid{T} where {T} 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[])), + copy(grid.concentrations), + copy(grid.alphaX), + copy(grid.alphaY), + copy(grid.alphaX_t), + copy(grid.alphaY_t), ) end @@ -199,7 +199,7 @@ Retrieves the alpha coefficients matrix in the X direction from the specified gr The `alphaX` matrix of the grid. """ function getAlphaX(grid::Grid{T})::Matrix{T} where {T} - grid.alphaX[] + grid.alphaX end """ @@ -219,7 +219,7 @@ function getAlphaY(grid::Grid{T})::Matrix{T} where {T} error("Grid is 1D, so there is no alphaY matrix!") end - grid.alphaY[] + grid.alphaY end """ @@ -234,7 +234,7 @@ Retrieves the transposed alpha coefficients matrix in the X direction from the s The transposed `alphaX_t` matrix of the grid. """ function getAlphaX_t(grid::Grid{T})::Matrix{T} where {T} - grid.alphaX_t[] + grid.alphaX_t end """ @@ -254,7 +254,7 @@ function getAlphaY_t(grid::Grid{T})::Matrix{T} where {T} error("Grid is 1D, so there is no alphaY_t matrix!") end - grid.alphaY_t[] + grid.alphaY_t end """ @@ -284,7 +284,7 @@ Retrieves the concentration matrix from the specified grid. The concentration matrix of the grid. """ function getConcentrations(grid::Grid{T})::Matrix{T} where {T} - grid.concentrations[] + grid.concentrations end """ @@ -361,7 +361,7 @@ Throws an error if the dimensions of the new matrix don't match the grid's dimen 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!", @@ -369,8 +369,8 @@ function setAlphaX!(grid::Grid{T}, new_alphaX::Matrix{T})::Nothing where {T} ) end - grid.alphaX[] = new_alphaX - grid.alphaX_t[] = new_alphaX' + grid.alphaX .= new_alphaX + grid.alphaX_t .= new_alphaX' return end @@ -391,7 +391,7 @@ function setAlphaY!(grid::Grid{T}, new_alphaY::Matrix{T})::Nothing where {T} if grid.dim == 1 error("Grid is 1D, so there is no alphaY matrix!") end - if size(new_alphaY) != size(grid.alphaY[]) + if size(new_alphaY) != size(grid.alphaY) throw( ArgumentError( "Given matrix of alpha coefficients mismatch with Grid dimensions!", @@ -399,8 +399,8 @@ function setAlphaY!(grid::Grid{T}, new_alphaY::Matrix{T})::Nothing where {T} ) end - grid.alphaY[] = new_alphaY - grid.alphaY_t[] = new_alphaY' + grid.alphaY .= new_alphaY + grid.alphaY_t .= new_alphaY' return end @@ -418,12 +418,12 @@ Throws an error if the dimensions of the new matrix don't match the grid's dimen 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!"), ) end - grid.concentrations[] = new_concentrations + grid.concentrations .= new_concentrations return end