245 lines
9.4 KiB
Julia
245 lines
9.4 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.
|
|
|
|
using LinearAlgebra
|
|
using SparseArrays
|
|
|
|
include("../Boundary.jl")
|
|
include("../Grid.jl")
|
|
|
|
# Helper functions and types
|
|
function calcAlphaIntercell(alpha1::T, alpha2::T, useHarmonic::Bool=true) where {T}
|
|
if useHarmonic
|
|
return 2 / ((1 / alpha1) + (1 / alpha2))
|
|
else
|
|
return 0.5 * (alpha1 + alpha2)
|
|
end
|
|
end
|
|
|
|
# calculates coefficient for boundary in constant case
|
|
function calcBoundaryCoeffConstant(alpha_center::T, alpha_side::T, sx::T) where {T}
|
|
centerCoeff = 1 + sx * (calcAlphaIntercell(alpha_center, alpha_side) + 2 * alpha_center)
|
|
sideCoeff = -sx * calcAlphaIntercell(alpha_center, alpha_side)
|
|
return (centerCoeff, sideCoeff)
|
|
end
|
|
|
|
# calculates coefficient for boundary in closed case
|
|
function calcBoundaryCoeffClosed(alpha_center::T, alpha_side::T, sx::T) where {T}
|
|
centerCoeff = 1 + sx * calcAlphaIntercell(alpha_center, alpha_side)
|
|
sideCoeff = -sx * calcAlphaIntercell(alpha_center, alpha_side)
|
|
return (centerCoeff, sideCoeff)
|
|
end
|
|
|
|
# creates coefficient matrix for next time step from alphas in x-direction
|
|
function createCoeffMatrix(alpha::Matrix{T}, bcLeft::Vector{BoundaryElement{T}}, bcRight::Vector{BoundaryElement{T}}, numCols::Int, rowIndex::Int, sx::T) where {T}
|
|
numCols = max(numCols, 2)
|
|
cm = spzeros(T, numCols, numCols)
|
|
|
|
# left column
|
|
if getType(bcLeft[rowIndex]) == CONSTANT
|
|
centerCoeffTop, rightCoeffTop = calcBoundaryCoeffConstant(alpha[rowIndex, 1], alpha[rowIndex, 2], sx)
|
|
cm[1, 1] = centerCoeffTop
|
|
cm[1, 2] = rightCoeffTop
|
|
elseif getType(bcLeft[rowIndex]) == CLOSED
|
|
centerCoeffTop, rightCoeffTop = calcBoundaryCoeffClosed(alpha[rowIndex, 1], alpha[rowIndex, 2], sx)
|
|
cm[1, 1] = centerCoeffTop
|
|
cm[1, 2] = rightCoeffTop
|
|
else
|
|
error("Undefined Boundary Condition Type somewhere on Left or Top!")
|
|
end
|
|
|
|
# inner columns
|
|
for i in 2:(numCols-1)
|
|
cm[i, i-1] = -sx * calcAlphaIntercell(alpha[rowIndex, i-1], alpha[rowIndex, i])
|
|
cm[i, i] = 1 + sx * (calcAlphaIntercell(alpha[rowIndex, i], alpha[rowIndex, i+1]) + calcAlphaIntercell(alpha[rowIndex, i-1], alpha[rowIndex, i]))
|
|
cm[i, i+1] = -sx * calcAlphaIntercell(alpha[rowIndex, i], alpha[rowIndex, i+1])
|
|
end
|
|
|
|
# right column
|
|
if getType(bcRight[rowIndex]) == CONSTANT
|
|
centerCoeffBottom, leftCoeffBottom = calcBoundaryCoeffConstant(alpha[rowIndex, numCols], alpha[rowIndex, numCols-1], sx)
|
|
cm[numCols, numCols-1] = leftCoeffBottom
|
|
cm[numCols, numCols] = centerCoeffBottom
|
|
elseif getType(bcRight[rowIndex]) == CLOSED
|
|
centerCoeffBottom, leftCoeffBottom = calcBoundaryCoeffClosed(alpha[rowIndex, numCols], alpha[rowIndex, numCols-1], sx)
|
|
cm[numCols, numCols-1] = leftCoeffBottom
|
|
cm[numCols, numCols] = centerCoeffBottom
|
|
else
|
|
error("Undefined Boundary Condition Type somewhere on Right or Bottom!")
|
|
end
|
|
|
|
return cm
|
|
end
|
|
|
|
|
|
|
|
|
|
function calcExplicitConcentrationsBoundaryClosed(conc_center::T, alpha_center::T, alpha_neighbor::T, sy::T) where {T}
|
|
sy * calcAlphaIntercell(alpha_center, alpha_neighbor) * conc_center +
|
|
(1 - sy * calcAlphaIntercell(alpha_center, alpha_neighbor)) * conc_center
|
|
end
|
|
|
|
|
|
function calcExplicitConcentrationsBoundaryConstant(conc_center::T, conc_bc::T, alpha_center::T, alpha_neighbor::T, sy::T) where {T}
|
|
sy * calcAlphaIntercell(alpha_center, alpha_neighbor) * conc_center +
|
|
(1 - sy * (calcAlphaIntercell(alpha_center, alpha_center) + 2 * alpha_center)) * conc_center +
|
|
sy * alpha_center * conc_bc
|
|
end
|
|
|
|
|
|
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 = Vector{T}(undef, length)
|
|
|
|
# Inner rows
|
|
if rowIndex > 1 && rowIndex < numRows
|
|
for i = 1:length
|
|
sv[i] = sy * calcAlphaIntercell(alphaY[rowIndex, i], alphaY[rowIndex+1, i]) * concentrations[rowIndex+1, i] +
|
|
(1 - sy * (calcAlphaIntercell(alphaY[rowIndex, i], alphaY[rowIndex+1, i]) + calcAlphaIntercell(alphaY[rowIndex-1, i], alphaY[rowIndex, i]))) * concentrations[rowIndex, i] +
|
|
sy * calcAlphaIntercell(alphaY[rowIndex-1, i], alphaY[rowIndex, i]) * concentrations[rowIndex-1, i]
|
|
end
|
|
end
|
|
|
|
# First row
|
|
if rowIndex == 1
|
|
for i = 1:length
|
|
if getType(bcTop[i]) == CONSTANT
|
|
sv[i] = calcExplicitConcentrationsBoundaryConstant(concentrations[rowIndex, i], getValue(bcTop[i]), alphaY[rowIndex, i], alphaY[rowIndex+1, i], sy)
|
|
elseif getType(bcTop[i]) == CLOSED
|
|
sv[i] = calcExplicitConcentrationsBoundaryClosed(concentrations[rowIndex, i], alphaY[rowIndex, i], alphaY[rowIndex+1, i], sy)
|
|
else
|
|
error("Undefined Boundary Condition Type somewhere on Left or Top!")
|
|
end
|
|
end
|
|
end
|
|
|
|
# Last row
|
|
if rowIndex == numRows
|
|
for i = 1:length
|
|
if getType(bcBottom[i]) == CONSTANT
|
|
sv[i] = calcExplicitConcentrationsBoundaryConstant(concentrations[rowIndex, i], getValue(bcBottom[i]), alphaY[rowIndex, i], alphaY[rowIndex-1, i], sy)
|
|
elseif getType(bcBottom[i]) == CLOSED
|
|
sv[i] = calcExplicitConcentrationsBoundaryClosed(concentrations[rowIndex, i], alphaY[rowIndex, i], alphaY[rowIndex-1, i], sy)
|
|
else
|
|
error("Undefined Boundary Condition Type somewhere on Right or Bottom!")
|
|
end
|
|
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
|
|
|
|
return sv
|
|
end
|
|
|
|
|
|
# solver for linear equation system; A corresponds to coefficient matrix, b to the solution vector
|
|
function LinearAlgebraAlgorithm(A::SparseMatrixCSC{T}, b::Vector{T}) where {T}
|
|
return A \ b
|
|
end
|
|
|
|
# BTCS solution for 1D grid
|
|
function BTCS_1D(grid::Grid{T}, bc::Boundary{T}, timestep::T, solverFunc::Function) where {T}
|
|
length = grid.cols
|
|
sx = timestep / (grid.deltaCol * grid.deltaCol)
|
|
|
|
b = Vector{T}(undef, length)
|
|
|
|
alpha = grid.alphaX[]
|
|
|
|
bcLeft = getBoundarySide(bc, LEFT)
|
|
bcRight = getBoundarySide(bc, RIGHT)
|
|
|
|
concentrations = grid.concentrations[]
|
|
rowIndex = 1
|
|
A = createCoeffMatrix(alpha, bcLeft, bcRight, length, rowIndex, sx)
|
|
for i in 1:length
|
|
b[i] = concentrations[1, i]
|
|
end
|
|
if getType(getBoundarySide(bc, LEFT)[1]) == CONSTANT
|
|
b[1] += 2 * sx * alpha[1, 1] * bcLeft[1].value
|
|
end
|
|
if getType(getBoundarySide(bc, RIGHT)[1]) == CONSTANT
|
|
b[length] += 2 * sx * alpha[1, length] * bcRight[1].value
|
|
end
|
|
|
|
concentrations_t1 = solverFunc(A, b)
|
|
|
|
for j in 1:length
|
|
concentrations[1, j] = concentrations_t1[j]
|
|
end
|
|
|
|
setConcentrations!(grid, concentrations)
|
|
end
|
|
|
|
# BTCS solution for 2D grid
|
|
function BTCS_2D(grid::Grid{T}, bc::Boundary{T}, timestep::T, solverFunc::Function, numThreads::Int) where {T}
|
|
rowMax = grid.rows
|
|
colMax = grid.cols
|
|
sx = timestep / (2 * grid.deltaCol * grid.deltaCol)
|
|
sy = timestep / (2 * grid.deltaRow * grid.deltaRow)
|
|
|
|
A = spzeros(T, rowMax, rowMax)
|
|
concentrations_t1 = zeros(T, rowMax, colMax)
|
|
row_t1 = Vector{T}(undef, colMax)
|
|
|
|
alphaX = grid.alphaX[]
|
|
alphaY = grid.alphaY[]
|
|
|
|
bcLeft = getBoundarySide(bc, LEFT)
|
|
bcRight = getBoundarySide(bc, RIGHT)
|
|
bcTop = getBoundarySide(bc, TOP)
|
|
bcBottom = getBoundarySide(bc, BOTTOM)
|
|
|
|
concentrations = grid.concentrations[]
|
|
|
|
for i = 1:rowMax
|
|
A = createCoeffMatrix(alphaX, bcLeft, bcRight, colMax, i, sx)
|
|
b = createSolutionVector(concentrations, alphaX, alphaY, bcLeft, bcRight, bcTop, bcBottom, colMax, i, sx, sy)
|
|
|
|
row_t1 = solverFunc(A, b)
|
|
|
|
concentrations_t1[i, :] = row_t1
|
|
end
|
|
|
|
concentrations_t1 = copy(transpose(concentrations_t1))
|
|
concentrations = copy(transpose(concentrations))
|
|
alphaX = copy(transpose(alphaX))
|
|
alphaY = copy(transpose(alphaY))
|
|
|
|
for i = 1:colMax
|
|
# Swap alphas, boundary conditions and sx/sy for column-wise calculation
|
|
A = createCoeffMatrix(alphaY, bcTop, bcBottom, rowMax, i, sy)
|
|
b = createSolutionVector(concentrations_t1, alphaY, alphaX, bcTop, bcBottom, bcLeft, bcRight, rowMax, i, sy, sx)
|
|
|
|
row_t1 = solverFunc(A, b)
|
|
|
|
concentrations[i, :] = row_t1
|
|
end
|
|
|
|
concentrations = copy(transpose(concentrations))
|
|
|
|
setConcentrations!(grid, concentrations)
|
|
end
|
|
|
|
|
|
# Entry point for EigenLU solver; differentiate between 1D and 2D grid
|
|
function BTCS_LU(grid::Grid{T}, bc::Boundary{T}, timestep::T, numThreads::Int=1) where {T}
|
|
if grid.dim == 1
|
|
BTCS_1D(grid, bc, timestep, LinearAlgebraAlgorithm)
|
|
elseif grid.dim == 2
|
|
BTCS_2D(grid, bc, timestep, LinearAlgebraAlgorithm, numThreads)
|
|
else
|
|
error("Error: Only 1- and 2-dimensional grids are defined!")
|
|
end
|
|
end
|