0eb96ed0ad
Modified the createSolutionVector function to use matrix operations Additionally added a function to create the entire solution matrix. This function, whilst currently active, is better suited for GPU usage. [skip ci]
356 lines
15 KiB
Julia
356 lines
15 KiB
Julia
# BTCS.jl
|
|
# Implementation of heterogenous BTCS (backward time-centered space)
|
|
# solution of diffusion equation in 1D and 2D space using the
|
|
# alternating-direction implicit (ADI) method.
|
|
# Translated from C++'s BTCS.hpp.
|
|
|
|
using LinearAlgebra
|
|
using SparseArrays
|
|
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
|
|
|
|
|
|
function calcBoundaryCoeffClosed(alpha_center::T, alpha_side::T, sx::T) where {T}
|
|
alpha = calcAlphaIntercell(alpha_center, alpha_side)
|
|
centerCoeff = 1 + sx * alpha
|
|
sideCoeff = -sx * alpha
|
|
return (centerCoeff, sideCoeff)
|
|
end
|
|
|
|
function calcBoundaryCoeffConstant(alpha_center::T, alpha_side::T, sx::T) where {T}
|
|
alpha = calcAlphaIntercell(alpha_center, alpha_side)
|
|
centerCoeff = 1 + sx * (alpha + 2 * alpha_center)
|
|
sideCoeff = -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::Matrix{T}, alpha_right::Matrix{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
|
|
error("Undefined Boundary Condition Type on Left!")
|
|
end
|
|
|
|
# Precompute boundary condition type check for efficiency
|
|
bcTypeRight = getType(bcRight[rowIndex])
|
|
|
|
# Determine right side boundary coefficients based on boundary condition
|
|
centerCoeffBottom, leftCoeffBottom = if bcTypeRight == CONSTANT
|
|
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[rowIndex, :]; leftCoeffBottom]
|
|
d = [centerCoeffTop; 1 .+ sx .* (alpha_right[rowIndex, :] + alpha_left[rowIndex, :]); centerCoeffBottom]
|
|
du = [rightCoeffTop; -sx .* alpha_right[rowIndex, :]]
|
|
alpha_diagonal = Tridiagonal(dl, d, du)
|
|
|
|
return alpha_diagonal
|
|
end
|
|
|
|
|
|
function calcExplicitConcentrationsBoundaryClosed(conc_center::T, alpha_center::T, alpha_neighbor::T, sy::T) where {T}
|
|
alpha = calcAlphaIntercell(alpha_center, alpha_neighbor)
|
|
|
|
(sy * alpha + (1 - sy * alpha)) * conc_center
|
|
end
|
|
|
|
function calcExplicitConcentrationsBoundaryConstant(conc_center::T, conc_bc::T, alpha_center::T, alpha_neighbor::T, sy::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)
|
|
|
|
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) where {T}
|
|
(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
|
|
|
|
# creates a solution vector for next time step from the current state of concentrations
|
|
function createSolutionVector(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}}, length::Int, rowIndex::Int, sx::T, sy::T) where {T}
|
|
numRows = size(concentrations, 1)
|
|
sv = zeros(T, length)
|
|
|
|
# Inner rows
|
|
if rowIndex > 1 && rowIndex < numRows
|
|
alpha_neighbour_below = calcAlphaIntercell(alphaY[rowIndex, :], alphaY[rowIndex+1, :])
|
|
alpha_neighbour_above = calcAlphaIntercell(alphaY[rowIndex, :], alphaY[rowIndex-1, :])
|
|
|
|
# Compute sv with Array Operations
|
|
sv = sy .* alpha_neighbour_below .* concentrations[rowIndex+1, :] +
|
|
(1 .- sy .* (alpha_neighbour_below .+ alpha_neighbour_above)) .* concentrations[rowIndex, :] +
|
|
sy .* alpha_neighbour_above .* concentrations[rowIndex-1, :]
|
|
end
|
|
|
|
# First row
|
|
if rowIndex == 1
|
|
alpha_center = calcAlphaIntercell(alphaY[rowIndex, :], alphaY[rowIndex, :])
|
|
alpha_neighour_right = calcAlphaIntercell(alphaY[rowIndex, :], alphaY[rowIndex+1, :])
|
|
|
|
# Apply vectorized operations based on boundary condition
|
|
if getType(bcTop[1]) == CONSTANT
|
|
sv = sy .* alpha_neighour_right .* concentrations[rowIndex, :] +
|
|
(1 .- sy .* (alpha_center .+ 2 .* alphaY[rowIndex, :])) .* concentrations[rowIndex, :] +
|
|
sy .* alphaY[rowIndex, :] .* getValue(bcTop)
|
|
elseif getType(bcTop[1]) == CLOSED
|
|
sv = (sy .* alpha_neighour_right .+ (1 .- sy .* alpha_neighour_right)) .* concentrations[rowIndex, :]
|
|
else
|
|
error("Undefined Boundary Condition Type somewhere on Left or Top!")
|
|
end
|
|
end
|
|
|
|
# Last row
|
|
if rowIndex == numRows
|
|
alpha_center = calcAlphaIntercell(alphaY[rowIndex, :], alphaY[rowIndex, :])
|
|
alpha_neighbour_left = calcAlphaIntercell(alphaY[rowIndex, :], alphaY[rowIndex-1, :])
|
|
|
|
# Apply vectorized operations based on boundary condition
|
|
if getType(bcBottom[1]) == CONSTANT
|
|
sv = sy .* alpha_neighbour_left .* concentrations[rowIndex, :] +
|
|
(1 .- sy .* (alpha_center .+ 2 .* alphaY[rowIndex, :])) .* concentrations[rowIndex, :] +
|
|
sy .* alphaY[rowIndex, :] .* getValue(bcBottom)
|
|
elseif getType(bcBottom[1]) == CLOSED
|
|
sv = (sy .* alpha_neighbour_left .+ (1 .- sy .* alpha_neighbour_left)) .* concentrations[rowIndex, :]
|
|
else
|
|
error("Undefined Boundary Condition Type somewhere on Right or Bottom!")
|
|
end
|
|
end
|
|
|
|
# Conditions for the first and last columns
|
|
is_first_column = (1:length) .== 1
|
|
is_last_column = (1:length) .== length
|
|
|
|
# Apply operations conditionally
|
|
if getType(bcLeft[rowIndex]) == CONSTANT
|
|
sv .+= (2 * sx * alphaX[rowIndex, 1] * getValue(bcLeft[rowIndex])) .* is_first_column
|
|
end
|
|
|
|
if getType(bcRight[rowIndex]) == CONSTANT
|
|
sv .+= (2 * sx * alphaX[rowIndex, end] * getValue(bcRight[rowIndex])) .* is_last_column
|
|
end
|
|
|
|
return sv
|
|
end
|
|
|
|
function createSolutionMatrix(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}}, sx::T, sy::T) where {T}
|
|
numRows, numCols = size(concentrations)
|
|
solutionMatrix = zeros(T, numRows, numCols)
|
|
|
|
# Vectorized Precomputation of Alphas
|
|
alpha_neighbour_below = calcAlphaIntercell(alphaY[1:end-1, :], alphaY[2:end, :])
|
|
alpha_neighbour_above = calcAlphaIntercell(alphaY[2:end, :], alphaY[1:end-1, :])
|
|
alpha_center = calcAlphaIntercell(alphaY, alphaY)
|
|
|
|
# Compute solutionMatrix for inner rows
|
|
solutionMatrix[2:end-1, :] = sy .* alpha_neighbour_below[1:end-1, :] .* concentrations[3:end, :] +
|
|
(1 .- sy .* (alpha_neighbour_below[1:end-1, :] .+ alpha_neighbour_above[2:end, :])) .* concentrations[2:end-1, :] +
|
|
sy .* alpha_neighbour_above[2:end, :] .* concentrations[1:end-2, :]
|
|
|
|
|
|
# Apply boundary conditions for first row
|
|
if getType(bcTop[1]) == CONSTANT
|
|
solutionMatrix[1, :] = sy .* alpha_center[1, :] .* concentrations[1, :] +
|
|
(1 .- sy .* (alpha_center[1, :] .+ 2 .* alphaY[1, :])) .* concentrations[1, :] +
|
|
sy .* alphaY[1, :] .* getValue(bcTop)
|
|
elseif getType(bcTop[1]) == CLOSED
|
|
solutionMatrix[1, :] = (sy .* alpha_center[1, :] .+ (1 .- sy .* alpha_center[1, :])) .* concentrations[1, :]
|
|
end
|
|
|
|
# Apply boundary conditions for last row
|
|
if getType(bcBottom[1]) == CONSTANT
|
|
solutionMatrix[end, :] = sy .* alpha_center[end, :] .* concentrations[end, :] +
|
|
(1 .- sy .* (alpha_center[end, :] .+ 2 .* alphaY[end, :])) .* concentrations[end, :] +
|
|
sy .* alphaY[end, :] .* getValue(bcBottom)
|
|
elseif getType(bcBottom[1]) == CLOSED
|
|
solutionMatrix[end, :] = (sy .* alpha_center[end, :] .+ (1 .- sy .* alpha_center[end, :])) .* concentrations[end, :]
|
|
end
|
|
|
|
# Apply boundary conditions for first and last columns
|
|
if getType(bcLeft[1]) == CONSTANT
|
|
solutionMatrix[:, 1] .+= 2 * sx * alphaX[:, 1] .* getValue(bcLeft)
|
|
end
|
|
if getType(bcRight[1]) == CONSTANT
|
|
solutionMatrix[:, end] .+= 2 * sx * alphaX[:, end] .* getValue(bcRight)
|
|
end
|
|
|
|
return solutionMatrix
|
|
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}
|
|
sx = timestep / (getDeltaCol(grid) * getDeltaCol(grid))
|
|
|
|
alpha = getAlphaX(grid)
|
|
bcLeft = getBoundarySide(bc, LEFT)
|
|
bcRight = getBoundarySide(bc, RIGHT)
|
|
length = getCols(grid)
|
|
|
|
concentrations::Matrix{T} = getConcentrations(grid)
|
|
|
|
A::Tridiagonal{T} = createCoeffMatrix(alpha, alpha_left, alpha_right, bcLeft, bcRight, length, 1, sx)
|
|
b = concentrations[1, :]
|
|
|
|
if getType(bcLeft[1]) == CONSTANT
|
|
b[1] += 2 * sx * alpha[1, 1] * getValue(bcLeft[1])
|
|
end
|
|
if getType(bcRight[1]) == CONSTANT
|
|
b[end] += 2 * sx * alpha[1, end] * getValue(bcRight[1])
|
|
end
|
|
|
|
concentrations[1, :] = A \ b
|
|
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}
|
|
rows = getRows(grid)
|
|
cols = getCols(grid)
|
|
sx = timestep / (2 * getDeltaCol(grid) * getDeltaCol(grid))
|
|
sy = timestep / (2 * getDeltaRow(grid) * getDeltaRow(grid))
|
|
|
|
alphaX = getAlphaX(grid)
|
|
alphaY = getAlphaY(grid)
|
|
|
|
alphaX_t = getAlphaX_t(grid)
|
|
alphaY_t = getAlphaY_t(grid)
|
|
|
|
concentrations = getConcentrations(grid)
|
|
concentrations_intermediate = similar(concentrations)
|
|
concentrations_t_task = Threads.@spawn copy(transpose(concentrations))
|
|
|
|
bcLeft = getBoundarySide(bc, LEFT)
|
|
bcRight = getBoundarySide(bc, RIGHT)
|
|
bcTop = getBoundarySide(bc, TOP)
|
|
bcBottom = getBoundarySide(bc, BOTTOM)
|
|
|
|
|
|
# Thread Based Computation:
|
|
# Threads.@threads for i = 1:rows
|
|
# A::Tridiagonal{T} = createCoeffMatrix(alphaX, alphaX_left[i, :], alphaX_right[i, :], bcLeft, bcRight, cols, i, sx)
|
|
# b = createSolutionVector(concentrations, alphaX, alphaY, bcLeft, bcRight, bcTop, bcBottom, cols, i, sx, sy)
|
|
|
|
# concentrations_intermediate[i, :] = A \ b
|
|
# end
|
|
|
|
# Compute solution vectors for all rows
|
|
b_matrix = createSolutionMatrix(concentrations, alphaX, alphaY, bcLeft, bcRight, bcTop, bcBottom, sx, sy)
|
|
|
|
# Process each row for the coefficient matrix and solve the linear system
|
|
Threads.@threads for i = 1:rows
|
|
A::Tridiagonal{T} = createCoeffMatrix(alphaX, alphaX_left, alphaX_right, bcLeft, bcRight, cols, i, sx)
|
|
|
|
# Extract the solution vector for the current row from b_matrix
|
|
b = b_matrix[i, :]
|
|
|
|
concentrations_intermediate[i, :] = A \ b
|
|
end
|
|
|
|
|
|
|
|
concentrations_intermediate = copy(transpose(concentrations_intermediate))
|
|
concentrations_t = fetch(concentrations_t_task)
|
|
|
|
|
|
# Swap alphas, boundary conditions and sx/sy for column-wise calculation
|
|
|
|
# Thread Based Computation:
|
|
# Threads.@threads for i = 1:cols
|
|
# A::Tridiagonal{T} = createCoeffMatrix(alphaY_t, alphaY_t_left[i, :], alphaY_t_right[i, :], bcTop, bcBottom, rows, i, sy)
|
|
# b = createSolutionVector(concentrations_intermediate, alphaY_t, alphaX_t, bcTop, bcBottom, bcLeft, bcRight, rows, i, sy, sx)
|
|
|
|
# concentrations_t[i, :] = A \ b
|
|
# end
|
|
|
|
# Compute solution vectors for all rows
|
|
b_matrix = createSolutionMatrix(concentrations_intermediate, alphaY_t, alphaX_t, bcTop, bcBottom, bcLeft, bcRight, sy, sx)
|
|
|
|
# Process each row for the coefficient matrix and solve the linear system
|
|
Threads.@threads for i = 1:cols
|
|
A::Tridiagonal{T} = createCoeffMatrix(alphaY_t, alphaY_t_left, alphaY_t_right, bcTop, bcBottom, rows, i, sy)
|
|
|
|
# Extract the solution vector for the current row from b_matrix
|
|
b = b_matrix[i, :]
|
|
|
|
concentrations_t[i, :] = A \ b
|
|
end
|
|
|
|
concentrations = copy(transpose(concentrations_t))
|
|
|
|
setConcentrations!(grid, concentrations)
|
|
end
|
|
|
|
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 = fetch(alpha_left_task)
|
|
alpha_right = fetch(alpha_right_task)
|
|
|
|
for _ in 1:(iterations)
|
|
BTCS_1D(grid, bc, alpha_left, alpha_right, timestep)
|
|
|
|
stepCallback()
|
|
end
|
|
elseif getDim(grid) == 2
|
|
rows = getRows(grid)
|
|
cols = getCols(grid)
|
|
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 = 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)
|
|
|
|
stepCallback()
|
|
end
|
|
else
|
|
error("Error: Only 1- and 2-dimensional grids are defined!")
|
|
end
|
|
|
|
end
|