# 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")

function calcAlphaIntercell(alpha1::T, alpha2::T) where {T}
    return 2 / ((1 / alpha1) + (1 / alpha2))
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

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

# 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
    @inbounds for i in 2:(numCols-1)
        alpha_left_here = calcAlphaIntercell(alpha[rowIndex, i-1], alpha[rowIndex, i])
        alpha_here_right = alpha[rowIndex, i-1] == alpha[rowIndex, i+1] ? alpha_left_here : calcAlphaIntercell(alpha[rowIndex, i], alpha[rowIndex, i+1]) # calcAlphaIntercell is symmetric, so we can use it for both directions

        cm[i, i-1] = -sx * alpha_left_here
        cm[i, i] = 1 + sx * (alpha_here_right + alpha_left_here)
        cm[i, i+1] = -sx * alpha_here_right
    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}
    alpha = calcAlphaIntercell(alpha_center, alpha_neighbor)
    sy * alpha * conc_center + (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 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
        @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
            sv[i] = sy * alpha_here_below * concentrations[rowIndex+1, i] +
                    (1 - sy * (alpha_here_below + alpha_here_above)) * concentrations[rowIndex, i] +
                    sy * alpha_here_above * concentrations[rowIndex-1, i]
        end
    end

    # First row
    if rowIndex == 1
        @inbounds 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
        @inbounds 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 = getAlphaX(grid)

    bcLeft = getBoundarySide(bc, LEFT)
    bcRight = getBoundarySide(bc, RIGHT)

    concentrations = grid.concentrations[]
    rowIndex = 1
    A = createCoeffMatrix(alpha, bcLeft, bcRight, length, rowIndex, sx)
    @inbounds 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)

    @inbounds 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)

    concentrations_t1 = zeros(T, rowMax, colMax)
    row_t1 = Vector{T}(undef, colMax)

    alphaX = getAlphaX(grid)
    alphaY = getAlphaY(grid)

    bcLeft = getBoundarySide(bc, LEFT)
    bcRight = getBoundarySide(bc, RIGHT)
    bcTop = getBoundarySide(bc, TOP)
    bcBottom = getBoundarySide(bc, BOTTOM)

    concentrations = grid.concentrations[]

    @inbounds 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 = getAlphaX_t(grid)
    alphaY = getAlphaY_t(grid)

    @inbounds 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

function BTCS_step(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