diff --git a/docs/make.jl b/docs/make.jl index 0e3db72..9a6f098 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -68,8 +68,7 @@ end function example_info_footer(subdir, exname) return "\n\n# ## Example on GitHub\n"* "# If you would like to run this example yourself, it can be downloaded from "* - "the Fimbul.jl GitHub repository [as a script](https://github.com/sintefmath/Fimbul.jl/blob/main/examples/$subdir/$exname.jl), "* - "or as a [Jupyter Notebook](https://github.com/sintefmath/Fimbul.jl/blob/gh-pages/dev/final_site/notebooks/$subdir/$exname.ipynb)" + "the Fimbul.jl GitHub repository [as a script](https://github.com/sintefmath/Fimbul.jl/blob/main/examples/$subdir/$exname.jl)." end function update_footer(content, subdir, exname) @@ -193,8 +192,12 @@ function build_fimbul_docs( "man/cases/cases.md", "man/cases/utils.md", ], + "References" => [ + "Bibliography" => "extras/refs.md" + ], ], "Examples" => examples_markdown, + ] # for (k, subpages) in build_pages # println("$k") diff --git a/docs/src/extras/refs.md b/docs/src/extras/refs.md new file mode 100644 index 0000000..e6f88f3 --- /dev/null +++ b/docs/src/extras/refs.md @@ -0,0 +1,2 @@ +```@bibliography +``` \ No newline at end of file diff --git a/examples/analytical/analytical_1d.jl b/examples/analytical/analytical_1d.jl index 54797ef..d85b47a 100644 --- a/examples/analytical/analytical_1d.jl +++ b/examples/analytical/analytical_1d.jl @@ -1,7 +1,10 @@ # # Heat equation in 1D -# This example goes through the the classical solution of the heat equation in -# 1D, and compares the analytical solution to the numerical solution obtained -# using JutulDary to verify the numerical scheme. +# This example demonstrates the classical solution of the heat equation in 1D +# and compares the analytical solution to the numerical solution obtained using +# JutulDarcy to verify the accuracy and convergence properties of the numerical +# scheme. + +# Add required modules and define convenience functions for unit conversion using Jutul, JutulDarcy, Fimbul using LinearAlgebra using HYPRE @@ -12,20 +15,21 @@ to_celsius = T -> convert_from_si(T, :Celsius) # ## The 1D heat equation # We consider a piece of solid rock with length $L$ and constant thermal -# donductivity $\lambda$, and heat capacity $C_p$ and density $\rho$, for -# which conservation of energy takes the form of the archetypical heat equation: +# conductivity $\lambda$, heat capacity $C_p$ and density $\rho$. Conservation +# of energy in this system takes the form of the archetypical heat equation: # # ``\dot{T} = \alpha \Delta T``, ``\alpha = \frac{\lambda}{\rho C_p}`` # -# where $T$ is the temperature and $\alpha$ is called the thermal diffusivity. -# Imposing equal and fixed boundary conditions $T(0, t) = T(L, t) = T_b$ and -# initial conditions T(0, x) = T_0(x), the solution to this equation can be -# found in any textbook on partial differential equations, and is given by +# where $T$ is the temperature and $\alpha$ is the thermal diffusivity. +# For this problem, we impose fixed boundary conditions $T(0, t) = T(L, t) = T_b$ +# and initial conditions $T(x, 0) = T_0(x)$. The analytical solution to this +# boundary value problem can be found in any textbook on partial differential +# equations, and is given by the Fourier series: # -# ``T(x, t) = T_b + \sum_{k = 1}^{\infty} C_k \exp\left(-\alpha \big(\frac{(k\pi}{L}\big)^2 t\right) \sin\big(\frac{k\pi}{L} x\big)`` +# ``T(x, t) = T_b + \sum_{k = 1}^{\infty} C_k \exp\left(-\alpha \big(\frac{k\pi}{L}\big)^2 t\right) \sin\big(\frac{k\pi}{L} x\big)`` # -# where the coefficients $C_k$ are determined by the initial condition and -# boundary conditions: +# where the Fourier coefficients $C_k$ are determined by the initial condition +# and boundary conditions through the orthogonality of sine functions: # # ``C_k = \frac{2}{L} \int_0^L (T_0(x) - T_b) \sin\big(\frac{n\pi}{L} x) dx`` # @@ -34,7 +38,7 @@ to_celsius = T -> convert_from_si(T, :Celsius) # We first consider a simple initial condition where the initial temperature # profile takes the form of a sine curve, # -# ``T_0(x) = T_\max\sin\big(\frac{2\pi}{L} x\big)`` +# ``T_0(x) = T_\max\sin\big(\frac{\pi}{L} x\big)`` L = 100.0 T_b = to_kelvin(10.0) @@ -45,10 +49,10 @@ case, sol, x, t = analytical_1d( results = simulate_reservoir(case, info_level=0) -# ### Plot the temperature profile -# We set up a function for plotting the numerical and analytical temperature -# profiles at a selected number of timesteps from the initial time up to a -# fraction of the total time. +# ### Visualize temperature evolution +# We set up a visualization function that plots both numerical and analytical +# temperature profiles at selected timesteps. This allows us to visually assess +# the agreement between the two solutions as thermal diffusion progresses. function plot_temperature_1d(case, sol_n, sol_a, x_n, t_n, n) fig = Figure(size=(800, 600), fontsize=20) ax = Axis(fig[1, 1]; xlabel="Distance (m)", ylabel="Temperature (°C)") @@ -82,8 +86,9 @@ end plot_temperature_1d(case, results, sol, x, t, 10) # ## Piecewise constant initial conditions -# We now consider a more complex initial condition where the initial temperature -# profile is piecewise constant, with four different constant values. +# We now consider a more challenging initial condition where the initial +# temperature profile is piecewise constant, with four different constant +# values. T_0 = x -> to_kelvin(100.0) .* (x < 25) + to_kelvin(20.0) .* (25 <= x < 50) + @@ -92,7 +97,7 @@ T_0 = x -> # We can still compute the analytical solution using the formula above, but the # sum will now be an infinite series. The function `analytical_1d` handles this -# pragmatically by cutting off the series when the contribution of the next term +# pragmatically by truncating the series when the contribution of the next term # is less than a 1e-6. case, sol, x, t = analytical_1d( L=L, temperature_boundary=T_b, initial_condition=T_0, @@ -101,13 +106,12 @@ case, sol, x, t = analytical_1d( results = simulate_reservoir(case, info_level=0) plot_temperature_1d(case, results, sol, x, t, 10) -# ## Convergence study -# Next, we perform a convergence study by simulating the same problem with -# increasing number of cells and timesteps, and comparing the numerical and -# analytical solutions. The default two-point flux apporximation (TPFA) scheme -# used in JutulDarcy reduces to a central finite difference scheme for the heat -# equation, which is second order accurate in space, whereas Backward Euler is -# first order accurate in time. +# ## Numerical convergence analysis +# We perform a comprehensive convergence study by systematically refining the +# spatial and temporal discretization. The default two-point flux approximation +# (TPFA) scheme used in JutulDarcy reduces to a central finite difference scheme +# for the heat equation, which is second-order accurate in space, whereas the +# implicit Backward Euler time integration is first-order accurate in time. # We use the same initial conditions as in the first example above T_0 = x -> to_kelvin(90.0) .* sin(π * x / L) .+ T_b; @@ -159,13 +163,12 @@ function convergence_1d(setup_fn, type=:space; Nx=2 .^ (range(3, 6)), Nt=2 .^ (r end -# ### Convergence in space -# We simulate the problem with increasing number of cells and compare the -# numerical and analytical solutions. To eliminate the error contribution from -# temporal discretization, we use a high number of timesteps (10 000). This -# number is arbitrarily, and must be increased if we want to study the spatial -# convergence at higher spatial resolutions. We see that the method converges as -# expected, with second order accuracy in space. +# ### Spatial convergence analysis +# We examine spatial convergence by systematically increasing the number of grid +# cells while using a very fine temporal discretization (10,000 timesteps) to +# eliminate temporal error contributions. The L2 norm of the error is computed +# over the entire space-time domain. We expect to observe second-order +# convergence as predicted by finite difference theory. Δx, _, err_space = convergence_1d(setup_case, :space) Δx ./= Δx[1] @@ -178,10 +181,11 @@ lines!(ax, Δx, err_space, linewidth=2, color=:black, label="Numerical") Legend(fig[1, 2], ax) fig -# ### Convergence in time -# Next, we study convergene in time, this time setting the number of cells to 10 -# 000 to mitigate the spatial error. Again, the method converges as expected, -# with first order accuracy in time. +# ### Temporal convergence analysis +# We now study temporal convergence by using a very fine spatial discretization +# (10,000 cells) to eliminate spatial error contributions while systematically +# reducing the timestep size. We expect to observe first-order convergence +# characteristic of the implicit Backward Euler time integration scheme. _, Δt, err_time = convergence_1d(setup_case, :time) Δt ./= Δt[1] diff --git a/examples/production/doublet.jl b/examples/production/doublet.jl index b223bbd..2b02564 100644 --- a/examples/production/doublet.jl +++ b/examples/production/doublet.jl @@ -1,16 +1,18 @@ -# # Simulation of geothermal energy production from well doublet +# # Geothermal energy production from a well doublet # This example demonstrates how to simulate and visualize the production of # geothermal energy from a doublet well system. The setup consists of a # producer well that extracts hot water from the reservoir, and an injector well -# that reinjects produced water at a lower temperature. +# that reinjects produced water at a lower temperature. + +# Add required modules to namespace using Jutul, JutulDarcy, Fimbul using HYPRE using GLMakie # ## Set up case -# The injector and producer wells are placed 100 m apart at the top, and runs -# parallel down to 800 m before they diverge to a distance of 1000 m at 2500 m -# depth. +# The injector and producer wells are placed 100 m apart at the top, and run +# parallel down to 800 m depth before they diverge to a distance of 1000 m at +# 2500 m depth. case, plot_args = geothermal_doublet(); # ## Inspect model @@ -36,7 +38,10 @@ plot_reservoir(case.model; plot_args...) # produce at a rate of 300 m^3/hour with a lower BHP limit of 1 bar, while the # injector is set to reinject the produced water, lowered to a temperature of 20 # °C. -results = simulate_reservoir(case; info_level = 0) + +# Note: this simulation can take a few minutes to run. Setting `info_level = +# 0` will show a progress bar while the simulation runs. +results = simulate_reservoir(case; info_level = -1) # ## Visualize results # We first plot the reservoir state interactively. You can notice how the @@ -44,13 +49,14 @@ results = simulate_reservoir(case; info_level = 0) plot_reservoir(case.model, results.states; colormap = :seaborn_icefire_gradient, key = :Temperature, plot_args...) -# Next, we plot the well output +# ### Plot well output +# Next, we plot the well output to examine the production rates and temperatures. plot_well_results(results.wells) # ### Well temperature evolution -# For a more detailed analysis of the preduction temperature evolution, we plot +# For a more detailed analysis of the production temperature evolution, we plot # the temperature inside the production well as a function of well depth for all -# the 200 years of production. We also highlight the temperature at selected +# 200 years of production. We also highlight the temperature at selected # timesteps (7, 21, 65, and 200 years). states = results.result.states times = convert_from_si.(cumsum(case.dt), :year) @@ -59,12 +65,13 @@ geo = tpfv_geometry(msh) fig = Figure(size = (1200, 800)) ax = Axis(fig[1, 1], xlabel = "Temperature (°C)", ylabel = "Depth (m)", yreversed = true) colors = cgrad(:seaborn_icefire_gradient, length(states), categorical = true) +# Plot temperature profiles for all timesteps with transparency for (n, state) in enumerate(states) T = convert_from_si.(state[:Producer][:Temperature], :Celsius) plot_mswell_values!(ax, case.model, :Producer, T; geo = geo, linewidth = 4, color = colors[n], alpha = 0.25) end - +# Highlight selected timesteps with solid lines and labels timesteps = [7, 21, 65, 200] for n in timesteps T = convert_from_si.(states[n][:Producer][:Temperature], :Celsius) @@ -75,7 +82,7 @@ axislegend(ax; position = :lt, fontsize = 20) fig # We can clearly see the footprint of the cold front in the aquifer (2400-2500 m -# depth) as it near the production well. +# depth) as it approaches the production well. # ### Reservoir state evolution # Another informative plot is the change in reservoir states over time. We @@ -94,14 +101,17 @@ end plot_reservoir(case.model, Δstates; colormap = :seaborn_icefire_gradient, key = :Temperature, plot_args...) -# Finally, we plot the change in temperature at the same timesteps higlighted in +# ### 3D visualization of temperature changes +# Finally, we plot the change in temperature at the same timesteps highlighted in # the production well temperature above. We cut away parts of the model for # better visualization. +# Define the temperature range and create a mask to cut away parts of the model T_min = minimum(Δstates[end][:Temperature]) T_max = maximum(Δstates[1][:Temperature]) cells = geo.cell_centroids[1, :] .> -1000.0/2 cells = cells .&& geo.cell_centroids[2, :] .> 50.0 cells = cells .|| geo.cell_centroids[3, :] .> 2475.0 +# Create subplots for each highlighted timestep fig = Figure(size = (800, 800)) for (i, n) in enumerate(timesteps) ax = Axis3(fig[(i-1)÷2+1, (i-1)%2+1]; title = "$(times[n]) years", diff --git a/examples/storage/btes.jl b/examples/storage/btes.jl index 668a7a5..b622888 100644 --- a/examples/storage/btes.jl +++ b/examples/storage/btes.jl @@ -1,6 +1,6 @@ # # Borehole Thermal Energy Storage (BTES) -# This example demonstrates how to set up and simulate a Borehole Thermal Energy -# Storage (BTES) system using Fimbul. The BTES system consists of multiple +# This example demonstrates how to set up and simulate seasonal Borehole Thermal +# Energy Storage (BTES) using Fimbul. The BTES system consists of multiple # boreholes coupled in series and parallel. using Jutul, JutulDarcy using Fimbul @@ -8,22 +8,24 @@ using HYPRE using GLMakie # ## Set up simulation case -# We consider a domain with 48 BTES wells that all reach 100 m depth. The BTES -# wells are charged during the summer at 90 °C and discharged during the winter -# months at 10 °C. This cycle is simulated for 4 years. +# We create a BTES system with 48 wells that extend to 100 m depth. The system +# operates on a seasonal cycle: wells are charged during summer months (April +# to September) with hot water at 90°C, and discharged during winter months +# (October to March) with cold water at 10°C. This seasonal operation is +# simulated over a 4-year period to study the thermal energy storage performance. case, sections = btes(num_wells = 48, depths = [0.0, 0.5, 100, 125], charge_months = ["April", "May", "June", "July", "August", "September"], discharge_months = ["October", "November", "December", "January", "February", "March"], num_years = 4, ); -# ### Visualize BTES system -# The 48 BTES wells are placed in a pattern so that all wells are approximately -# the same distance apart, and divided into 6 sections of 8 wells each. During -# charging, water is injected into the inner-most well of each section at 0.5 -# l/s, with water flowing from the inner-most well to the outer-most well. -# During discharging, water flows in the opposite direction, from the outer-most -# well to the inner-most well. +# ### Visualize BTES system layout +# The 48 BTES wells are arranged in a pattern that ensures approximately equal +# spacing between wells. The wells are organized into 6 sections of 8 wells +# each. During charging, water is injected at 0.5 l/s into the innermost well of +# each section and flows sequentially through all wells in series to the +# outermost well. During discharging, the flow direction reverses, with water +# flowing from the outermost well back to the innermost well. msh = physical_representation(reservoir_model(case.model).data_domain) fig = Figure(size = (800, 800)) ax = Axis3(fig[1, 1]; zreversed = true, aspect = :data, @@ -45,6 +47,8 @@ Legend(fig[1, 2], lns, labels); fig # ## Set up reservoir simulator +# We configure the simulator with specific tolerances and timestep controls to +# handle the challenging thermal transients in the BTES system. simulator, config = setup_reservoir_simulator(case; tol_cnv = 1e-2, tol_mb = 1e-6, @@ -52,9 +56,10 @@ simulator, config = setup_reservoir_simulator(case; initial_dt = 5.0, relaxation = true); -# Changing from charging to discharging results in a thermal shock that is -# challenging to resolve for the nonlinear solver. We therefore use a timestep -# selector that reduces the timestep to 5 seconds when the control changes. +# The transition from charging to discharging creates a thermal shock that is +# numerically challenging for the nonlinear solver. We use a specialized +# timestep selector that reduces the timestep to 5 seconds during control +# changes to maintain numerical stability and convergence. sel = JutulDarcy.ControlChangeTimestepSelector( case.model, 0.0, convert_to_si(5.0, :second)) push!(config[:timestep_selectors], sel) @@ -63,12 +68,19 @@ for ws in well_symbols(case.model) config[:tolerances][ws][:default] = 1e-2 end -# ## Simulate the case +# ## Simulate the BTES system +# Run the full 4-year simulation with the configured solver settings. + +# Note that this simulation can take a few minutes to run. Setting `info_level = +# 0` will show a progress bar while the simulation runs.` results = simulate_reservoir(case; -simulator=simulator, config=config, info_level=0); +simulator=simulator, config=config, info_level=-1); -# ### Interactive visualization -# We plot the reservoir state and the well output interactively. +# ### Interactive visualization of results +# We first plot the final temperature distribution in the reservoir and well +# performance over time. The reservoir plot shows the thermal plumes developed +# around each well section, while the well results show the temperature and flow +# rate evolution throughout the simulation. plot_reservoir(case.model, results.states; well_fontsize = 0, key = :Temperature, step = length(case.dt), colormap = :seaborn_icefire_gradient) @@ -76,9 +88,11 @@ colormap = :seaborn_icefire_gradient) plot_well_results(results.wells, field = :temperature) # ## Temperature evolution in the reservoir -# We plot the temperature evolution in the reservoir after each of the four -# charging stages. To enhance the visualization, we only plot the cells below 50 -# m depth that have a temperature above 15 °C. +# Next, we visualize how thermal energy storage develops in the subsurface over +# time. We examine the temperature distribution after each of the four charging +# cycles to see how the thermal plume grows. For better visualization, we only +# display cells below 50 m depth with temperatures above 15°C (warmer than the +# natural ground temperature). dstep = Int64.(length(case.dt)/(4*2)) steps = dstep:2*dstep:length(case.dt) fig = Figure(size = (1000, 750)) @@ -101,17 +115,18 @@ colormap = :seaborn_icefire_gradient, colorrange = (T_min, T_max), label = "T ≥ 15.0 °C", vertical = true) fig -# ## Temperature evolution in the wells -# Designing a BTES system requires in-depth understanding of how the temperature -# evolves as water runs through the pipes. We plot the temperature in the first -# section for each report step during the first charge and discharge stages. -# This kind of visualization can be used when designing a BTES system to -# determine e.g., how many wells are needed for a given depth, and how wells -# should be coupled in series/parallel. +# ## Temperature evolution within the boreholes +# Understanding the temperature profiles within individual boreholes is crucial +# for BTES system design. We analyze how temperature varies with depth along the +# wellbore during both charging and discharging phases. This information can +# help determine optimal well depth, number of wells needed, and the most +# effective series/parallel coupling configuration for a given thermal storage +# capacity requirement. -# We set up a function that plots the temperature evolution in a given section -# of the BTES system for a given set of timesteps. The function plots the -# temperature as a function of well depth for all wells in the section. +# To this end, we define a utility function to plot temperature evolution in a +# BTES section. This function visualizes temperature as a function of well depth +# for all wells in a given section at specified timesteps, helping to understand +# thermal propagation through the wellbore network. function plot_btes_temperature(ax, section, timesteps) colors = cgrad(:ice, length(timesteps), categorical = true) msh = physical_representation(reservoir_model(case.model).data_domain) @@ -134,9 +149,10 @@ end fig = Figure(size = (1000, 750)) section = 1 -time = convert_from_si.(cumsum(case.dt), :year) -# Temperature evolution during the first charging +# ### Temperature evolution during the first charging cycle +# We first visualize how thermal energy propagates through the wellbore network +# as hot water flows from the innermost to the outermost well in the section. well = sections[Symbol("S$section")][1] T_ch = convert_to_si(90.0, :Celsius) is_charge = [f[:Facility].control[well].temperature == T_ch for f in case.forces] @@ -146,7 +162,9 @@ xlabel = "T [°C]", ylabel = "Depth [m]", yreversed = true) plot_btes_temperature(ax, section, 1:stop) Legend(fig[1,2], ax; fontsize = 20) -# Temperature evolution during the first discharging +# ### Temperature evolution during the first discharging cycle +# Next, we show how stored thermal energy is extracted as flow direction +# reverses, with water flowing from the outermost to the innermost well. well = sections[Symbol("S$section")][end-1] T_dch = convert_to_si(10.0, :Celsius) is_discharge = [f[:Facility].control[well].temperature == T_dch for f in case.forces]