perf: removed solution matrix computation approach
[skip ci]
This commit is contained in:
nebmit
parent
0eb96ed0ad
commit
f538658e82
@ -39,7 +39,7 @@ function calcBoundaryCoeffConstant(alpha_center::T, alpha_side::T, sx::T) where
|
||||
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}
|
||||
function createCoeffMatrix(alpha::Matrix{T}, alpha_left::Vector{T}, alpha_right::Vector{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])
|
||||
|
||||
@ -64,9 +64,9 @@ function createCoeffMatrix(alpha::Matrix{T}, alpha_left::Matrix{T}, alpha_right:
|
||||
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, :]]
|
||||
dl = [-sx .* alpha_left; leftCoeffBottom]
|
||||
d = [centerCoeffTop; 1 .+ sx .* (alpha_right + alpha_left); centerCoeffBottom]
|
||||
du = [rightCoeffTop; -sx .* alpha_right]
|
||||
alpha_diagonal = Tridiagonal(dl, d, du)
|
||||
|
||||
return alpha_diagonal
|
||||
@ -99,114 +99,58 @@ function calcExplicitConcentrationsBoundaryConstant(conc_center::Vector{T}, conc
|
||||
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}
|
||||
function writeSolutionVector!(sv::Vector{T}, 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}}, rowIndex::Int, sx::T, sy::T) where {T}
|
||||
numRows = size(concentrations, 1)
|
||||
sv = zeros(T, length)
|
||||
length = size(sv, 1)
|
||||
|
||||
# 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, :]
|
||||
@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
|
||||
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!")
|
||||
@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
|
||||
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!")
|
||||
@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
|
||||
|
||||
# Conditions for the first and last columns
|
||||
is_first_column = (1:length) .== 1
|
||||
is_last_column = (1:length) .== length
|
||||
|
||||
# Apply operations conditionally
|
||||
# First column - additional fixed concentration change from perpendicular dimension in constant BC case
|
||||
if getType(bcLeft[rowIndex]) == CONSTANT
|
||||
sv .+= (2 * sx * alphaX[rowIndex, 1] * getValue(bcLeft[rowIndex])) .* is_first_column
|
||||
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 .+= (2 * sx * alphaX[rowIndex, end] * getValue(bcRight[rowIndex])) .* is_last_column
|
||||
sv[end] += 2 * sx * alphaX[rowIndex, end] * getValue(bcRight[rowIndex])
|
||||
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
|
||||
@ -220,7 +164,7 @@ function BTCS_1D(grid::Grid{T}, bc::Boundary{T}, alpha_left::Matrix{T}, alpha_ri
|
||||
|
||||
concentrations::Matrix{T} = getConcentrations(grid)
|
||||
|
||||
A::Tridiagonal{T} = createCoeffMatrix(alpha, alpha_left, alpha_right, bcLeft, bcRight, length, 1, sx)
|
||||
A::Tridiagonal{T} = createCoeffMatrix(alpha, alpha_left[1, :], alpha_right[1, :], bcLeft, bcRight, length, 1, sx)
|
||||
b = concentrations[1, :]
|
||||
|
||||
if getType(bcLeft[1]) == CONSTANT
|
||||
@ -256,54 +200,28 @@ function BTCS_2D(grid::Grid{T}, bc::Boundary{T}, alphaX_left::Matrix{T}, alphaX_
|
||||
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
|
||||
localBs = [zeros(T, cols) for _ in 1:Threads.nthreads()]
|
||||
Threads.@threads for i = 1:rows
|
||||
A::Tridiagonal{T} = createCoeffMatrix(alphaX, alphaX_left, alphaX_right, bcLeft, bcRight, cols, i, sx)
|
||||
localB = localBs[Threads.threadid()]
|
||||
A::Tridiagonal{T} = createCoeffMatrix(alphaX, alphaX_left[i, :], alphaX_right[i, :], bcLeft, bcRight, cols, i, sx)
|
||||
writeSolutionVector!(localB, concentrations, alphaX, alphaY, bcLeft, bcRight, bcTop, bcBottom, i, sx, sy)
|
||||
|
||||
# Extract the solution vector for the current row from b_matrix
|
||||
b = b_matrix[i, :]
|
||||
|
||||
concentrations_intermediate[i, :] = A \ b
|
||||
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
|
||||
|
||||
# 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
|
||||
localBs = [zeros(T, rows) for _ in 1:Threads.nthreads()]
|
||||
Threads.@threads for i = 1:cols
|
||||
A::Tridiagonal{T} = createCoeffMatrix(alphaY_t, alphaY_t_left, alphaY_t_right, bcTop, bcBottom, rows, i, sy)
|
||||
localB = localBs[Threads.threadid()]
|
||||
A::Tridiagonal{T} = createCoeffMatrix(alphaY_t, alphaY_t_left[i, :], alphaY_t_right[i, :], bcTop, bcBottom, rows, i, sy)
|
||||
writeSolutionVector!(localB, concentrations_intermediate, alphaY_t, alphaX_t, bcTop, bcBottom, bcLeft, bcRight, i, sy, sx)
|
||||
|
||||
# Extract the solution vector for the current row from b_matrix
|
||||
b = b_matrix[i, :]
|
||||
|
||||
concentrations_t[i, :] = A \ b
|
||||
concentrations_t[i, :] = A \ localB
|
||||
end
|
||||
|
||||
concentrations = copy(transpose(concentrations_t))
|
||||
@ -349,7 +267,7 @@ function runBTCS(grid::Grid{T}, bc::Boundary{T}, timestep::T, iterations::Int, s
|
||||
stepCallback()
|
||||
end
|
||||
else
|
||||
error("Error: Only 1- and 2-dimensional grids are defined!")
|
||||
error("BTCS only implemented for 1D and 2D grids")
|
||||
end
|
||||
|
||||
end
|
||||
|
||||
Loading…
x
Reference in New Issue
Block a user