refactor!: organized and added getters/setters

!Removed solver parameter from simulation.

[skip ci]
This commit is contained in:
nebmit 2023-11-22 13:48:09 +01:00
parent cfeb935c93
commit 1cdeb8d7a7
No known key found for this signature in database
4 changed files with 106 additions and 77 deletions

View File

@ -1,5 +1,7 @@
# Boundary.jl
# Julia implementation of Boundary types, translating from C++'s Boundary.hpp
# API of Boundary class, that holds all information for each boundary
# condition at the edges of the diffusion grid.
# Translated from C++'s Boundary.hpp.
@enum TYPE CLOSED CONSTANT
@enum SIDE LEFT = 1 RIGHT = 2 TOP = 3 BOTTOM = 4
@ -17,12 +19,6 @@ struct BoundaryElement{T}
BoundaryElement{T}(value::T) where {T} = new{T}(CONSTANT, value)
end
# Helper functions for setting type and value
function setType!(be::BoundaryElement{T}, type::Symbol) where {T}
be.type = type
end
function getType(be::BoundaryElement{T})::TYPE where {T}
be.type
end
@ -31,6 +27,10 @@ function getValue(be::BoundaryElement{T})::T where {T}
be.value
end
function setType!(be::BoundaryElement{T}, type::Symbol) where {T}
be.type = type
end
function setValue!(be::BoundaryElement{T}, value::T) where {T}
if value < 0
throw(ArgumentError("No negative concentration allowed."))
@ -38,9 +38,11 @@ function setValue!(be::BoundaryElement{T}, value::T) where {T}
if be.type == BC_TYPE_CLOSED
throw(ArgumentError("No constant boundary concentrations can be set for closed boundaries. Please change type first."))
end
be.value = value
end
# Boundary class
struct Boundary{T}
dim::UInt8
@ -71,6 +73,14 @@ struct Boundary{T}
end
end
function getBoundarySide(boundary::Boundary{T}, side::SIDE)::Vector{BoundaryElement{T}} where {T}
if boundary.dim == 1 && (side == BOTTOM || side == TOP)
throw(ArgumentError("For the one-dimensional case, only the left and right borders exist."))
end
boundary.boundaries[Int(side)]
end
function setBoundarySideClosed!(boundary::Boundary{T}, side::SIDE) where {T}
if boundary.dim == 1 && (side == BOTTOM || side == TOP)
throw(ArgumentError("For the one-dimensional case, only the left and right borders exist."))
@ -92,11 +102,3 @@ function setBoundarySideConstant!(boundary::Boundary{T}, side::SIDE, value::T) w
boundary.boundaries[Int(side)] = [BoundaryElement{T}(value) for _ in 1:n]
end
function getBoundarySide(boundary::Boundary{T}, side::SIDE)::Vector{BoundaryElement{T}} where {T}
if boundary.dim == 1 && (side == BOTTOM || side == TOP)
throw(ArgumentError("For the one-dimensional case, only the left and right borders exist."))
end
boundary.boundaries[Int(side)]
end

View File

