Implement midpoint approximation (#16)

* Implement MidPointApproximation

* Update tests for medium

* Restructure code for responses

* Remove old code

* Bump version

Project.toml

@@ -1,11 +1,10 @@
 name = "BoreholeNetworksSimulator"
 uuid = "b11dda51-a240-44f0-a43d-fea4e1309b86"
 authors = ["Marc Basquens <[email protected]>", "Alberto Lazzarotto <[email protected]>"]
-version = "0.1.13"
+version = "0.1.14"
 
 [deps]
 Bessels = "0e736298-9ec6-45e8-9647-e4fc86a2fe38"
-BoreholeResponseFunctions = "40984d4a-ff4b-5bb1-8527-cee9e1eb1c87"
 CSV = "336ed68f-0bac-5ca0-87d4-7b16caf5d00b"
 CoolProp = "e084ae63-2819-5025-826e-f8e611a84251"
 DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0"
@@ -34,8 +33,8 @@ Aqua = "4c88cf16-eb10-579e-8560-4a9242c79595"
 CSV = "336ed68f-0bac-5ca0-87d4-7b16caf5d00b"
 Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
 Literate = "98b081ad-f1c9-55d3-8b20-4c87d4299306"
-Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 Tables = "bd369af6-aec1-5ad0-b16a-f7cc5008161c"
+Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [targets]
 docs = ["Documenter", "Literate"]

docs/src/api.md

@@ -179,3 +179,7 @@ NoBoundary
 ```@docs
 DirichletBoundaryCondition
 ```
+
+```@docs
+AdiabaticBoundaryCondition
+```

src/BoreholeNetworksSimulator.jl

@@ -18,7 +18,6 @@ using LegendrePolynomials
 using SpecialFunctions
 using RequiredInterfaces
 
-using BoreholeResponseFunctions
 using FiniteLineSource
 
 include("utils.jl")
@@ -35,7 +34,7 @@ export Borefield, EqualBoreholesBorefield
 export Medium, GroundMedium, FlowInPorousMedium
 export Constraint, TotalHeatLoadConstraint, HeatLoadConstraint, InletTempConstraint, constant_HeatLoadConstraint, uniform_HeatLoadConstraint, constant_InletTempConstraint, uniform_InletTempConstraint
 export TimeSuperpositionMethod, ConvolutionMethod, NonHistoryMethod
-export BoundaryCondition, NoBoundary, DirichletBoundaryCondition
+export BoundaryCondition, NoBoundary, DirichletBoundaryCondition, AdiabaticBoundaryCondition
 export Approximation, MidPointApproximation, MeanApproximation
 export SimulationOptions, SimulationContainers, BoreholeNetwork
 export BoreholeOperation, Operator, SimpleOperator, operate

src/modular/approximations/MidPointApproximation.jl

@@ -3,4 +3,4 @@
 
 Option to use the midpoint temperature as an approximation for the wall temperature of a segment. 
 """
-#struct MidPointApproximation <: Approximation end
+struct MidPointApproximation <: Approximation end

src/modular/boundary_conditions/AdiabaticBoundaryCondition.jl

@@ -0,0 +1,6 @@
+"""
+    AdiabaticBoundaryCondition <: BoundaryCondition
+
+Option to enforce zero heat flow at the surface plane `z=0`.
+"""
+struct AdiabaticBoundaryCondition <: BoundaryCondition end

src/modular/core/core.jl

@@ -32,6 +32,8 @@ BoreholeOperation(::Nothing) = BoreholeOperation(BoreholeNetwork([], 0), @view o
 
 """
     struct SimulationOptions{
+                    N <: Number,
+                    Tol <: Number,
                     TSM <: TimeSuperpositionMethod,
                     C <: Constraint,
                     B <: Borefield, 
@@ -50,23 +52,28 @@ BoreholeOperation(::Nothing) = BoreholeOperation(BoreholeNetwork([], 0), @view o
         configurations::Vector{BoreholeNetwork}
         Δt
         Nt::Int
+        atol::Tol = 0.
+        rtol::Tol = sqrt(eps())
     )
 
 Specifies all the options for the simulation.
 
 - `method`: time superposition method used to compute the response. Available options: `ConvolutionMethod`, `NonHistoryMethod`.
-- `constraint`: constraint that the system must satisfy. Can be variable with time. Available options: `HeatLoadConstraint`, `InletTempConstraint`.
+- `constraint`: constraint that the system must satisfy. Can be variable with time. Available options: `HeatLoadConstraint`, `InletTempConstraint`, `TotalHeatLoadConstraint`.
 - `borefield`: describes the geometrical properties and the boreholes of the borefield on which the simulation will be performed. Available options: `EqualBoreholesBorefield`.
 - `medium`: properties of the ground where the `borefield` is places. Available options: `GroundMedium`, `FlowInPorousMedium`.
-- `boundary_condition`: boundary condition of the domain where the simulation is performed. Available options: `NoBoundary`, `DirichletBoundaryCondition`.
-- `approximation`: determines how the approximate value for each segment is computed.
+- `boundary_condition`: boundary condition of the domain where the simulation is performed. Available options: `NoBoundary`, `DirichletBoundaryCondition`, `AdiabaticBoundaryCondition`.
+- `approximation`: determines how the approximate value for each segment is computed. Available options: `MeanApproximation`, `MidPointApproximation`.
 - `fluid`: properties of the fluid flowing through the hydraulic system.
 - `configurations`: possible hydraulic topologies possible in the system, including reverse flow.
 - `Δt`: time step used in the simulation.
 - `Nt`: total amount of time steps of the simulation.
+- `atol`: absolute tolerance used for the adaptive integration methods.
+- `rtol`: relative tolerance used for the adaptive integration methods.
 """
 @with_kw struct SimulationOptions{
                     N <: Number,
+                    Tol <: Number,
                     TSM <: TimeSuperpositionMethod,
                     C <: Constraint,
                     B <: Borefield, 
@@ -90,6 +97,8 @@ Specifies all the options for the simulation.
     Tmax = Δt * Nt
     t = Δt:Δt:Tmax
     configurations::Vector{BoreholeNetwork}
+    atol::Tol = 0.
+    rtol::Tol = sqrt(eps())
 end
 
 """

src/modular/interfaces/Medium.jl

@@ -8,15 +8,11 @@ Required functions:
 - `get_λ(::Medium)`: Return the thermal conductivity of the medium.
 - `get_α(::Medium)`: Return the thermal diffusivity of the medium.
 - `get_T0(::Medium)`: Return the initial temperature of the medium.
-- `compute_response!(g, ::Medium, borefield::Borefield, boundary_condition::BoundaryCondition, t)`: 
-    Compute inplace in `g` the thermal responses between boreholes in `borefield`, 
-    imposing the boundary condition `boundary_condition`, for all times in `t`.
 """
 abstract type Medium end
 
 @required Medium begin
     get_λ(::Medium)
     get_α(::Medium)
     get_T0(::Medium)
-    compute_response!(g, ::Medium, borefield, boundary_condition, t) 
 end

src/modular/methods/ConvolutionMethod.jl

@@ -16,11 +16,12 @@ end
 ConvolutionMethod() = ConvolutionMethod(zeros(0,0,0), zeros(0,0), zeros(0))
 
 function precompute_auxiliaries!(method::ConvolutionMethod, options)
-    @unpack Nb, Nt, t, borefield, medium, boundary_condition = options
+    @unpack Nb, Nt, t, borefield, medium, boundary_condition, approximation = options
+    @unpack λ, α = medium
     method.g = zeros(Nb, Nb, Nt)
     method.q = zeros(Nb, Nt)
     method.aux = zeros(Nb)
-    compute_response!(method.g, medium, borefield, boundary_condition, t)
+    compute_response!(method.g, medium, options)
     return method
 end
 

src/modular/methods/NonHistoryMethod.jl

@@ -24,23 +24,15 @@ mutable struct NonHistoryMethod{T} <: TimeSuperpositionMethod
 end
 NonHistoryMethod(;n_disc::Int=20) = NonHistoryMethod(zeros(0, 0), zeros(0), zeros(0, 0), zeros(0), n_disc, zeros(0))
 
-function get_boreholes_distance(borefield, i, j)
-    x1, y1, D1, H1 = segment_coordinates(borefield, i)
-    x2, y2, D2, H2 = segment_coordinates(borefield, j)
-
-    σ = i == j ? get_rb(borefield, i) : sqrt((x1 - x2)^2 + (y1 - y2)^2)
-    return SegmentToSegment(σ=σ, D1=D1, D2=D2, H1=H1, H2=H2)
-end
-
-function distances(borefield, boundary_condition)
-    map = Dict{SegmentToSegment{Float64}, Int}()
+function distances(borefield, boundary_condition, approximation, medium)
+    map = Dict{setup_type(approximation, medium), Int}()
     buffers = get_buffers(boundary_condition)
 
     k = 1
     boreholes = 1:n_boreholes(borefield)
     for i in boreholes
         for j in boreholes
-            s = get_boreholes_distance(borefield, i, j)
+            s = setup(approximation, medium, borefield, i, j)
             if !haskey(map, s)
                 map[s] = k
                 k += 1
@@ -54,7 +46,7 @@ function distances(borefield, boundary_condition)
 end
 
 function precompute_auxiliaries!(method::NonHistoryMethod, options)
-    @unpack Nb, Nt, Ns, Δt, borefield, medium, boundary_condition, approximation = options
+    @unpack Nb, Nt, Ns, Δt, borefield, medium, boundary_condition, approximation, atol, rtol = options
     @unpack n_disc = method
     α = get_α(medium)
     rb = get_rb(borefield, 1) 
@@ -64,14 +56,14 @@ function precompute_auxiliaries!(method::NonHistoryMethod, options)
     b = ceil(erfcinv(ϵ / sqrt(π/Δt̃)) / sqrt(Δt̃))
 
     constants = Constants(Δt=Δt, α=α, rb=rb, kg=kg, b=b)
-    _, _, segments = quadgk_segbuf(f_guess(setup(approximation, borefield, 1, 1), constants), 0., b)
+    _, _, segments = quadgk_segbuf(f_guess(setup(approximation, medium, borefield, 1, 1), constants), 0., b, atol=atol, rtol=rtol)
     xt, w = gausslegendre(n_disc+1)  
     dps = [make_DiscretizationParameters(s.a, s.b, n_disc, xt=xt, w=w) for s in segments]
     ζ = reduce(vcat, (dp.x for dp in dps)) 
     expΔt = @. exp(-ζ^2 * Δt̃)
 
-    distances_map, quadgk_buffers = distances(borefield, boundary_condition)
-    disc_map, containers = initialize_containers(setup(approximation, borefield, 1, 1), dps)    
+    distances_map, quadgk_buffers = distances(borefield, boundary_condition, approximation, medium)
+    disc_map, containers = initialize_containers(setup(approximation, medium, borefield, 1, 1), dps)    
 
     n = length(ζ)
     w = zeros(n, Ns*Ns)
@@ -80,12 +72,12 @@ function precompute_auxiliaries!(method::NonHistoryMethod, options)
     for (key, value) in pairs(distances_map)
         for (k, dp) in enumerate(dps)
             range = (n_disc+1)*(k-1)+1:(n_disc+1)*k
-            w_buffer[range, value] .= weights(boundary_condition, key, constants, dp, containers[disc_map[k]], quadgk_buffers[value])
+            w_buffer[range, value] .= weights(boundary_condition, key, constants, dp, containers[disc_map[k]], quadgk_buffers[value], atol=atol, rtol=rtol)
         end
     end
     for i in 1:Ns
         for j in 1:Ns
-            k = distances_map[get_boreholes_distance(borefield, i, j)]
+            k = distances_map[setup(approximation, medium, borefield, i, j)]
             @views @. w[:, (i-1)*Ns+j] = w_buffer[:, k]
         end
     end
@@ -112,11 +104,10 @@ function method_coeffs!(M, method::NonHistoryMethod, options)
     @unpack borefield, medium, boundary_condition, approximation = options
     Nb = n_boreholes(borefield)
     Ns = n_segments(borefield)
-    λ = get_λ(medium)
 
     for i in 1:Ns
         for j in 1:Ns
-            M[i, 3Nb+j] = q_coef(boundary_condition, medium, method, setup(approximation, borefield, i, j), λ, (i-1)*Ns+j) 
+            M[i, 3Nb+j] = q_coef(boundary_condition, medium, method, setup(approximation, medium, borefield, i, j), (i-1)*Ns+j) 
         end
     end
 

src/modular/responses/convolution/responses.jl

@@ -0,0 +1,77 @@
+
+function compute_response!(g, medium::Medium, options) 
+    @unpack λ, α = medium
+    @unpack borefield, boundary_condition, approximation, atol, rtol, t = options
+    Nb = n_boreholes(borefield)
+
+    for (k, tt) in enumerate(t)
+        for j in 1:Nb
+            for i in 1:Nb
+                rb = get_rb(borefield, i)
+                s = setup(approximation, medium, borefield, i, j)
+                g[i, j, k] = response(boundary_condition, s, Constants(α=α, kg=λ, rb=rb), tt, atol=atol, rtol=rtol)
+            end
+        end
+    end
+end
+
+function response(::NoBoundary, setup, params::Constants, t; atol, rtol)
+    step_response(setup, params, t, atol=atol, rtol=rtol)
+end
+
+function response(::DirichletBoundaryCondition, setup, params::Constants, t; atol, rtol)
+    setup_image = image(setup)
+    Ip = step_response(setup, params, t, atol=atol, rtol=rtol)
+    In = step_response(setup_image, params, t, atol=atol, rtol=rtol)
+    Ip - In
+end
+
+function response(::AdiabaticBoundaryCondition, setup, params::Constants, t; atol, rtol)
+    setup_image = image(setup)
+    Ip = step_response(setup, params, t, atol=atol, rtol=rtol)
+    In = step_response(setup_image, params, t, atol=atol, rtol=rtol)
+    Ip + In
+end
+
+step_response(s::SegmentToPoint, params::Constants, t; atol, rtol) = stp_response(s, params, t, atol=atol, rtol=rtol)
+step_response(s::SegmentToSegment, params::Constants, t; atol, rtol) = sts_response(s, params, t, atol=atol, rtol=rtol)
+
+step_response(s::MovingSegmentToPoint, params::Constants, t; atol, rtol) = mstp_response(s, params, t, atol=atol, rtol=rtol)
+step_response(s::MovingSegmentToSegment, params::Constants, t; atol, rtol) = msts_response(s, params, t, atol=atol, rtol=rtol)
+
+function stp_response(s::SegmentToPoint, params::Constants, t; atol, rtol)
+    @unpack σ, H, D, z = s
+    @unpack α, kg = params
+    return fls_step_response(t, σ, z, D, D+H, α, kg, atol=atol, rtol=rtol)
+end
+
+function sts_response(s::SegmentToSegment, params::Constants, t; atol, rtol)
+    @unpack α, kg = params
+    params = FiniteLineSource.MeanSegToSegEvParams(s)
+    r_min, r_max = FiniteLineSource.h_mean_lims(params)
+    f(r) = FiniteLineSource.h_mean_sts(r, params) * point_step_response(t, r, α, kg)
+    x, w = FiniteLineSource.adaptive_nodes_and_weights(f, r_min, r_max, atol=atol, rtol=rtol)
+    return dot(f.(x), w)
+end
+
+function mstp_response(s::MovingSegmentToPoint, params::Constants, t; atol, rtol)
+    @unpack x, y, H, D, z, v = s
+    @unpack α, kg = params
+    return mfls_step_response(t, x, y, z, v, D, D+H, α, kg, atol=atol, rtol=rtol)
+end
+
+function msts_response(s::MovingSegmentToSegment, params::Constants, t; atol, rtol)
+    @unpack α, kg = params
+    @unpack x, v, D1, H1, D2, H2 = s
+    params = FiniteLineSource.MeanSegToSegEvParams(s)
+    r_min, r_max = FiniteLineSource.h_mean_lims(params)
+    f(r) = FiniteLineSource.h_mean_sts(r, params) * moving_point_step_response(t, r, x, v, α, kg)
+    X, W = FiniteLineSource.adaptive_nodes_and_weights(f, r_min, r_max, atol=atol, rtol=rtol)
+    return dot(f.(X), W)
+end
+
+point_step_response(t, r, α, kg) = erfc(r/(2*sqrt(t*α))) / (4*π*r*kg)   
+moving_point_step_response(t, r, x, v, α, kg) = exp(-v * (r-x)/ (2α)) * (erfc( (r-t*v) / sqrt(4t*α)) + erfc((r+t*v) / sqrt(4t*α)) * exp(v*r/α) ) / (8π*r*kg)
+
+fls_step_response(t, σ, z, a, b, α, kg; atol, rtol) = quadgk(ζ -> point_step_response(t, sqrt(σ^2 + (z-ζ)^2), α, kg), a, b, atol=atol, rtol=rtol)[1]
+mfls_step_response(t, x, y, z, v, a, b, α, kg; atol, rtol) = quadgk(ζ -> moving_point_step_response(t, sqrt(x^2 + y^2 + (z-ζ)^2), x, v, α, kg), a, b, atol=atol, rtol=rtol)[1]