From 51705e3eefa84559efbe5823184bde1bffd6a4f4 Mon Sep 17 00:00:00 2001 From: nebmit <76664673+nebmit@users.noreply.github.com> Date: Thu, 30 Nov 2023 16:09:03 +0100 Subject: [PATCH] perf: optimize FTCS calculation for enhanced performance Refactored calculation functions to combine and simplify logic. Revised functions to accept direct parameters instead of operating on the entire grid. Resulted in significant performance improvements through enhanced compiler optimization. [skip ci] --- julia/tug/Core/FTCS.jl | 285 +++++++++++++++++++---------------------- 1 file changed, 131 insertions(+), 154 deletions(-) diff --git a/julia/tug/Core/FTCS.jl b/julia/tug/Core/FTCS.jl index 87d9fd5..4afd953 100644 --- a/julia/tug/Core/FTCS.jl +++ b/julia/tug/Core/FTCS.jl @@ -7,152 +7,88 @@ include("../Boundary.jl") include("../Grid.jl") include("Utils.jl") -function calcHorizontalChange(grid::Grid{T}, row::Int, col::Int) where {T} - return calcAlphaIntercell(getAlphaX(grid)[row, col+1], - getAlphaX(grid)[row, col]) * - getConcentrations(grid)[row, col+1] - - (calcAlphaIntercell(getAlphaX(grid)[row, col+1], - getAlphaX(grid)[row, col]) + - calcAlphaIntercell(getAlphaX(grid)[row, col-1], - getAlphaX(grid)[row, col])) * - getConcentrations(grid)[row, col] + - calcAlphaIntercell(getAlphaX(grid)[row, col-1], - getAlphaX(grid)[row, col]) * - getConcentrations(grid)[row, col-1] +function calcHorizontalChange(alphaX_next::T, alphaX_prev::T, alphaX_current::T, + concentration_next::T, concentration_prev::T, concentration_current::T) where {T} + + intercellAlpha_next = calcAlphaIntercell(alphaX_next, alphaX_current) + intercellAlpha_prev = calcAlphaIntercell(alphaX_prev, alphaX_current) + + return intercellAlpha_next * concentration_next - + (intercellAlpha_next + intercellAlpha_prev) * concentration_current + + intercellAlpha_prev * concentration_prev end -function calcHorizontalChangeLeftBoundary(grid::Grid{T}, bc::Boundary{T}, row::Int, col::Int) where {T} - if getBoundaryElementType(bc, LEFT, row) == CONSTANT - return calcHorizontalChangeLeftBoundaryConstant(grid, bc, row, col) - elseif getBoundaryElementType(bc, LEFT, row) == CLOSED - return calcHorizontalChangeLeftBoundaryClosed(grid, row, col) +function calcHorizontalChangeLeftBoundary(boundaryType::TYPE, alphaX_next::T, alphaX_current::T, + concentration_next::T, concentration_current::T, boundaryValue::T) where {T} + intercellAlpha = calcAlphaIntercell(alphaX_next, alphaX_current) + + if boundaryType == CONSTANT + return intercellAlpha * concentration_next - + (intercellAlpha + 2 * alphaX_current) * concentration_current + + 2 * alphaX_current * boundaryValue + elseif boundaryType == CLOSED + return intercellAlpha * (concentration_next - concentration_current) else error("Undefined Boundary Condition Type!") end end -function calcHorizontalChangeLeftBoundaryConstant(grid::Grid{T}, bc::Boundary{T}, row::Int, col::Int) where {T} - return calcAlphaIntercell(getAlphaX(grid)[row, col+1], - getAlphaX(grid)[row, col]) * - getConcentrations(grid)[row, col+1] - - (calcAlphaIntercell(getAlphaX(grid)[row, col+1], - getAlphaX(grid)[row, col]) + - 2 * getAlphaX(grid)[row, col]) * - getConcentrations(grid)[row, col] + - 2 * getAlphaX(grid)[row, col] * - getBoundaryElementValue(bc, LEFT, row) -end +function calcHorizontalChangeRightBoundary(boundaryType::TYPE, alphaX_prev::T, alphaX_current::T, + concentration_prev::T, concentration_current::T, boundaryValue::T) where {T} + intercellAlpha = calcAlphaIntercell(alphaX_prev, alphaX_current) -function calcHorizontalChangeLeftBoundaryClosed(grid::Grid{T}, row::Int, col::Int) where {T} - return calcAlphaIntercell(getAlphaX(grid)[row, col+1], - getAlphaX(grid)[row, col]) * - (getConcentrations(grid)[row, col+1] - - getConcentrations(grid)[row, col]) -end - - -function calcHorizontalChangeRightBoundary(grid::Grid{T}, bc::Boundary{T}, row::Int, col::Int) where {T} - if getBoundaryElementType(bc, RIGHT, row) == CONSTANT - return calcHorizontalChangeRightBoundaryConstant(grid, bc, row, col) - elseif getBoundaryElementType(bc, RIGHT, row) == CLOSED - return calcHorizontalChangeRightBoundaryClosed(grid, row, col) + if boundaryType == CONSTANT + return 2 * alphaX_current * boundaryValue - + (intercellAlpha + 2 * alphaX_current) * concentration_current + + intercellAlpha * concentration_prev + elseif boundaryType == CLOSED + return -intercellAlpha * (concentration_current - concentration_prev) else error("Undefined Boundary Condition Type!") end end -function calcHorizontalChangeRightBoundaryConstant(grid::Grid{T}, bc::Boundary{T}, row::Int, col::Int) where {T} - return 2 * getAlphaX(grid)[row, col] * - getBoundaryElementValue(bc, RIGHT, row) - - (calcAlphaIntercell(getAlphaX(grid)[row, col-1], - getAlphaX(grid)[row, col]) + - 2 * getAlphaX(grid)[row, col]) * - getConcentrations(grid)[row, col] + - calcAlphaIntercell(getAlphaX(grid)[row, col-1], - getAlphaX(grid)[row, col]) * - getConcentrations(grid)[row, col-1] +function calcVerticalChange(alphaY_above::T, alphaY_below::T, alphaY_current::T, + concentration_above::T, concentration_below::T, concentration_current::T) where {T} + intercellAlpha_above = calcAlphaIntercell(alphaY_above, alphaY_current) + intercellAlpha_below = calcAlphaIntercell(alphaY_below, alphaY_current) + + return intercellAlpha_above * concentration_above - + (intercellAlpha_above + intercellAlpha_below) * concentration_current + + intercellAlpha_below * concentration_below end -function calcHorizontalChangeRightBoundaryClosed(grid::Grid{T}, row::Int, col::Int) where {T} - return -(calcAlphaIntercell(getAlphaX(grid)[row, col-1], - getAlphaX(grid)[row, col]) * - (getConcentrations(grid)[row, col] - - getConcentrations(grid)[row, col-1])) -end +function calcVerticalChangeTopBoundary(boundaryType::TYPE, alphaY_above::T, alphaY_current::T, + concentration_above::T, concentration_current::T, boundaryValue::T) where {T} + intercellAlpha = calcAlphaIntercell(alphaY_above, alphaY_current) - -function calcVerticalChange(grid::Grid{T}, row::Int, col::Int) where {T} - return calcAlphaIntercell(getAlphaY(grid)[row+1, col], - getAlphaY(grid)[row, col]) * - getConcentrations(grid)[row+1, col] - - (calcAlphaIntercell(getAlphaY(grid)[row+1, col], - getAlphaY(grid)[row, col]) + - calcAlphaIntercell(getAlphaY(grid)[row-1, col], - getAlphaY(grid)[row, col])) * - getConcentrations(grid)[row, col] + - calcAlphaIntercell(getAlphaY(grid)[row-1, col], - getAlphaY(grid)[row, col]) * - getConcentrations(grid)[row-1, col] -end - -function calcVerticalChangeTopBoundary(grid::Grid{T}, bc::Boundary{T}, row::Int, col::Int) where {T} - if getBoundaryElementType(bc, TOP, col) == CONSTANT - return calcVerticalChangeTopBoundaryConstant(grid, bc, row, col) - elseif getBoundaryElementType(bc, TOP, col) == CLOSED - return calcVerticalChangeTopBoundaryClosed(grid, row, col) + if boundaryType == CONSTANT + return intercellAlpha * concentration_above - + (intercellAlpha + 2 * alphaY_current) * concentration_current + + 2 * alphaY_current * boundaryValue + elseif boundaryType == CLOSED + return intercellAlpha * (concentration_above - concentration_current) else error("Undefined Boundary Condition Type!") end end -function calcVerticalChangeTopBoundaryConstant(grid::Grid{T}, bc::Boundary{T}, row::Int, col::Int) where {T} - return calcAlphaIntercell(getAlphaY(grid)[row+1, col], - getAlphaY(grid)[row, col]) * - getConcentrations(grid)[row+1, col] - - (calcAlphaIntercell(getAlphaY(grid)[row+1, col], - getAlphaY(grid)[row, col]) + - 2 * getAlphaY(grid)[row, col]) * - getConcentrations(grid)[row, col] + - 2 * getAlphaY(grid)[row, col] * - getBoundaryElementValue(bc, TOP, col) -end -function calcVerticalChangeTopBoundaryClosed(grid::Grid{T}, row::Int, col::Int) where {T} - return calcAlphaIntercell(getAlphaY(grid)[row+1, col], - getAlphaY(grid)[row, col]) * - (getConcentrations(grid)[row+1, col] - - getConcentrations(grid)[row, col]) -end +function calcVerticalChangeBottomBoundary(boundaryType::TYPE, alphaY_current::T, alphaY_below::T, + concentration_current::T, concentration_below::T, boundaryValue::T) where {T} + intercellAlpha = calcAlphaIntercell(alphaY_current, alphaY_below) -function calcVerticalChangeBottomBoundary(grid::Grid{T}, bc::Boundary{T}, row::Int, col::Int) where {T} - if getBoundaryElementType(bc, BOTTOM, col) == CONSTANT - return calcVerticalChangeBottomBoundaryConstant(grid, bc, row, col) - elseif getBoundaryElementType(bc, BOTTOM, col) == CLOSED - return calcVerticalChangeBottomBoundaryClosed(grid, row, col) + if boundaryType == CONSTANT + return 2 * alphaY_current * boundaryValue - + (intercellAlpha + 2 * alphaY_current) * concentration_current + + intercellAlpha * concentration_below + elseif boundaryType == CLOSED + return -intercellAlpha * (concentration_current - concentration_below) else error("Undefined Boundary Condition Type!") end end -function calcVerticalChangeBottomBoundaryConstant(grid::Grid{T}, bc::Boundary{T}, row::Int, col::Int) where {T} - return 2 * getAlphaY(grid)[row, col] * - getBoundaryElementValue(bc, BOTTOM, col) - - (calcAlphaIntercell(getAlphaY(grid)[row, col], - getAlphaY(grid)[row-1, col]) + - 2 * getAlphaY(grid)[row, col]) * - getConcentrations(grid)[row, col] + - calcAlphaIntercell(getAlphaY(grid)[row, col], - getAlphaY(grid)[row-1, col]) * - getConcentrations(grid)[row-1, col] -end - -function calcVerticalChangeBottomBoundaryClosed(grid::Grid{T}, row::Int, col::Int) where {T} - return -(calcAlphaIntercell(getAlphaY(grid)[row, col], - getAlphaY(grid)[row-1, col]) * - (getConcentrations(grid)[row, col] - - getConcentrations(grid)[row-1, col])) -end - function FTCS_1D(grid::Grid{T}, bc::Boundary{T}, timestep::T) where {T} @@ -160,78 +96,119 @@ function FTCS_1D(grid::Grid{T}, bc::Boundary{T}, timestep::T) where {T} sx = timestep / (getDeltaCol(grid)^2) concentrations = getConcentrations(grid) + alphaX = getAlphaX(grid) concentrations_t1 = copy(concentrations) - # only one row in 1D case -> row constant at index 1 row = 1 # inner cells - # independent of boundary condition type for col = 2:colMax-1 - concentrations_t1[row, col] += sx * calcHorizontalChange(grid, row, col) + concentrations_t1[row, col] += sx * calcHorizontalChange(alphaX[row, col+1], alphaX[row, col-1], alphaX[row, col], + concentrations[row, col+1], concentrations[row, col-1], concentrations[row, col]) end - # left boundary; hold column constant at index 1 - concentrations_t1[row, 1] += sx * calcHorizontalChangeLeftBoundary(grid, bc, row, 1) + # left boundary + leftBoundaryType = getBoundaryElementType(bc, LEFT, row) + leftBoundaryValue = getBoundaryElementValue(bc, LEFT, row) + concentrations_t1[row, 1] += sx * calcHorizontalChangeLeftBoundary(leftBoundaryType, alphaX[row, 2], alphaX[row, 1], + concentrations[row, 2], concentrations[row, 1], leftBoundaryValue) - # right boundary; hold column constant at max index - concentrations_t1[row, colMax] += sx * calcHorizontalChangeRightBoundary(grid, bc, row, colMax) + # right boundary + rightBoundaryType = getBoundaryElementType(bc, RIGHT, row) + rightBoundaryValue = getBoundaryElementValue(bc, RIGHT, row) + concentrations_t1[row, colMax] += sx * calcHorizontalChangeRightBoundary(rightBoundaryType, alphaX[row, colMax-1], alphaX[row, colMax], + concentrations[row, colMax-1], concentrations[row, colMax], rightBoundaryValue) - # overwrite obsolete concentrations setConcentrations!(grid, concentrations_t1) end + function FTCS_2D(grid::Grid{T}, bc::Boundary{T}, timestep::T) where {T} rowMax = getRows(grid) colMax = getCols(grid) sx = timestep / (getDeltaCol(grid)^2) sy = timestep / (getDeltaRow(grid)^2) + alphaX = getAlphaX(grid) + alphaY = getAlphaY(grid) concentrations = getConcentrations(grid) concentrations_t1 = copy(concentrations) + # Currently boundary sides can only be set for the whole side, not for each cell individually + leftBoundaryType = getBoundaryElementType(bc, LEFT, 1) + rightBoundaryType = getBoundaryElementType(bc, RIGHT, 1) + topBoundaryType = getBoundaryElementType(bc, TOP, 1) + bottomBoundaryType = getBoundaryElementType(bc, BOTTOM, 1) + Threads.@threads for row = 2:rowMax-1 - # inner cells for col = 2:colMax-1 - concentrations_t1[row, col] += sy * calcVerticalChange(grid, row, col) + - sx * calcHorizontalChange(grid, row, col) + concentrations_t1[row, col] += sy * calcVerticalChange(alphaY[row+1, col], alphaY[row-1, col], alphaY[row, col], + concentrations[row+1, col], concentrations[row-1, col], concentrations[row, col]) + + sx * calcHorizontalChange(alphaX[row, col+1], alphaX[row, col-1], alphaX[row, col], + concentrations[row, col+1], concentrations[row, col-1], concentrations[row, col]) end - # left boundary without corners - concentrations_t1[row, 1] += sx * calcHorizontalChangeLeftBoundary(grid, bc, row, 1) + - sy * calcVerticalChange(grid, row, 1) + # Boundary conditions for each row + # Left boundary without corners + leftBoundaryValue = getBoundaryElementValue(bc, LEFT, row) + concentrations_t1[row, 1] += sx * calcHorizontalChangeLeftBoundary(leftBoundaryType, alphaX[row, 2], alphaX[row, 1], + concentrations[row, 2], concentrations[row, 1], leftBoundaryValue) + + sy * calcVerticalChange(alphaY[row+1, 1], alphaY[row-1, 1], alphaY[row, 1], + concentrations[row+1, 1], concentrations[row-1, 1], concentrations[row, 1]) - # right boundary without corners - concentrations_t1[row, colMax] += sx * calcHorizontalChangeRightBoundary(grid, bc, row, colMax) + - sy * calcVerticalChange(grid, row, colMax) + # Right boundary without corners + rightBoundaryValue = getBoundaryElementValue(bc, RIGHT, row) + concentrations_t1[row, colMax] += sx * calcHorizontalChangeRightBoundary(rightBoundaryType, alphaX[row, colMax-1], alphaX[row, colMax], + concentrations[row, colMax-1], concentrations[row, colMax], rightBoundaryValue) + + sy * calcVerticalChange(alphaY[row+1, colMax], alphaY[row-1, colMax], alphaY[row, colMax], + concentrations[row+1, colMax], concentrations[row-1, colMax], concentrations[row, colMax]) end - # top / bottom without corners / looping over columns + # Handle top/bottom boundaries Threads.@threads for col = 2:colMax-1 - # top - concentrations_t1[1, col] += sy * calcVerticalChangeTopBoundary(grid, bc, 1, col) + - sx * calcHorizontalChange(grid, 1, col) + # Top boundary + topBoundaryValue = getBoundaryElementValue(bc, TOP, col) + concentrations_t1[1, col] += sy * calcVerticalChangeTopBoundary(topBoundaryType, alphaY[2, col], alphaY[1, col], + concentrations[2, col], concentrations[1, col], topBoundaryValue) + + sx * calcHorizontalChange(alphaX[1, col+1], alphaX[1, col-1], alphaX[1, col], + concentrations[1, col+1], concentrations[1, col-1], concentrations[1, col]) - # bottom - concentrations_t1[rowMax, col] += sy * calcVerticalChangeBottomBoundary(grid, bc, rowMax, col) + - sx * calcHorizontalChange(grid, rowMax, col) + # Bottom boundary + bottomBoundaryValue = getBoundaryElementValue(bc, BOTTOM, col) + concentrations_t1[rowMax, col] += sy * calcVerticalChangeBottomBoundary(bottomBoundaryType, alphaY[rowMax, col], alphaY[rowMax-1, col], + concentrations[rowMax, col], concentrations[rowMax-1, col], bottomBoundaryValue) + + sx * calcHorizontalChange(alphaX[rowMax, col+1], alphaX[rowMax, col-1], alphaX[rowMax, col], + concentrations[rowMax, col+1], concentrations[rowMax, col-1], concentrations[rowMax, col]) end - # corner top left - concentrations_t1[1, 1] += sx * calcHorizontalChangeLeftBoundary(grid, bc, 1, 1) + - sy * calcVerticalChangeTopBoundary(grid, bc, 1, 1) + # Handle corners + # Top left corner + topBoundaryValue = getBoundaryElementValue(bc, TOP, 1) + concentrations_t1[1, 1] += sx * calcHorizontalChange(alphaX[1, 2], alphaX[1, 1], alphaX[1, 1], + concentrations[1, 2], concentrations[1, 1], concentrations[1, 1]) + + sy * calcVerticalChangeTopBoundary(topBoundaryType, alphaY[2, 1], alphaY[1, 1], + concentrations[2, 1], concentrations[1, 1], topBoundaryValue) - # corner top right - concentrations_t1[1, colMax] += sx * calcHorizontalChangeRightBoundary(grid, bc, 1, colMax) + - sy * calcVerticalChangeTopBoundary(grid, bc, 1, colMax) + # Top right corner + topBoundaryValue = getBoundaryElementValue(bc, TOP, colMax) + concentrations_t1[1, colMax] += sx * calcHorizontalChange(alphaX[1, colMax-1], alphaX[1, colMax], alphaX[1, colMax], + concentrations[1, colMax-1], concentrations[1, colMax], concentrations[1, colMax]) + + sy * calcVerticalChangeTopBoundary(topBoundaryType, alphaY[2, colMax], alphaY[1, colMax], + concentrations[2, colMax], concentrations[1, colMax], topBoundaryValue) - # corner bottom left - concentrations_t1[rowMax, 1] += sx * calcHorizontalChangeLeftBoundary(grid, bc, rowMax, 1) + - sy * calcVerticalChangeBottomBoundary(grid, bc, rowMax, 1) + # Bottom left corner + bottomBoundaryValue = getBoundaryElementValue(bc, BOTTOM, 1) + concentrations_t1[rowMax, 1] += sx * calcHorizontalChange(alphaX[rowMax, 2], alphaX[rowMax, 1], alphaX[rowMax, 1], + concentrations[rowMax, 2], concentrations[rowMax, 1], concentrations[rowMax, 1]) + + sy * calcVerticalChangeBottomBoundary(bottomBoundaryType, alphaY[rowMax, 1], alphaY[rowMax-1, 1], + concentrations[rowMax, 1], concentrations[rowMax-1, 1], bottomBoundaryValue) - # corner bottom right - concentrations_t1[rowMax, colMax] += sx * calcHorizontalChangeRightBoundary(grid, bc, rowMax, colMax) + - sy * calcVerticalChangeBottomBoundary(grid, bc, rowMax, colMax) + # Bottom right corner + bottomBoundaryValue = getBoundaryElementValue(bc, BOTTOM, colMax) + concentrations_t1[rowMax, colMax] += sx * calcHorizontalChange(alphaX[rowMax, colMax-1], alphaX[rowMax, colMax], alphaX[rowMax, colMax], + concentrations[rowMax, colMax-1], concentrations[rowMax, colMax], concentrations[rowMax, colMax]) + + sy * calcVerticalChangeBottomBoundary(bottomBoundaryType, alphaY[rowMax, colMax], alphaY[rowMax-1, colMax], + concentrations[rowMax, colMax], concentrations[rowMax-1, colMax], bottomBoundaryValue) setConcentrations!(grid, concentrations_t1) end