@ -2,21 +2,29 @@
# Implementation of heterogenous BTCS (backward time-centered space)
# solution of diffusion equation in 1D and 2D space using the
# alternating-direction implicit (ADI) method.
# Translated from C++'s BTCS.hpp.
using LinearAlgebra
using SparseArrays
using Base.Threads
using CUDA
include("../Boundary.jl")
include("../Grid.jl")
function calcAlphaIntercell(alpha1::T, alpha2::T) where {T}
return 2 / ((1 / alpha1) + (1 / alpha2))
2 / ((1 / alpha1) + (1 / alpha2))
end
function calcAlphaIntercell(alpha1::Matrix{T}, alpha2::Matrix{T}) where {T}
return 2 ./ ((1 ./ alpha1) .+ (1 ./ alpha2))
2 ./ ((1 ./ alpha1) .+ (1 ./ alpha2))
end
function calcBoundaryCoeffClosed(alpha_center::T, alpha_side::T, sx::T) where {T}
alpha = calcAlphaIntercell(alpha_center, alpha_side)
centerCoeff = 1 + sx * alpha
sideCoeff = -sx * alpha
return (centerCoeff, sideCoeff)
end
function calcBoundaryCoeffConstant(alpha_center::T, alpha_side::T, sx::T) where {T}
@ -26,13 +34,6 @@ function calcBoundaryCoeffConstant(alpha_center::T, alpha_side::T, sx::T) where
return (centerCoeff, sideCoeff)
end
function calcBoundaryCoeffClosed(alpha_center::T, alpha_side::T, sx::T) where {T}
alpha = calcAlphaIntercell(alpha_center, alpha_side)
centerCoeff = 1 + sx * alpha
sideCoeff = -sx * alpha
return (centerCoeff, sideCoeff)
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}
# Precompute boundary condition type check for efficiency
@ -70,18 +71,20 @@ end
function calcExplicitConcentrationsBoundaryClosed(conc_center::T, alpha_center::T, alpha_neighbor::T, sy::T) where {T}
alpha = calcAlphaIntercell(alpha_center, alpha_neighbor)
sy * alpha * conc_center + (1 - sy * alpha) * conc_center
end
(sy * alpha + (1 - sy * alpha)) * conc_center
end
function calcExplicitConcentrationsBoundaryConstant(conc_center::T, conc_bc::T, alpha_center::T, alpha_neighbor::T, sy::T) where {T}
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
# 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}
numRows = size(concentrations, 1)
length = size(sv, 1)
@ -134,9 +137,10 @@ function writeSolutionVector!(sv::Vector{T}, concentrations::Matrix{T}, alphaX::
end
end
# BTCS solution for 1D grid
function BTCS_1D(grid::Grid{T}, bc::Boundary{T}, alpha_left::Matrix{T}, alpha_right::Matrix{T}, timestep::T) where {T}
sx = timestep / (grid.deltaCol * grid.deltaCol)
sx = timestep / (getDeltaCol(grid) * getDeltaCol(grid))
alpha = getAlphaX(grid)
bcLeft = getBoundarySide(bc, LEFT)
@ -163,8 +167,8 @@ end
function BTCS_2D(grid::Grid{T}, bc::Boundary{T}, alphaX_left::Matrix{T}, alphaX_right::Matrix{T}, alphaY_t_left::Matrix{T}, alphaY_t_right::Matrix{T}, timestep::T) where {T}
rows = getRows(grid)
cols = getCols(grid)
sx = timestep / (2 * grid.deltaCol * grid.deltaCol)
sy = timestep / (2 * grid.deltaRow * grid.deltaRow)
sx = timestep / (2 * getDeltaCol(grid) * getDeltaCol(grid))
sy = timestep / (2 * getDeltaRow(grid) * getDeltaRow(grid))
alphaX = getAlphaX(grid)
alphaY = getAlphaY(grid)
@ -209,7 +213,7 @@ function BTCS_2D(grid::Grid{T}, bc::Boundary{T}, alphaX_left::Matrix{T}, alphaX_
end
function runBTCS(grid::Grid{T}, bc::Boundary{T}, timestep::T, iterations::Int, stepCallback::Function) where {T}
if grid.dim == 1
if getDim(grid) == 1
length = getCols(grid)
alpha = getAlphaX(grid)
@ -224,7 +228,7 @@ function runBTCS(grid::Grid{T}, bc::Boundary{T}, timestep::T, iterations::Int, s
stepCallback()
end
elseif grid.dim == 2
elseif getDim(grid) == 2
rows = getRows(grid)
cols = getCols(grid)
alphaX = getAlphaX(grid)

View File

@ -1,5 +1,11 @@
# Grid.jl
# API of Grid class, that holds a matrix with concenctrations and
# a respective matrix/matrices of alpha coefficients.
# Translated from C++'s Grid.hpp.
using LinearAlgebra
# Grid class
struct Grid{T}
cols::Int
rows::Int
@ -35,6 +41,7 @@ struct Grid{T}
error("Given matrices of alpha coefficients mismatch with Grid dimensions!")
end
# Precompute alphaX_t and alphaY_t
alphaX_t = alphaX'
alphaY_t = alphaY'
@ -58,18 +65,30 @@ function getAlphaY_t(grid::Grid{T})::Matrix{T} where {T}
grid.alphaY_t
end
function getCols(grid::Grid{T})::Int where {T}
grid.cols
end
function getConcentrations(grid::Grid{T})::Matrix{T} where {T}
grid.concentrations[]
end
function setConcentrations!(grid::Grid{T}, new_concentrations::Matrix{T}) where {T}
grid.concentrations[] = new_concentrations
function getDeltaCol(grid::Grid{T})::T where {T}
grid.deltaCol
end
function getCols(grid::Grid{T})::Int where {T}
grid.cols
function getDeltaRow(grid::Grid{T})::T where {T}
grid.deltaRow
end
function getDim(grid::Grid{T})::Int where {T}
grid.dim
end
function getRows(grid::Grid{T})::Int where {T}
grid.rows
end
function setConcentrations!(grid::Grid{T}, new_concentrations::Matrix{T}) where {T}
grid.concentrations[] = new_concentrations
end

View File

@ -1,3 +1,9 @@
# Simulation.jl
# API of Simulation class, that holds all information regarding a
# specific simulation run like its timestep, number of iterations and output
# options. Simulation object also holds a predefined Grid and Boundary object.
# Translated from C++'s Simulation.hpp.
using Printf
include("Grid.jl")
@ -5,16 +11,14 @@ include("Boundary.jl")
include("Core/BTCS.jl")
@enum APPROACH BTCS
@enum SOLVER EIGEN_LU_SOLVER
@enum CONSOLE_OUTPUT CONSOLE_OUTPUT_OFF CONSOLE_OUTPUT_ON CONSOLE_OUTPUT_VERBOSE
@enum CSV_OUTPUT CSV_OUTPUT_OFF CSV_OUTPUT_ON CSV_OUTPUT_VERBOSE CSV_OUTPUT_XTREME
# Create the Simulation class
struct Simulation{T,approach,solver}
# Simulation class
struct Simulation{T}
grid::Grid{T}
bc::Boundary{T}
approach::APPROACH
solver::SOLVER
iterations::Int
timestep::T
@ -22,36 +26,36 @@ struct Simulation{T,approach,solver}
consoleOutput::CONSOLE_OUTPUT
csvOutput::CSV_OUTPUT
function Simulation(grid::Grid{T}, bc::Boundary{T}, approach::APPROACH=BTCS,
solver::SOLVER=EIGEN_LU_SOLVER, iterations::Int=1, timestep::T=0.1,
# Constructor
function Simulation(grid::Grid{T}, bc::Boundary{T}, approach::APPROACH=BTCS, iterations::Int=1, timestep::T=0.1,
consoleOutput::CONSOLE_OUTPUT=CONSOLE_OUTPUT_OFF, csvOutput::CSV_OUTPUT=CSV_OUTPUT_OFF) where {T}
new{T,APPROACH,SOLVER}(grid, bc, approach, solver, iterations, timestep, consoleOutput, csvOutput)
new{T}(grid, bc, approach, iterations, timestep, consoleOutput, csvOutput)
end
end
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
col = simulation.grid.cols
function _createCSVfile(simulation::Simulation{T})::IOStream where {T}
approachString = string(simulation.approach)
rows = getRows(simulation.grid)
cols = getCols(simulation.grid)
numIterations = simulation.iterations
filename = string(approachString, "_", row, "_", col, "_", numIterations, ".csv")
filename = string(approachString, "_", rows, "_", cols, "_", numIterations, ".csv")
appendIdent = 0
while isfile(filename)
appendIdent += 1
filename = string(approachString, "_", row, "_", col, "_", numIterations, "-", appendIdent, ".csv")
filename = string(approachString, "_", rows, "_", cols, "_", numIterations, "-", appendIdent, ".csv")
end
# Write boundary conditions if required
if simulation.csvOutput == CSV_OUTPUT_XTREME
if simulation.csvOutput >= CSV_OUTPUT_XTREME
open(filename, "w") do file
writeBoundarySideValues(file, simulation.bc, LEFT)
writeBoundarySideValues(file, simulation.bc, RIGHT)
_writeBoundarySideValues(file, simulation.bc, LEFT)
_writeBoundarySideValues(file, simulation.bc, RIGHT)
if simulation.grid.dim == 2
writeBoundarySideValues(file, simulation.bc, TOP)
writeBoundarySideValues(file, simulation.bc, BOTTOM)
if getDim(simulation.grid) == 2
_writeBoundarySideValues(file, simulation.bc, TOP)
_writeBoundarySideValues(file, simulation.bc, BOTTOM)
end
write(file, "\n\n")
@ -62,40 +66,40 @@ function createCSVfile(simulation::Simulation{T,approach,solver})::IOStream wher
return file
end
function writeBoundarySideValues(file, bc::Boundary{T}, side) where {T}
function _writeBoundarySideValues(file::IOStream, bc::Boundary{T}, side) where {T}
values::Vector{BoundaryElement} = getBoundarySide(bc, side)
formatted_values = join(map(getValue, values), " ")
write(file, formatted_values, "\n")
end
function printConcentrationsCSV(simulation::Simulation{T,approach,solver}, file::IOStream) where {T,approach,solver}
concentrations = simulation.grid.concentrations[]
function _printConcentrationsCSV(file::IOStream, simulation::Simulation{T}) where {T}
concentrations = getConcentrations(simulation.grid)
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)
println(file)
end
function printConcentrations(simulation::Simulation{T,approach,solver}) where {T,approach,solver}
println(simulation.grid.concentrations[])
function _printConcentrations(simulation::Simulation{T}) where {T}
println(getConcentrations(simulation.grid))
end
function run(simulation::Simulation{T,approach,solver}) where {T,approach,solver}
function run(simulation::Simulation{T}) where {T}
file = nothing
try
if simulation.csvOutput > CSV_OUTPUT_OFF
file = createCSVfile(simulation)
file = _createCSVfile(simulation)
end
function simulationStepCallback()
if simulation.consoleOutput >= CONSOLE_OUTPUT_VERBOSE
printConcentrations(simulation)
_printConcentrations(simulation)
end
if simulation.csvOutput >= CSV_OUTPUT_VERBOSE
printConcentrationsCSV(simulation, file)
_printConcentrationsCSV(file, simulation)
end
end
@ -106,11 +110,11 @@ function run(simulation::Simulation{T,approach,solver}) where {T,approach,solver
end
if simulation.consoleOutput >= CONSOLE_OUTPUT_ON
printConcentrations(simulation)
_printConcentrations(simulation)
end
if simulation.csvOutput >= CSV_OUTPUT_ON
printConcentrationsCSV(simulation, file)
_printConcentrationsCSV(file, simulation)
end
finally
@ -120,18 +124,18 @@ function run(simulation::Simulation{T,approach,solver}) where {T,approach,solver
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.iterations, timestep, simulation.consoleOutput, simulation.csvOutput)
function setIterations(simulation::Simulation{T}, iterations::Int)::Simulation{T} where {T}
return Simulation(simulation.grid, simulation.bc, simulation.approach, iterations, simulation.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, iterations, simulation.timestep, simulation.consoleOutput, simulation.csvOutput)
function setOutputConsole(simulation::Simulation{T}, consoleOutput::CONSOLE_OUTPUT)::Simulation{T} where {T}
return Simulation(simulation.grid, simulation.bc, simulation.approach, simulation.iterations, simulation.timestep, 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.iterations, simulation.timestep, consoleOutput, simulation.csvOutput)
function setOutputCSV(simulation::Simulation{T}, csvOutput::CSV_OUTPUT)::Simulation{T} where {T}
return Simulation(simulation.grid, simulation.bc, simulation.approach, simulation.iterations, simulation.timestep, simulation.consoleOutput, csvOutput)
end
function setOutputCSV(simulation::Simulation{T,approach,solver}, csvOutput::CSV_OUTPUT) where {T,approach,solver}
return Simulation(simulation.grid, simulation.bc, simulation.approach, simulation.solver, simulation.iterations, simulation.timestep, simulation.consoleOutput, csvOutput)
function setTimestep(simulation::Simulation{T}, timestep::T)::Simulation{T} where {T}
return Simulation(simulation.grid, simulation.bc, simulation.approach, simulation.iterations, timestep, simulation.consoleOutput, simulation.csvOutput)
end