From 957f73bb8356ebccd0d332170c09d2e75aadfe59 Mon Sep 17 00:00:00 2001 From: nebmit <76664673+nebmit@users.noreply.github.com> Date: Mon, 20 Nov 2023 12:16:15 +0100 Subject: [PATCH] refactor!: structural changes Improved julia structs and removed redundant calculations [skip ci] --- ...{BTCS_1_45_750.cpp => BTCS_1_45_75000.cpp} | 2 +- julia/tests/julia_bench/BTCS_100_100_1000.jl | 9 +- julia/tests/julia_bench/BTCS_1_20_100.jl | 7 +- .../{BTCS_1_45_750.jl => BTCS_1_45_75000.jl} | 13 ++- julia/tests/julia_bench/BTCS_20_20_500.jl | 9 +- julia/tests/julia_bench/BTCS_450_670_100.jl | 17 ++-- julia/tests/test.py | 31 +++--- julia/tug/Core/BTCS.jl | 84 ++++++++-------- julia/tug/Grid.jl | 56 ++++++----- julia/tug/Simulation.jl | 95 ++++++++++--------- 10 files changed, 165 insertions(+), 158 deletions(-) rename julia/tests/cpp_bench/{BTCS_1_45_750.cpp => BTCS_1_45_75000.cpp} (96%) rename julia/tests/julia_bench/{BTCS_1_45_750.jl => BTCS_1_45_75000.jl} (86%) diff --git a/julia/tests/cpp_bench/BTCS_1_45_750.cpp b/julia/tests/cpp_bench/BTCS_1_45_75000.cpp similarity index 96% rename from julia/tests/cpp_bench/BTCS_1_45_750.cpp rename to julia/tests/cpp_bench/BTCS_1_45_75000.cpp index b98b1f8..5514c0f 100644 --- a/julia/tests/cpp_bench/BTCS_1_45_750.cpp +++ b/julia/tests/cpp_bench/BTCS_1_45_75000.cpp @@ -27,7 +27,7 @@ int main(int argc, char *argv[]) // **** SIMULATION **** Simulation simulation = Simulation(grid, bc); simulation.setTimestep(1.23); - simulation.setIterations(750); + simulation.setIterations(75000); simulation.setOutputCSV(CSV_OUTPUT_VERBOSE); simulation.setOutputConsole(CONSOLE_OUTPUT_OFF); diff --git a/julia/tests/julia_bench/BTCS_100_100_1000.jl b/julia/tests/julia_bench/BTCS_100_100_1000.jl index ed818e8..671c85f 100644 --- a/julia/tests/julia_bench/BTCS_100_100_1000.jl +++ b/julia/tests/julia_bench/BTCS_100_100_1000.jl @@ -4,17 +4,16 @@ function main() # **** GRID **** rows::Int = 100 cols::Int = 100 - grid::Grid = Grid{Float64}(rows, cols) + + alphaX = fill(1.0, rows, cols) + alphaY = fill(1.0, rows, cols) + grid::Grid = Grid{Float64}(rows, cols, alphaX, alphaY) concentrations = fill(0.0, rows, cols) concentrations[11, 11] = 2000 concentrations[91, 91] = 2000 setConcentrations!(grid, concentrations) - alphaX = fill(1.0, rows, cols) - alphaY = fill(1.0, rows, cols) - setAlpha!(grid, alphaX, alphaY) - # **** BOUNDARY **** bc::Boundary = Boundary(grid) setBoundarySideConstant!(bc, LEFT, 1.0) diff --git a/julia/tests/julia_bench/BTCS_1_20_100.jl b/julia/tests/julia_bench/BTCS_1_20_100.jl index 2267425..d413617 100644 --- a/julia/tests/julia_bench/BTCS_1_20_100.jl +++ b/julia/tests/julia_bench/BTCS_1_20_100.jl @@ -3,15 +3,14 @@ include("../../tug/Simulation.jl") function main() # **** GRID **** cells::Int = 20 - grid::Grid = Grid{Float64}(cells) + + alpha = fill(1.0, 1, cells) + grid::Grid = Grid{Float64}(cells, alpha) concentrations = fill(0.0, 1, cells) concentrations[1] = 2000 setConcentrations!(grid, concentrations) - alpha = fill(1.0, 1, cells) - setAlpha!(grid, alpha) - # **** BOUNDARY **** bc::Boundary = Boundary(grid) setBoundarySideConstant!(bc, LEFT, 0.0) diff --git a/julia/tests/julia_bench/BTCS_1_45_750.jl b/julia/tests/julia_bench/BTCS_1_45_75000.jl similarity index 86% rename from julia/tests/julia_bench/BTCS_1_45_750.jl rename to julia/tests/julia_bench/BTCS_1_45_75000.jl index 5f38221..baeb5b8 100644 --- a/julia/tests/julia_bench/BTCS_1_45_750.jl +++ b/julia/tests/julia_bench/BTCS_1_45_75000.jl @@ -3,16 +3,15 @@ include("../../tug/Simulation.jl") function main() # **** GRID **** cells::Int = 45 - grid::Grid = Grid{Float64}(cells) - - concentrations = fill(10.0, 1, cells) - concentrations[6] = 2000 - setConcentrations!(grid, concentrations) alpha = fill(1.0, 1, cells) alpha[1:15] .= 0.5 alpha[31:45] .= 1.5 - setAlpha!(grid, alpha) + grid::Grid = Grid{Float64}(cells, alpha) + + concentrations = fill(10.0, 1, cells) + concentrations[6] = 2000 + setConcentrations!(grid, concentrations) # **** BOUNDARY **** bc::Boundary = Boundary(grid) @@ -22,7 +21,7 @@ function main() # **** SIMULATION **** simulation::Simulation = Simulation(grid, bc) simulation = setTimestep(simulation, 1.23) - simulation = setIterations(simulation, 750) + simulation = setIterations(simulation, 75000) simulation = setOutputConsole(simulation, CONSOLE_OUTPUT_OFF) simulation = setOutputCSV(simulation, CSV_OUPUT_VERBOSE) diff --git a/julia/tests/julia_bench/BTCS_20_20_500.jl b/julia/tests/julia_bench/BTCS_20_20_500.jl index 5ec0e18..262fb01 100644 --- a/julia/tests/julia_bench/BTCS_20_20_500.jl +++ b/julia/tests/julia_bench/BTCS_20_20_500.jl @@ -4,16 +4,15 @@ function main() # **** GRID **** rows::Int = 20 cols::Int = 20 - grid::Grid = Grid{Float64}(rows, cols) + + alphaX = fill(1.0, rows, cols) + alphaY = fill(1.0, rows, cols) + grid::Grid = Grid{Float64}(rows, cols, alphaX, alphaY) concentrations = fill(0.0, rows, cols) concentrations[11, 11] = 2000 setConcentrations!(grid, concentrations) - alphaX = fill(1.0, rows, cols) - alphaY = fill(1.0, rows, cols) - setAlpha!(grid, alphaX, alphaY) - # **** BOUNDARY **** bc::Boundary = Boundary(grid) setBoundarySideClosed!(bc, LEFT) diff --git a/julia/tests/julia_bench/BTCS_450_670_100.jl b/julia/tests/julia_bench/BTCS_450_670_100.jl index cb2dfba..eec3f37 100644 --- a/julia/tests/julia_bench/BTCS_450_670_100.jl +++ b/julia/tests/julia_bench/BTCS_450_670_100.jl @@ -4,7 +4,14 @@ function main() # **** GRID **** rows::Int = 450 cols::Int = 670 - grid::Grid = Grid{Float64}(rows, cols) + + alphaX = fill(1.0, rows, cols) + alphaY = fill(1.0, rows, cols) + alphaX[1:100, :] .= 0.5 + alphaX[101:200, :] .= 0.8 + alphaY[:, 1:200] .= 0.6 + alphaY[:, 201:400] .= 0.9 + grid::Grid = Grid{Float64}(rows, cols, alphaX, alphaY) concentrations = fill(0.0, rows, cols) concentrations[11, 11] = 1500 @@ -14,14 +21,6 @@ function main() concentrations[221, 336] = 1500 setConcentrations!(grid, concentrations) - alphaX = fill(1.0, rows, cols) - alphaY = fill(1.0, rows, cols) - alphaX[1:100, :] .= 0.5 - alphaX[101:200, :] .= 0.8 - alphaY[:, 1:200] .= 0.6 - alphaY[:, 201:400] .= 0.9 - setAlpha!(grid, alphaX, alphaY) - # **** BOUNDARY **** bc::Boundary = Boundary(grid) setBoundarySideClosed!(bc, LEFT) diff --git a/julia/tests/test.py b/julia/tests/test.py index 0642ae3..9241072 100644 --- a/julia/tests/test.py +++ b/julia/tests/test.py @@ -28,7 +28,14 @@ def get_max_name_length(directory): return max_length def format_difference(diff): - return '{:.5f}'.format(diff).rjust(8) if diff != 0 else '0'.rjust(8) + threshold = 1e-5 + if diff != 0: + if abs(diff) < threshold: + return '{:.2e}'.format(diff).rjust(9) # Scientific notation for small values + else: + return '{:.5f}'.format(diff).rjust(9) # Fixed-point notation for larger values + else: + return '0'.rjust(9) def run_benchmark(command, runs): total_time = 0 @@ -38,11 +45,11 @@ def run_benchmark(command, runs): total_time += time.time() - start_time return total_time / runs -def main(tolerance, runs, silent): +def main(tolerance, runs, silent, no_clean): BENCHMARK_DIR = "./cpp_bench" JULIA_DIR = "./julia_bench" COMPILER = "g++" - CFLAGS = ["-O3", "-fopenmp", "-I", "/usr/local/include", "-I", "../../include/", "-I", "/usr/include/eigen3"] + CFLAGS = ["-O3", "-fopenmp", "-I", "../../include/", "-I", "/usr/include/eigen3"] BIN_DIR = "./cpp_bin_temp" OUTPUT_DIR = "./csv_temp" @@ -104,12 +111,13 @@ def main(tolerance, runs, silent): # Clean up - remove_non_empty_dir(BIN_DIR) - remove_non_empty_dir(OUTPUT_DIR) - - for file in os.listdir('.'): - if file.endswith('.csv'): - os.remove(file) + if not no_clean: + remove_non_empty_dir(BIN_DIR) + remove_non_empty_dir(OUTPUT_DIR) + + for file in os.listdir('.'): + if file.endswith('.csv'): + os.remove(file) # Print results if not silent: print("\n----- Benchmark Results -----") @@ -124,8 +132,9 @@ def main(tolerance, runs, silent): if __name__ == "__main__": parser = argparse.ArgumentParser(description='Benchmark and Compare Script') - parser.add_argument('--tolerance', type=float, default=0.005, help='Tolerance for CSV comparison') + parser.add_argument('--tolerance', type=float, default=0, help='Tolerance for CSV comparison') parser.add_argument('--runs', type=int, default=1, help='Number of benchmark runs') parser.add_argument('--silent', action='store_true', help='Run in silent mode without printing details') + parser.add_argument('--no-clean', action='store_true', help='Do not clean up temporary files') args = parser.parse_args() - main(args.tolerance, args.runs, args.silent) + main(args.tolerance, args.runs, args.silent, args.no_clean) diff --git a/julia/tug/Core/BTCS.jl b/julia/tug/Core/BTCS.jl index 1046b77..b84950b 100644 --- a/julia/tug/Core/BTCS.jl +++ b/julia/tug/Core/BTCS.jl @@ -9,26 +9,21 @@ using SparseArrays include("../Boundary.jl") include("../Grid.jl") -# Helper functions and types -function calcAlphaIntercell(alpha1::T, alpha2::T, useHarmonic::Bool=true) where {T} - if useHarmonic - return 2 / ((1 / alpha1) + (1 / alpha2)) - else - return 0.5 * (alpha1 + alpha2) - end +function calcAlphaIntercell(alpha1::T, alpha2::T) where {T} + return 2 / ((1 / alpha1) + (1 / alpha2)) end -# calculates coefficient for boundary in constant case function calcBoundaryCoeffConstant(alpha_center::T, alpha_side::T, sx::T) where {T} - centerCoeff = 1 + sx * (calcAlphaIntercell(alpha_center, alpha_side) + 2 * alpha_center) - sideCoeff = -sx * calcAlphaIntercell(alpha_center, alpha_side) + alpha = calcAlphaIntercell(alpha_center, alpha_side) + centerCoeff = 1 + sx * (alpha + 2 * alpha_center) + sideCoeff = -sx * alpha return (centerCoeff, sideCoeff) end -# calculates coefficient for boundary in closed case function calcBoundaryCoeffClosed(alpha_center::T, alpha_side::T, sx::T) where {T} - centerCoeff = 1 + sx * calcAlphaIntercell(alpha_center, alpha_side) - sideCoeff = -sx * calcAlphaIntercell(alpha_center, alpha_side) + alpha = calcAlphaIntercell(alpha_center, alpha_side) + centerCoeff = 1 + sx * alpha + sideCoeff = -sx * alpha return (centerCoeff, sideCoeff) end @@ -51,10 +46,13 @@ function createCoeffMatrix(alpha::Matrix{T}, bcLeft::Vector{BoundaryElement{T}}, end # inner columns - for i in 2:(numCols-1) - cm[i, i-1] = -sx * calcAlphaIntercell(alpha[rowIndex, i-1], alpha[rowIndex, i]) - cm[i, i] = 1 + sx * (calcAlphaIntercell(alpha[rowIndex, i], alpha[rowIndex, i+1]) + calcAlphaIntercell(alpha[rowIndex, i-1], alpha[rowIndex, i])) - cm[i, i+1] = -sx * calcAlphaIntercell(alpha[rowIndex, i], alpha[rowIndex, i+1]) + @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 @@ -74,37 +72,38 @@ function createCoeffMatrix(alpha::Matrix{T}, bcLeft::Vector{BoundaryElement{T}}, end - - function calcExplicitConcentrationsBoundaryClosed(conc_center::T, alpha_center::T, alpha_neighbor::T, sy::T) where {T} - sy * calcAlphaIntercell(alpha_center, alpha_neighbor) * conc_center + - (1 - sy * calcAlphaIntercell(alpha_center, alpha_neighbor)) * conc_center + 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} - sy * calcAlphaIntercell(alpha_center, alpha_neighbor) * conc_center + - (1 - sy * (calcAlphaIntercell(alpha_center, alpha_center) + 2 * alpha_center)) * conc_center + + 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 - for i = 1:length - sv[i] = sy * calcAlphaIntercell(alphaY[rowIndex, i], alphaY[rowIndex+1, i]) * concentrations[rowIndex+1, i] + - (1 - sy * (calcAlphaIntercell(alphaY[rowIndex, i], alphaY[rowIndex+1, i]) + calcAlphaIntercell(alphaY[rowIndex-1, i], alphaY[rowIndex, i]))) * concentrations[rowIndex, i] + - sy * calcAlphaIntercell(alphaY[rowIndex-1, i], alphaY[rowIndex, i]) * concentrations[rowIndex-1, i] + @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 - for i = 1:length + @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 @@ -117,7 +116,7 @@ function createSolutionVector(concentrations::Matrix{T}, alphaX::Matrix{T}, alph # Last row if rowIndex == numRows - for i = 1:length + @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 @@ -141,7 +140,6 @@ function createSolutionVector(concentrations::Matrix{T}, alphaX::Matrix{T}, alph 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 @@ -154,7 +152,7 @@ function BTCS_1D(grid::Grid{T}, bc::Boundary{T}, timestep::T, solverFunc::Functi b = Vector{T}(undef, length) - alpha = grid.alphaX[] + alpha = getAlphaX(grid) bcLeft = getBoundarySide(bc, LEFT) bcRight = getBoundarySide(bc, RIGHT) @@ -162,9 +160,10 @@ function BTCS_1D(grid::Grid{T}, bc::Boundary{T}, timestep::T, solverFunc::Functi concentrations = grid.concentrations[] rowIndex = 1 A = createCoeffMatrix(alpha, bcLeft, bcRight, length, rowIndex, sx) - for i in 1:length + @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 @@ -174,7 +173,7 @@ function BTCS_1D(grid::Grid{T}, bc::Boundary{T}, timestep::T, solverFunc::Functi concentrations_t1 = solverFunc(A, b) - for j in 1:length + @inbounds for j in 1:length concentrations[1, j] = concentrations_t1[j] end @@ -188,12 +187,11 @@ function BTCS_2D(grid::Grid{T}, bc::Boundary{T}, timestep::T, solverFunc::Functi sx = timestep / (2 * grid.deltaCol * grid.deltaCol) sy = timestep / (2 * grid.deltaRow * grid.deltaRow) - A = spzeros(T, rowMax, rowMax) concentrations_t1 = zeros(T, rowMax, colMax) row_t1 = Vector{T}(undef, colMax) - alphaX = grid.alphaX[] - alphaY = grid.alphaY[] + alphaX = getAlphaX(grid) + alphaY = getAlphaY(grid) bcLeft = getBoundarySide(bc, LEFT) bcRight = getBoundarySide(bc, RIGHT) @@ -202,7 +200,7 @@ function BTCS_2D(grid::Grid{T}, bc::Boundary{T}, timestep::T, solverFunc::Functi concentrations = grid.concentrations[] - for i = 1:rowMax + @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) @@ -213,10 +211,10 @@ function BTCS_2D(grid::Grid{T}, bc::Boundary{T}, timestep::T, solverFunc::Functi concentrations_t1 = copy(transpose(concentrations_t1)) concentrations = copy(transpose(concentrations)) - alphaX = copy(transpose(alphaX)) - alphaY = copy(transpose(alphaY)) + alphaX = getAlphaX_t(grid) + alphaY = getAlphaY_t(grid) - for i = 1:colMax + @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) @@ -231,9 +229,7 @@ function BTCS_2D(grid::Grid{T}, bc::Boundary{T}, timestep::T, solverFunc::Functi setConcentrations!(grid, concentrations) end - -# Entry point for EigenLU solver; differentiate between 1D and 2D grid -function BTCS_LU(grid::Grid{T}, bc::Boundary{T}, timestep::T, numThreads::Int=1) where {T} +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 diff --git a/julia/tug/Grid.jl b/julia/tug/Grid.jl index 1c3d476..63f4461 100644 --- a/julia/tug/Grid.jl +++ b/julia/tug/Grid.jl @@ -9,25 +9,36 @@ struct Grid{T} deltaCol::T deltaRow::T concentrations::Ref{Matrix{T}} - alphaX::Ref{Matrix{T}} - alphaY::Ref{Matrix{T}} + alphaX::Matrix{T} + alphaY::Union{Matrix{T},Nothing} + alphaX_t::Union{Matrix{T},Nothing} + alphaY_t::Union{Matrix{T},Nothing} # Constructor for 1D-Grid - function Grid{T}(length::Int) where {T} + function Grid{T}(length::Int, alpha::Matrix{T}) where {T} if length <= 3 throw(ArgumentError("Given grid length too small. Must be greater than 3.")) end + if size(alpha, 1) != 1 || size(alpha, 2) != length + error("Given matrix of alpha coefficients mismatch with Grid dimensions!") + end - new{T}(length, 1, 1, T(length), 0, T(1), 0, Ref(fill(T(0), 1, length)), Ref(fill(T(0), 1, length)), Ref(fill(T(0), 1, length))) + new{T}(length, 1, 1, T(length), 0, T(1), 0, Ref(fill(T(0), 1, length)), alpha, nothing, nothing, nothing) end # Constructor for 2D-Grid - function Grid{T}(row::Int, col::Int) where {T} - if row <= 3 || col <= 3 + function Grid{T}(rows::Int, cols::Int, alphaX::Matrix{T}, alphaY::Matrix{T}) where {T} + if rows <= 3 || cols <= 3 throw(ArgumentError("Given grid dimensions too small. Must each be greater than 3.")) end + if size(alphaX) != (rows, cols) || size(alphaY) != (rows, cols) + error("Given matrices of alpha coefficients mismatch with Grid dimensions!") + end - new{T}(col, row, 2, T(col), T(row), T(1), T(1), Ref(fill(T(0), row, col)), Ref(fill(T(0), row, col)), Ref(fill(T(0), row, col))) + alphaX_t = alphaX' + alphaY_t = alphaY' + + new{T}(cols, rows, 2, T(cols), T(rows), T(1), T(1), Ref(fill(T(0), rows, cols)), alphaX, alphaY, alphaX_t, alphaY_t) end end @@ -35,25 +46,18 @@ function setConcentrations!(grid::Grid{T}, new_concentrations::Matrix{T}) where grid.concentrations[] = new_concentrations end -function setAlpha!(grid::Grid{T}, alpha::Matrix{T}) where {T} - if grid.dim != 1 - error("Grid is not one dimensional, you should probably use the 2D setter function!") - end - if size(alpha, 1) != 1 || size(alpha, 2) != grid.cols - error("Given matrix of alpha coefficients mismatch with Grid dimensions!") - end - - grid.alphaX[] = alpha +function getAlphaX(grid::Grid{T})::Matrix{T} where {T} + grid.alphaX end -function setAlpha!(grid::Grid{T}, alphaX::Matrix{T}, alphaY::Matrix{T}) where {T} - if grid.dim != 2 - error("Grid is not two dimensional, you should probably use the 1D setter function!") - end - if size(alphaX) != (grid.rows, grid.cols) || size(alphaY) != (grid.rows, grid.cols) - error("Given matrices of alpha coefficients mismatch with Grid dimensions!") - end - - grid.alphaX[] = alphaX - grid.alphaY[] = alphaY +function getAlphaY(grid::Grid{T})::Matrix{T} where {T} + grid.alphaY +end + +function getAlphaX_t(grid::Grid{T})::Matrix{T} where {T} + grid.alphaX_t +end + +function getAlphaY_t(grid::Grid{T})::Matrix{T} where {T} + grid.alphaY_t end diff --git a/julia/tug/Simulation.jl b/julia/tug/Simulation.jl index 70fdc9a..ce440f8 100644 --- a/julia/tug/Simulation.jl +++ b/julia/tug/Simulation.jl @@ -16,7 +16,6 @@ struct Simulation{T,approach,solver} approach::APPROACH solver::SOLVER - innerIterations::Int iterations::Int numThreads::Int timestep::T @@ -25,16 +24,14 @@ struct Simulation{T,approach,solver} csvOutput::CSV_OUPUT function Simulation(grid::Grid{T}, bc::Boundary{T}, approach::APPROACH=BTCS, - solver::SOLVER=EIGEN_LU_SOLVER, innerIterations::Int=1, iterations::Int=1, - numThreads::Int=1, timestep::T=0.1, + solver::SOLVER=EIGEN_LU_SOLVER, iterations::Int=1, numThreads::Int=1, timestep::T=0.1, consoleOutput::CONSOLE_OUTPUT=CONSOLE_OUTPUT_OFF, csvOutput::CSV_OUPUT=CSV_OUPUT_OFF) where {T} - new{T,APPROACH,SOLVER}(grid, bc, approach, solver, innerIterations, iterations, numThreads, timestep, consoleOutput, csvOutput) + new{T,APPROACH,SOLVER}(grid, bc, approach, solver, iterations, numThreads, timestep, consoleOutput, csvOutput) end end - -function createCSVfile(simulation::Simulation{T,approach,solver}) where {T,approach,solver} +function createCSVfile(simulation::Simulation{T,approach,solver})::IOStream where {T,approach,solver} appendIdent = 0 approachString = (simulation.approach == BTCS) ? "BTCS" : "UNKNOWN" # Add other approaches as needed row = simulation.grid.rows @@ -47,9 +44,9 @@ function createCSVfile(simulation::Simulation{T,approach,solver}) where {T,appro filename = string(approachString, "_", row, "_", col, "_", numIterations, "-", appendIdent, ".csv") end - open(filename, "w") do file - # Write boundary conditions if required - if simulation.csvOutput == CSV_OUPUT_XTREME + # Write boundary conditions if required + if simulation.csvOutput == CSV_OUPUT_XTREME + open(filename, "w") do file writeBoundarySideValues(file, simulation.bc, LEFT) writeBoundarySideValues(file, simulation.bc, RIGHT) @@ -62,7 +59,8 @@ function createCSVfile(simulation::Simulation{T,approach,solver}) where {T,appro end end - return filename + file = open(filename, "a") + return file end function writeBoundarySideValues(file, bc::Boundary{T}, side) where {T} @@ -71,17 +69,15 @@ function writeBoundarySideValues(file, bc::Boundary{T}, side) where {T} write(file, formatted_values, "\n") end - -function printConcentrationsCSV(simulation::Simulation{T,approach,solver}, filename::String) where {T,approach,solver} +function printConcentrationsCSV(simulation::Simulation{T,approach,solver}, file::IOStream) where {T,approach,solver} concentrations = simulation.grid.concentrations[] - open(filename, "a") do file # Open file in append mode - for row in eachrow(concentrations) - println(file, join(row, " ")) - end - println(file) # Add extra newlines for separation - println(file) + for row in eachrow(concentrations) + formatted_row = [Printf.@sprintf("%.6g", x) for x in row] # Format each element like is done in the C++ version using Eigen3 + println(file, join(formatted_row, " ")) end + println(file) # Add extra newlines for separation + println(file) end function printConcentrations(simulation::Simulation{T,approach,solver}) where {T,approach,solver} @@ -89,52 +85,59 @@ function printConcentrations(simulation::Simulation{T,approach,solver}) where {T end function run(simulation::Simulation{T,approach,solver}) where {T,approach,solver} - filename::String = "" - if simulation.csvOutput > CSV_OUPUT_OFF - filename = createCSVfile(simulation) - end + file = nothing + try + if simulation.csvOutput > CSV_OUPUT_OFF + file = createCSVfile(simulation) + end - if simulation.approach == BTCS - if simulation.solver == EIGEN_LU_SOLVER - for i in 1:(simulation.iterations*simulation.innerIterations) - if simulation.consoleOutput == CONSOLE_OUTPUT_VERBOSE - printConcentrations(simulation) + if simulation.approach == BTCS + if simulation.solver == EIGEN_LU_SOLVER + for _ in 1:(simulation.iterations) + if simulation.consoleOutput >= CONSOLE_OUTPUT_VERBOSE + printConcentrations(simulation) + end + + if simulation.csvOutput >= CSV_OUPUT_VERBOSE + printConcentrationsCSV(simulation, file) + end + + BTCS_step(simulation.grid, simulation.bc, simulation.timestep, simulation.numThreads) end - - if simulation.csvOutput >= CSV_OUPUT_VERBOSE - printConcentrationsCSV(simulation, filename) - end - - BTCS_LU(simulation.grid, simulation.bc, simulation.timestep, simulation.numThreads) + else + error("Undefined solver!") end else - error("Undefined solver!") + error("Undefined approach!") end - else - error("Undefined approach!") - end - if simulation.consoleOutput == CONSOLE_OUTPUT_ON || simulation.consoleOutput == CONSOLE_OUTPUT_VERBOSE - printConcentrations(simulation) - end + if simulation.consoleOutput == CONSOLE_OUTPUT_ON || simulation.consoleOutput == CONSOLE_OUTPUT_VERBOSE + printConcentrations(simulation) + end - if simulation.csvOutput == CSV_OUPUT_ON || simulation.csvOutput == CSV_OUPUT_VERBOSE || simulation.csvOutput == CSV_OUPUT_XTREME - printConcentrationsCSV(simulation, filename) + if simulation.csvOutput == CSV_OUPUT_ON || simulation.csvOutput == CSV_OUPUT_VERBOSE || simulation.csvOutput == CSV_OUPUT_XTREME + printConcentrationsCSV(simulation, file) + end + + finally + if file !== nothing + close(file) + end end end function setTimestep(simulation::Simulation{T,approach,solver}, timestep::T) where {T,approach,solver} - return Simulation(simulation.grid, simulation.bc, simulation.approach, simulation.solver, simulation.innerIterations, simulation.iterations, simulation.numThreads, timestep, simulation.consoleOutput, simulation.csvOutput) + return Simulation(simulation.grid, simulation.bc, simulation.approach, simulation.solver, simulation.iterations, simulation.numThreads, timestep, simulation.consoleOutput, simulation.csvOutput) end function setIterations(simulation::Simulation{T,approach,solver}, iterations::Int) where {T,approach,solver} - return Simulation(simulation.grid, simulation.bc, simulation.approach, simulation.solver, simulation.innerIterations, iterations, simulation.numThreads, simulation.timestep, simulation.consoleOutput, simulation.csvOutput) + return Simulation(simulation.grid, simulation.bc, simulation.approach, simulation.solver, iterations, simulation.numThreads, simulation.timestep, simulation.consoleOutput, simulation.csvOutput) end function setOutputConsole(simulation::Simulation{T,approach,solver}, consoleOutput::CONSOLE_OUTPUT) where {T,approach,solver} - return Simulation(simulation.grid, simulation.bc, simulation.approach, simulation.solver, simulation.innerIterations, simulation.iterations, simulation.numThreads, simulation.timestep, consoleOutput, simulation.csvOutput) + return Simulation(simulation.grid, simulation.bc, simulation.approach, simulation.solver, simulation.iterations, simulation.numThreads, simulation.timestep, consoleOutput, simulation.csvOutput) end function setOutputCSV(simulation::Simulation{T,approach,solver}, csvOutput::CSV_OUPUT) where {T,approach,solver} - return Simulation(simulation.grid, simulation.bc, simulation.approach, simulation.solver, simulation.innerIterations, simulation.iterations, simulation.numThreads, simulation.timestep, simulation.consoleOutput, csvOutput) + return Simulation(simulation.grid, simulation.bc, simulation.approach, simulation.solver, simulation.iterations, simulation.numThreads, simulation.timestep, simulation.consoleOutput, csvOutput) end