Add boundary conditions

marcbasquensmunoz · Jul 19, 2024 · 4b06176 · 4b06176
1 parent 5c02482
commit 4b06176
Show file tree

Hide file tree

Showing 13 changed files with 83 additions and 37 deletions.
diff --git a/Project.toml b/Project.toml
@@ -32,4 +32,4 @@ Aqua = "4c88cf16-eb10-579e-8560-4a9242c79595"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [targets]
-test = ["Test", "Aqua"]
+test = ["Test", "Aqua"]
diff --git a/README.md b/README.md
@@ -15,7 +15,7 @@ Simply clone the project using e.g. git clone, cd to the project directory, ente
 ## Example
 The example considers a borehole field of 48 boreholes connected according to scheme utilized in the installation in Braedsturp, Denmark. 
 
-<img src="./examples/plots/results/configuration.png" width="400" height="400" />
+<img src="./examples/old/results/configuration.png" width="400" height="400" />
 
 ## Run the example
 ```
@@ -45,13 +45,13 @@ julia > include("examples/plots/sim1.jl")
 
 Inlet borehole temperatures and heat flows for boreholes along two branches in the borehole field. The time series are color coded according to the previous configuration plot above. In addition to the inlet temperature, the output temperature from the branch (grey dot), and the mean output temperature from the field (black dot) are displayed.
 
-<img src="./examples/plots/results/sym1/branch1_test1.png" width="600" height="300" />
-<img src="./examples/plots/results/sym1/branch2_test1.png" width="600" height="300" />
+<img src="./examples/old/results/sym1/branch1_test1.png" width="600" height="300" />
+<img src="./examples/old/results/sym1/branch2_test1.png" width="600" height="300" />
 
 
 Finally we can display the heatmap of the temperature field in the borehole region during the 10th year of operation
 
-<img src="./examples/plots/results/sym1/heatmap_test1.png" width="600" height="300" />
+<img src="./examples/old/results/sym1/heatmap_test1.png" width="600" height="300" />
 
 
 ## Running the code in Python

diff --git a/examples/modular/test.jl b/examples/modular/test.jl
@@ -2,27 +2,29 @@ using BoreholeNetworksSimulator
 using FiniteLineSource
 
 Δt = 8760*3600/12.
-Nt = 10
-network = [
+Nt = 1
+configurations = [
     BoreholeNetwork([[1]])
 ]
-Nbr = length(network)
 
 α = 1e-6
 λ = 3.
 rb = 0.1
 
+D = 0.
+H = 1.
+
 positions = [(0., i) for i in 0:0]
 
 T0 = 0.
-borehole = SingleUPipeBorehole(H=1., D=0., rb=rb)
+borehole = SingleUPipeBorehole(H=H, D=D, rb=rb)
 medium = GroundMedium(α=α, λ=λ)
 borefield = EqualBoreholesBorefield(borehole_prototype=borehole, positions=positions, medium=medium, T0 = T0)
-constraint = InletTempConstraint(10*ones(1))
-#constraint = HeatLoadConstraint([1.])
+#constraint = InletTempConstraint(10*ones(1))
+constraint = HeatLoadConstraint([1.])
 fluid = Fluid(cpf = 4182., name = "INCOMP::MEA-20%")
 function operator(i, Tin, Tout, Tb, q)
-    BoreholeOperation(network[1], 0.1 .* ones(1))
+    BoreholeOperation(configurations[1], 0.1 .* ones(1))
 end
 
 
@@ -31,6 +33,7 @@ options_c = SimulationOptions(
     constraint = constraint,
     borefield = borefield,
     fluid = fluid,
+   #boundary_condition=NoBoundary(),
     Δt = Δt,
     Nt = Nt
 )
@@ -45,6 +48,7 @@ options_nh = SimulationOptions(
     constraint = constraint,
     borefield = borefield,
     fluid = fluid,
+    #boundary_condition=NoBoundary(),
     Δt = Δt,
     Nt = Nt
 )
@@ -57,12 +61,18 @@ X1-X2
 
 
 ### FLS
-q = [1 for i=1:10]
+q = [1 for i=1:Nt]
 I = zeros(length(q))
-
-setup = FiniteLineSource.SegmentToSegment(D1=0., H1=1., D2=0., H2=1., σ=.1)
 params = FiniteLineSource.Constants(Δt=Δt, rb=rb, b=2.)
+
+setup = FiniteLineSource.SegmentToSegment(D1=D, H1=H, D2=D, H2=H, σ=rb)
 precomp = FiniteLineSource.precompute_parameters(setup, params=params)
 compute_integral_throught_history!(setup, I=I, q=q, precomp=precomp, params=params)
 
-BoreholeNetworksSimulator.sts(setup, params, tstep)
+I2 = zeros(length(q))
+setup2 = FiniteLineSource.SegmentToSegment(D1=-D, H1=-H, D2=D, H2=H, σ=rb)
+precomp2 = FiniteLineSource.precompute_parameters(setup2, params=params)
+compute_integral_throught_history!(setup2, I=I2, q=q, precomp=precomp2, params=params)
+
+I-I2
+
diff --git a/src/BoreholeNetworksSimulator.jl b/src/BoreholeNetworksSimulator.jl
@@ -30,6 +30,7 @@ export Borefield, EqualBoreholesBorefield
 export Medium, GroundMedium, FlowInPorousMedium
 export Constraint, HeatLoadConstraint, InletTempConstraint
 export TimeSuperpositionMethod, ConvolutionMethod, NonHistoryMethod
+export BoundaryCondition, NoBoundary, DirichletBoundaryCondition
 export SimulationOptions, SimulationContainers, BoreholeOperation, BoreholeNetwork, Fluid
 export simulate, compute_parameters, load_cache!, save_cache, initialize
 

diff --git a/src/modular/boundary_conditions/DirichletBoundaryCondition.jl b/src/modular/boundary_conditions/DirichletBoundaryCondition.jl
@@ -0,0 +1 @@
+struct DirichletBoundaryCondition <: BoundaryCondition end
diff --git a/src/modular/boundary_conditions/NoBoundary.jl b/src/modular/boundary_conditions/NoBoundary.jl
@@ -0,0 +1 @@
+struct NoBoundary <: BoundaryCondition end
diff --git a/src/modular/core/core.jl b/src/modular/core/core.jl
@@ -2,7 +2,6 @@
 @with_kw struct BoreholeOperation
     network         
     mass_flows:: Vector          # Mass flow for each branch in network
-    #cpf::T = 4182.                  # specific heat capacity
 end
 BoreholeOperation(::Nothing) = BoreholeOperation([[]], [0.])
 
@@ -23,6 +22,7 @@ first_boreholes(network::BoreholeNetwork) = map(first, network.branches)
     constraint::Constraint
     borefield::Borefield
     fluid::Fluid
+    boundary_condition::BoundaryCondition = DirichletBoundaryCondition()
     Δt
     Nt
     Nb = n_boreholes(borefield)

diff --git a/src/modular/interfaces/BoundaryCondition.jl b/src/modular/interfaces/BoundaryCondition.jl
@@ -0,0 +1 @@
+abstract type BoundaryCondition end
diff --git a/src/modular/methods/ConvolutionMethod.jl b/src/modular/methods/ConvolutionMethod.jl
@@ -6,10 +6,10 @@ end
 ConvolutionMethod() = ConvolutionMethod(zeros(0,0,0), zeros(0,0))
 
 function precompute_auxiliaries!(model::ConvolutionMethod; options::SimulationOptions)
-    @unpack Nb, Nt, t, borefield = options
+    @unpack Nb, Nt, t, borefield, boundary_condition = options
     model.g = zeros(Nb, Nb, Nt)
     model.q = zeros(Nb, Nt)
-    compute_response!(model.g, borefield.medium, borefield, t)
+    compute_response!(model.g, borefield.medium, borefield, boundary_condition, t)
     return model
 end
 
@@ -18,7 +18,7 @@ function update_auxiliaries!(method::ConvolutionMethod, X, borefield::Borefield,
     method.q[:, step] = @view X[3Nb+1:end, step] 
 end
 
-function method_coeffs!(M, method::ConvolutionMethod, borefield::Borefield)
+function method_coeffs!(M, method::ConvolutionMethod, borefield::Borefield, ::BoundaryCondition)
     Nb = n_boreholes(borefield)
     Ns = n_segments(borefield)
     M[1:Ns, 3Nb+1:3Nb+Ns] = @view method.g[:,:,1]

diff --git a/src/modular/methods/NonHistoryMethod.jl b/src/modular/methods/NonHistoryMethod.jl
@@ -1,3 +1,5 @@
+using .FiniteLineSource: SegmentToSegment, Constants, adaptive_gk_segments, discretization_parameters, f_guess, precompute_coefficients
+
 mutable struct NonHistoryMethod{T} <: TimeSuperpositionMethod 
     F::Matrix{T}
     ζ::Vector{T}
@@ -6,10 +8,11 @@ mutable struct NonHistoryMethod{T} <: TimeSuperpositionMethod
 end
 NonHistoryMethod() = NonHistoryMethod(zeros(0, 0), zeros(0), zeros(0, 0), zeros(0))
 
-FiniteLineSource.SegmentToSegment(s::MeanSegToSegEvParams) = FiniteLineSource.SegmentToSegment(D1=s.D1, H1=s.H1, D2=s.D2, H2=s.H2, σ=s.σ)
+SegmentToSegment(s::MeanSegToSegEvParams) = SegmentToSegment(D1=s.D1, H1=s.H1, D2=s.D2, H2=s.H2, σ=s.σ)
+image(s::SegmentToSegment) = SegmentToSegment(D1=-s.D1, H1=-s.H1, D2=s.D2, H2=s.H2, σ=s.σ)
 
 function precompute_auxiliaries!(model::NonHistoryMethod; options::SimulationOptions)
-    @unpack Nb, Nt, Ns, Δt, borefield = options
+    @unpack Nb, Nt, Ns, Δt, borefield, boundary_condition = options
     b = 2.
     α = get_α(borefield.medium)
     rb = get_rb(borefield, 1) 
@@ -18,9 +21,9 @@ function precompute_auxiliaries!(model::NonHistoryMethod; options::SimulationOpt
 
     n_segment = 20
 
-    constants = FiniteLineSource.Constants(Δt=Δt, α=α, rb=rb, kg=kg, b=b)
-    segments = FiniteLineSource.adaptive_gk_segments(FiniteLineSource.f_guess(FiniteLineSource.SegmentToSegment(get_sts(borefield, 1, 1)), constants), 0., b)
-    dps = @views [FiniteLineSource.discretization_parameters(s.a, s.b, n_segment) for s in segments]
+    constants = Constants(Δt=Δt, α=α, rb=rb, kg=kg, b=b)
+    segments = adaptive_gk_segments(f_guess(SegmentToSegment(get_sts(borefield, 1, 1)), constants), 0., b)
+    dps = @views [discretization_parameters(s.a, s.b, n_segment) for s in segments]
     ζ = reduce(vcat, (dp.x for dp in dps)) 
     expΔt = @. exp(-ζ^2 * Δt̃)
 
@@ -30,8 +33,8 @@ function precompute_auxiliaries!(model::NonHistoryMethod; options::SimulationOpt
     for i in 1:Ns
         for j in 1:Ns
             rb = get_rb(borefield, (i-1)*Ns+j) # Get the right rb
-            sts = get_sts(borefield, i, j)
-            w[:, (i-1)*Ns+j] = reduce(vcat, [FiniteLineSource.precompute_coefficients(FiniteLineSource.SegmentToSegment(sts), params=constants, dp=dp) for dp in dps])
+            sts = SegmentToSegment(get_sts(borefield, i, j))
+            w[:, (i-1)*Ns+j] = reduce(vcat, [coefficients_sts(boundary_condition, sts, constants, dp) for dp in dps])
         end
     end
 
@@ -43,6 +46,17 @@ function precompute_auxiliaries!(model::NonHistoryMethod; options::SimulationOpt
     model.expΔt = expΔt[perm]
 end
 
+function coefficients_sts(::NoBoundary, sts::SegmentToSegment, params::Constants, dp)
+    precompute_coefficients(sts, params=params, dp=dp)
+end
+
+function coefficients_sts(::DirichletBoundaryCondition, sts::SegmentToSegment, params::Constants, dp)
+    image_sts = image(sts)
+    w1 = precompute_coefficients(sts, params=params, dp=dp)
+    w2 = precompute_coefficients(image_sts, params=params, dp=dp)
+    w1 - w2
+end
+
 function get_sts(borefield::Borefield, i, j)
     xi, yi, Di, Hi = segment_coordinates(borefield, i)
     xj, yj, Dj, Hj = segment_coordinates(borefield, j)
@@ -59,16 +73,15 @@ function update_auxiliaries!(method::NonHistoryMethod, X, borefield::Borefield,
     end
 end
 
-function method_coeffs!(M, method::NonHistoryMethod, borefield::Borefield)
+function method_coeffs!(M, method::NonHistoryMethod, borefield::Borefield, boundary_condition::BoundaryCondition)
     Nb = n_boreholes(borefield)
     Ns = n_segments(borefield)
     λ = get_λ(borefield.medium)
 
     for i in 1:Ns
         for j in 1:Ns
-            sts = get_sts(borefield, i, j)
-            rb = get_rb(borefield, i) 
-            M[i, 3Nb+j] = - q_coef(borefield.medium, method, sts, λ, (i-1)*Ns+j)
+            sts = SegmentToSegment(get_sts(borefield, i, j))
+            M[i, 3Nb+j] = - q_coef(boundary_condition, borefield.medium, method, sts, λ, (i-1)*Ns+j) 
         end
     end
 
@@ -88,6 +101,14 @@ function method_b!(b, method::NonHistoryMethod, borefield::Borefield, step)
     end
 end
 
+function q_coef(::NoBoundary, medium, method, sts, λ, i)
+    q_coef(medium, method, sts, λ, i)
+end
+
+function q_coef(::DirichletBoundaryCondition, medium, method, sts, λ, i)
+    q_coef(medium, method, sts, λ, i) - q_coef(medium, method, image(sts), λ, i)
+end
+
 function q_coef(::GroundMedium, method, sts, λ, i)
     @unpack D1, H1, D2, H2, σ = sts
     @unpack expΔt, w, ζ = method

diff --git a/src/modular/methods/ground_response.jl b/src/modular/methods/ground_response.jl
@@ -1,5 +1,5 @@
 
-function compute_response!(g, medium::GroundMedium, borefield::Borefield, t) 
+function compute_response!(g, medium::GroundMedium, borefield::Borefield, boundary_condition::BoundaryCondition, t) 
     @unpack λ, α = medium
     Ns = n_segments(borefield)
 
@@ -13,7 +13,7 @@ function compute_response!(g, medium::GroundMedium, borefield::Borefield, t)
                 rb = get_rb(borefield, i)
                 σ = i == j ? rb : sqrt((x2-x1)^2 + (y2-y1)^2)
                 s = SegmentToSegment(D1=D1, H1=H1, D2=D2, H2=H2, σ=σ)
-                g[i, j, k] = sts(s, Constants(α=α, kg=λ, rb=rb), tt)
+                g[i, j, k] = sts(boundary_condition, s, Constants(α=α, kg=λ, rb=rb), tt)
             end
         end
     end
@@ -41,11 +41,22 @@ function fls(t, x, y, z, H, D, α, kg, atol = 1e-8)
     return fls_step_response(t, σ, z, D, D+H, α, kg, atol=atol)
 end
 
-function sts(s:: SegmentToSegment, params::Constants, t)
+function sts_response(s::SegmentToSegment, params::Constants, t)
     @unpack α, kg = params
     params = FiniteLineSource.MeanSegToSegEvParams(s)
     h_mean_sts, r_min, r_max = FiniteLineSource.mean_sts_evaluation(params)
     f(r) = h_mean_sts(r) * point_step_response(t, r, α, kg)
     x, w = FiniteLineSource.adaptive_gk(f, r_min, r_max)
     dot(f.(x), w)
 end
+
+function sts(::NoBoundary, s::SegmentToSegment, params::Constants, t)
+    sts_response(s, params, t)
+end
+
+function sts(::DirichletBoundaryCondition, s::SegmentToSegment, params::Constants, t)
+    s_image = SegmentToSegment(D1=-s.D1, H1=-s.H1, D2=s.D2, H2=s.H2, σ=s.σ)
+    Ip = sts_response(s, params, t)
+    In = sts_response(s_image, params, t)
+    Ip - In
+end
diff --git a/src/modular/methods/ground_water_response.jl b/src/modular/methods/ground_water_response.jl
@@ -1,5 +1,5 @@
 
-function compute_response!(g, medium::FlowInPorousMedium, borefield::Borefield, t) 
+function compute_response!(g, medium::FlowInPorousMedium, borefield::Borefield, boundary_condition::BoundaryCondition, t) 
     @unpack θ, λ, vt, α = medium
     Ns = segment_amount(borefield)
 

diff --git a/src/modular/simulate.jl b/src/modular/simulate.jl
@@ -11,7 +11,7 @@
 #   (4)  Heat balance equations
 
 function simulate(;operator, options::SimulationOptions, containers::SimulationContainers)
-    @unpack method, constraint, borefield, fluid = options
+    @unpack method, constraint, borefield, fluid, boundary_condition = options
     @unpack Nb, Ns, Nt, Ts = options
     @unpack M, b, X = containers 
 
@@ -45,7 +45,7 @@ function simulate(;operator, options::SimulationOptions, containers::SimulationC
         end
         @views constraints_coeffs!(M[constraints_eqs, :], constraint, operation)
         if i == Ts
-            @views method_coeffs!(M[method_eqs, :], method, borefield)
+            @views method_coeffs!(M[method_eqs, :], method, borefield, boundary_condition)
         end
         if last_operation.mass_flows != operation.mass_flows
             @views heat_balance_coeffs!(M[balance_eqs, :], borefield, operation, fluid)
Original file line number	Diff line number	Diff line change
		@@ -0,0 +1 @@
		struct DirichletBoundaryCondition <: BoundaryCondition end