Improve example text

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")

docs/src/extras/refs.md

@@ -0,0 +1,2 @@
+```@bibliography
+```

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]

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,21 +38,25 @@ 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
 # cold front propagates from the injector well by filtering out high values.
 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",