Optimizations (#24)

* Accelerate midpoint for convolution

* Move code

* Make convolution sts fast

* Rename Adiabatic to Neumann

* Rewrite integrals for the step responses in the convolution case

* Fix allocations bug

* Remove unnecessary dependencies

docs/src/api.md

@@ -228,5 +228,5 @@ DirichletBoundaryCondition
 ```
 
 ```@docs
-AdiabaticBoundaryCondition
+NeumannBoundaryCondition
 ```

src/BoreholeNetworksSimulator.jl

@@ -35,7 +35,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, AdiabaticBoundaryCondition
+export BoundaryCondition, NoBoundary, DirichletBoundaryCondition, NeumannBoundaryCondition
 export Approximation, MidPointApproximation, MeanApproximation
 export SimulationOptions, SimulationContainers
 export BoreholeNetwork, Valve, BoreholeOperation

src/modular/borefields/EqualBoreholesBorefield.jl

@@ -10,10 +10,11 @@ Note that the length of `positions` determines the amount of boreholes in the fi
 @with_kw struct EqualBoreholesBorefield{T <: Borehole, S <: Real} <: Borefield
     borehole_prototype::T
     positions::Vector{Tuple{S, S}}
+    Nb::Int = length(positions)
 end
 
-n_boreholes(bf::EqualBoreholesBorefield) = length(bf.positions)
-n_segments(bf::EqualBoreholesBorefield) = length(bf.positions)
+n_boreholes(bf::EqualBoreholesBorefield) = bf.Nb
+n_segments(bf::EqualBoreholesBorefield) = bf.Nb
 get_H(bf::EqualBoreholesBorefield, i) = get_H(bf.borehole_prototype)
 get_h(bf::EqualBoreholesBorefield, i) = get_h(bf.borehole_prototype)
 get_rb(bf::EqualBoreholesBorefield, i) = get_rb(bf.borehole_prototype)

src/modular/boundary_conditions/AdiabaticBoundaryCondition.jl

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

src/modular/core/core.jl

@@ -31,7 +31,7 @@ Specifies all the options for the simulation.
 - `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`, `AdiabaticBoundaryCondition`.
+- `boundary_condition`: boundary condition of the domain where the simulation is performed. Available options: `NoBoundary`, `DirichletBoundaryCondition`, `NeumannBoundaryCondition`.
 - `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.
@@ -63,8 +63,8 @@ Specifies all the options for the simulation.
     Nb::Int = n_boreholes(borefield)
     Ns::Int = n_segments(borefield)
     Ts::Int = 1
-    Tmax = Δt * Nt
-    t = Δt:Δt:Tmax
+    Tmax::N = Δt * Nt
+    t::Vector{N} = collect(Δt:Δt:Tmax)
     configurations::Vector{BoreholeNetwork}
     atol::Tol = 0.
     rtol::Tol = sqrt(eps())

src/modular/methods/NonHistoryMethod.jl

@@ -57,8 +57,8 @@ function precompute_auxiliaries!(method::NonHistoryMethod, options)
 
     constants = Constants(Δt=Δt, α=α, rb=rb, kg=kg, b=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]
+    xt, wt = gausslegendre(n_disc+1)  
+    dps = [make_DiscretizationParameters(s.a, s.b, n_disc, xt=xt, w=wt) for s in segments]
     ζ = reduce(vcat, (dp.x for dp in dps)) 
     expΔt = @. exp(-ζ^2 * Δt̃)
 
@@ -78,7 +78,7 @@ function precompute_auxiliaries!(method::NonHistoryMethod, options)
     for i in 1:Ns
         for j in 1:Ns
             k = distances_map[setup(approximation, medium, borefield, i, j)]
-            @views @. w[:, (i-1)*Ns+j] = w_buffer[:, k]
+            @inbounds @views @. w[:, (i-1)*Ns+j] = w_buffer[:, k]
         end
     end
 

src/modular/responses/convolution/responses.jl

@@ -15,63 +15,40 @@ function compute_response!(g, medium::Medium, options)
     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
+# TODO: Might be able to improve performance by explicitly writing the integrand in each case
+integrand(::NoBoundary, setup, s) = integrand(setup, s)
+integrand(::DirichletBoundaryCondition, setup, s) = integrand(setup, s) - integrand(image(setup), s)
+integrand(::NeumannBoundaryCondition, setup, s) = integrand(setup, s) + integrand(image(setup), s)
 
-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
+integrand(setup::SegmentToPoint, s) = I_stp(s, setup.D, setup.H, setup.z)
+integrand(setup::SegmentToSegment, s) = I_sts(s, setup.D1, setup.D2, setup.H1, setup.H2) 
+integrand(setup::MovingSegmentToPoint, s) = I_stp(s, setup.D, setup.H, setup.z)
+integrand(setup::MovingSegmentToSegment, s) = I_sts(s, setup.D1, setup.D2, setup.H1, setup.H2) 
 
-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
+function response(boundary_condition, setup::SegmentToPoint, params::Constants, t; atol, rtol)
     @unpack α, kg = params
-    return fls_step_response(t, σ, z, D, D+H, α, kg, atol=atol, rtol=rtol)
+    @unpack σ, H, D, z = setup
+    1 / (4π * kg) * quadgk(s -> exp(-σ^2*s^2) / s * integrand(boundary_condition, setup, s), 1/sqrt(4*α*t), Inf, rtol = rtol, atol = atol)[1]
 end
 
-function sts_response(s::SegmentToSegment, params::Constants, t; atol, rtol)
+function response(boundary_condition, setup::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)
+    @unpack σ, H2 = setup
+    1 / (4π * H2 * kg) * quadgk(s -> exp(-σ^2*s^2) / s^2 * integrand(boundary_condition, setup, s), 1/sqrt(4*α*t), Inf, rtol = rtol, atol = atol)[1]
 end
 
-function mstp_response(s::MovingSegmentToPoint, params::Constants, t; atol, rtol)
-    @unpack x, y, H, D, z, v = s
+function response(boundary_condition, setup::MovingSegmentToPoint, params::Constants, t; atol, rtol)
     @unpack α, kg = params
-    return mfls_step_response(t, x, y, z, v, D, D+H, α, kg, atol=atol, rtol=rtol)
+    @unpack x, y, H, D, z, v = setup
+    1 / (4π * kg) * quadgk(s -> exp(-((x-v/(4α*s^2))^2 + y^2)*s^2) / s * integrand(boundary_condition, setup, s), 1/sqrt(4*α*t), Inf, rtol = rtol, atol = atol)[1]
 end
 
-function msts_response(s::MovingSegmentToSegment, params::Constants, t; atol, rtol)
+function response(boundary_condition, setup::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)
+    @unpack x, y, v, D1, H1, D2, H2 = setup
+    1 / (4π * H2 * kg) * quadgk(s -> exp(-s^2*((x - v/(4α*s^2))^2 + y^2)) / s^2 * integrand(boundary_condition, setup, s), 1/sqrt(4*α*t), Inf, rtol = rtol, atol = atol)[1]
 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]
