From f538658e82f6fb74a6530e6e02b6902c28368b3f Mon Sep 17 00:00:00 2001 From: nebmit <76664673+nebmit@users.noreply.github.com> Date: Tue, 28 Nov 2023 12:20:09 +0100 Subject: [PATCH] perf: removed solution matrix computation approach [skip ci] --- julia/tug/Core/BTCS.jl | 172 +++++++++++------------------------------ 1 file changed, 45 insertions(+), 127 deletions(-) diff --git a/julia/tug/Core/BTCS.jl b/julia/tug/Core/BTCS.jl index fd394fe..38e63d6 100644 --- a/julia/tug/Core/BTCS.jl +++ b/julia/tug/Core/BTCS.jl @@ -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