From 0eb96ed0ada3fbc1387755831f4a42ef29ab62bd Mon Sep 17 00:00:00 2001 From: nebmit <76664673+nebmit@users.noreply.github.com> Date: Wed, 22 Nov 2023 17:43:52 +0100 Subject: [PATCH] perf: solution vector creation using matrices 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] --- julia/tug/Boundary.jl | 4 + julia/tug/Core/BTCS.jl | 197 +++++++++++++++++++++++++++++++---------- 2 files changed, 152 insertions(+), 49 deletions(-) diff --git a/julia/tug/Boundary.jl b/julia/tug/Boundary.jl index cb6cb60..d4e87ca 100644 --- a/julia/tug/Boundary.jl +++ b/julia/tug/Boundary.jl @@ -27,6 +27,10 @@ function getValue(be::BoundaryElement{T})::T where {T} be.value end +function getValue(be::Vector{BoundaryElement{T}})::Vector{T} where {T} + [b.value for b in be] +end + function setType!(be::BoundaryElement{T}, type::Symbol) where {T} be.type = type end diff --git a/julia/tug/Core/BTCS.jl b/julia/tug/Core/BTCS.jl index cb51fd0..fd394fe 100644 --- a/julia/tug/Core/BTCS.jl +++ b/julia/tug/Core/BTCS.jl @@ -15,6 +15,10 @@ 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 @@ -35,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::Vector{T}, alpha_right::Vector{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::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]) @@ -60,9 +64,9 @@ function createCoeffMatrix(alpha::Matrix{T}, alpha_left::Vector{T}, alpha_right: error("Undefined Boundary Condition Type on Right!") end - dl = [-sx .* alpha_left; leftCoeffBottom] - du = [rightCoeffTop; -sx .* alpha_right] - d = [centerCoeffTop; 1 .+ sx .* (alpha_right + alpha_left); centerCoeffBottom] + 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 @@ -84,57 +88,124 @@ function calcExplicitConcentrationsBoundaryConstant(conc_center::T, conc_bc::T, sy * alpha_center * conc_bc end -# creates a solution vector for next time step from the current state of concentrations inplace -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} +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) - length = size(sv, 1) + sv = zeros(T, 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 + 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 - @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 + 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 - @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 + 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 - # First column - additional fixed concentration change from perpendicular dimension in constant BC case + # 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[1] += 2 * sx * alphaX[rowIndex, 1] * getValue(bcLeft[rowIndex]) + sv .+= (2 * sx * alphaX[rowIndex, 1] * getValue(bcLeft[rowIndex])) .* is_first_column 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]) + 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 @@ -149,8 +220,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[1, :], alpha_right[1, :], bcLeft, bcRight, length, 1, sx) - + A::Tridiagonal{T} = createCoeffMatrix(alpha, alpha_left, alpha_right, bcLeft, bcRight, length, 1, sx) b = concentrations[1, :] if getType(bcLeft[1]) == CONSTANT @@ -185,26 +255,55 @@ function BTCS_2D(grid::Grid{T}, bc::Boundary{T}, alphaX_left::Matrix{T}, alphaX_ bcTop = getBoundarySide(bc, TOP) bcBottom = getBoundarySide(bc, BOTTOM) - localBs = [zeros(T, cols) for _ in 1:Threads.nthreads()] - Threads.@threads for i = 1:rows - 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) - concentrations_intermediate[i, :] = A \ localB + # 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) - localBs = [zeros(T, rows) for _ in 1:Threads.nthreads()] - Threads.@threads for i = 1:cols - localB = localBs[Threads.threadid()] - # Swap alphas, boundary conditions and sx/sy for column-wise calculation - 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) - concentrations_t[i, :] = (A \ localB) + # 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))