+ierf(x) = x*erf(x) - (1 - exp(-x^2)) / sqrt(π)
+I_sts(s, D1, D2, H1, H2) = ierf((D2 - D1 + H2)*s) + ierf((D2 - D1 - H1)*s) - ierf((D2 - D1)*s) - ierf((D2 - D1 + H2 - H1)*s)
+I_stp(s, D, H, z) = erf(s * (z-D)) - erf(s * (z-D-H))

src/modular/responses/nonhistory/buffers.jl

@@ -15,7 +15,7 @@ ImageQuadGKBuffer(s, rb::T) where {T <: Number} = ImageQuadGKBuffer(Buffer(initi
 
 get_buffers(::NoBoundary) = Vector{SimpleQuadGKBuffer}()
 get_buffers(::DirichletBoundaryCondition) = Vector{ImageQuadGKBuffer}()
-get_buffers(::AdiabaticBoundaryCondition) = Vector{ImageQuadGKBuffer}()
+get_buffers(::NeumannBoundaryCondition) = Vector{ImageQuadGKBuffer}()
 
 function add_buffer!(buffers, ::NoBoundary, s, rb)
     push!(buffers, SimpleQuadGKBuffer(s, rb))
@@ -25,6 +25,6 @@ function add_buffer!(buffers, ::DirichletBoundaryCondition, s, rb)
     push!(buffers, ImageQuadGKBuffer(s, rb))
 end
 
-function add_buffer!(buffers, ::AdiabaticBoundaryCondition, s, rb)
+function add_buffer!(buffers, ::NeumannBoundaryCondition, s, rb)
     push!(buffers, ImageQuadGKBuffer(s, rb))
 end

src/modular/responses/nonhistory/q_coef.jl

@@ -1,16 +1,16 @@
 
 function q_coef(::NoBoundary, medium, method, setup, i)
-    constant_integral(medium, method, setup, i) + constant_coef(method, i)
+    constant_integral(medium, setup, i) + constant_coef(method, i)
 end
 
 function q_coef(::DirichletBoundaryCondition, medium, method, setup, i)
     @unpack expΔt, w, ζ = method
-    constant_integral(medium, method, setup, i) - constant_integral(medium, method, image(setup), i) + constant_coef(method, i)
+    constant_integral(medium, setup, i) - constant_integral(medium, image(setup), i) + constant_coef(method, i)
 end
 
-function q_coef(::AdiabaticBoundaryCondition, medium, method, setup, i)
+function q_coef(::NeumannBoundaryCondition, medium, method, setup, i)
     @unpack expΔt, w, ζ = method
-    constant_integral(medium, method, setup, i) + constant_integral(medium, method, image(setup), i) + constant_coef(method, i)
+    constant_integral(medium, setup, i) + constant_integral(medium, image(setup), i) + constant_coef(method, i)
 end
 
 function constant_coef(method::NonHistoryMethod, i)
@@ -19,17 +19,16 @@ function constant_coef(method::NonHistoryMethod, i)
     @views -dot(w[:, i], aux)
 end
 
-function constant_integral(medium::GroundMedium, method, setup::SegmentToSegment, i)
+function constant_integral(medium::GroundMedium, setup::SegmentToSegment, i)
     @unpack D1, H1, D2, H2, σ = setup
     @unpack λ = medium
 
     β(d) = sqrt(σ^2 + d^2) + d*log(sqrt(σ^2 + d^2) - d)
     1/(4π*λ*H2) * (β(D1+H1-D2-H2) + β(D1-D2) - β(D1+H1-D2) - β(D1-D2-H2))
 end
 
-function constant_integral(medium::GroundMedium, method, setup::SegmentToPoint, i)
+function constant_integral(medium::GroundMedium, setup::SegmentToPoint, i)
     @unpack D, H, z, σ = setup
-    @unpack expΔt, w, ζ = method
     @unpack λ = medium
 
     1/(4π*λ) * log((z-D+sqrt(σ^2+(z-D)^2)) / (z-D-H+sqrt(σ^2+(z-D-H)^2)))

src/modular/responses/nonhistory/weights.jl

@@ -10,7 +10,7 @@ function weights(::DirichletBoundaryCondition, setup, params::Constants, dp, con
     w1 - w2
 end
 
-function weights(::AdiabaticBoundaryCondition, setup, params::Constants, dp, containers, buffer::ImageQuadGKBuffer; atol, rtol)
+function weights(::NeumannBoundaryCondition, setup, params::Constants, dp, containers, buffer::ImageQuadGKBuffer; atol, rtol)
     image_setup = image(setup)
     w1 = precompute_coefficients(setup, params=params, dp=dp, containers=containers, buffer=buffer.buffer, atol=atol, rtol=rtol)
     w2 = precompute_coefficients(image_setup, params=params, dp=dp, containers=containers, buffer=buffer.image_buffer, atol=atol, rtol=rtol)

test/api/api.jl

@@ -14,7 +14,7 @@ end
 end
 
 @testset "test_Boundary_Conditions" begin
-    @test AdiabaticBoundaryCondition() isa Any
+    @test NeumannBoundaryCondition() isa Any
     @test DirichletBoundaryCondition() isa Any
     @test NoBoundary() isa Any
 end