Refactor to new JutulDarcy release

Project.toml

@@ -1,7 +1,7 @@
 name = "Fimbul"
 uuid = "cac51d78-27de-4708-bd2b-bd13a82f652b"
 authors = ["Øystein Klemetsdal"]
-version = "0.2.1"
+version = "0.3.0"
 
 [deps]
 Dates = "ade2ca70-3891-5945-98fb-dc099432e06a"
@@ -23,7 +23,7 @@ Dates = "1.11.0"
 Gmsh = "0.3.1"
 Integrals = "4.5.0"
 Jutul = "0.4.5"
-JutulDarcy = "0.2.47"
+JutulDarcy = "0.3.0"
 LinearAlgebra = "1.11.0"
 Statistics = "1.11.1"
 julia = "1.7"

README.md

@@ -38,6 +38,24 @@ The first time you run this code, Julia will compile the packages, which may tak
 >[!NOTE]
 >Interactive plotting requires `GLMakie`, which may not work if you are running Julia over SSH.
 
+## Citing
+
+The current main work describing `Fimbul.jl` is [*Fimbul.jl – Fast, Flexible, Robust, and Differentiable Geothermal Energy Simulation in Julia*, available through EarthDoc](https://doi.org/10.3997/2214-4609.202521164):
+
+```bibtex
+@inproceedings{Klemetsdal2025,
+  title = {Fimbul.jl – Fast,  Flexible,  Robust,  and Differentiable Geothermal Energy Simulation in Julia},
+  DOI = {10.3997/2214-4609.202521164},
+  booktitle = {Sixth EAGE Global Energy Transition Conference & Exhibition (GET 2025)},
+  publisher = {European Association of Geoscientists & Engineers},
+  author = {Klemetsdal,  Ø. and Andersen,  O. and Møyner,  O.},
+  year = {2025},
+  pages = {1–5}
+}<>
+```
+
+[DOI link to extended abstract.](https://doi.org/10.3997/2214-4609.202521164) If you use Fimbul in your work, please cite this work.
+
 ## License
 
 Fimbul.jl is licensed under the MIT License. See [LICENSE](LICENSE) for details.

docs/make.jl

@@ -208,7 +208,7 @@ function build_fimbul_docs(
             modules = [Fimbul],
             authors = "Øystein Klemetsdal <oystein.klemetsdal@sintef.no> and contributors",
             repo = "https://github.com/sintefmath/Fimbul.jl/blob/{commit}{path}#{line}",
-            warnonly = true,
+            warnonly = false,
             sitename = "Fimbul.jl",
             checkdocs = :exports,
             plugins = [bib],

docs/src/index.md

@@ -53,4 +53,21 @@ The first time you run this code, Julia will compile the packages, which may tak
 >Interactive plotting requires `GLMakie`, which may not work if you are running Julia over SSH.
 
 ## Working with Fimbul
-If you plan to use Fimbul extensively in your work, we strongly recommend that you read the documentation of [JutulDarcy](https://sintefmath.github.io/JutulDarcy.jl/dev/), in particular on [getting started](https://sintefmath.github.io/JutulDarcy.jl/dev/man/intro). This also covers the basics of installing Julia and creating a Julia environments, written for users who may not already be familiar with Julia package management.
+If you plan to use Fimbul extensively in your work, we strongly recommend that you read the documentation of [JutulDarcy](https://sintefmath.github.io/JutulDarcy.jl/dev/), in particular on [getting started](https://sintefmath.github.io/JutulDarcy.jl/dev/man/intro). This also covers the basics of installing Julia and creating a Julia environments, written for users who may not already be familiar with Julia package management.
+
+## Citing
+The current main work describing `Fimbul.jl` is [*Fimbul.jl – Fast, Flexible, Robust, and Differentiable Geothermal Energy Simulation in Julia*, available through EarthDoc](https://doi.org/10.3997/2214-4609.202521164):
+
+```bibtex
+@inproceedings{Klemetsdal2025,
+  title = {Fimbul.jl – Fast,  Flexible,  Robust,  and Differentiable Geothermal Energy Simulation in Julia},
+  DOI = {10.3997/2214-4609.202521164},
+  booktitle = {Sixth EAGE Global Energy Transition Conference & Exhibition (GET 2025)},
+  publisher = {European Association of Geoscientists & Engineers},
+  author = {Klemetsdal,  Ø. and Andersen,  O. and Møyner,  O.},
+  year = {2025},
+  pages = {1–5}
+}<>
+```
+
+[DOI link to extended abstract.](https://doi.org/10.3997/2214-4609.202521164) If you use Fimbul in your work, please cite this work.

docs/src/man/cases/cases.md

@@ -12,12 +12,15 @@ egg_geothermal
 
 ## Geothermal energy production
 ```@docs
-egg_geothermal_doublet
 geothermal_doublet
+egs
+ags
+egg_geothermal_doublet
 ```
 
 ## Underground thermal energy storage
 ```@docs
-egg_ates
+ates
 btes
+egg_ates
 ```

docs/src/man/cases/utils.md

@@ -1,5 +1,6 @@
 # Utilities
 
 ```@docs
+make_schedule
 make_utes_schedule
 ```

examples/production/ags_demo.jl

@@ -0,0 +1,200 @@
+# # Advanced Geothermal System (AGS)
+# This example demonstrates simulation and analysis of energy production from an
+# Advanced Geothermal System (AGS). AGS technology utilizes closed-loop
+# circulation systems to extract geothermal energy from deep hot rock
+# formations, offering a novel approach to geothermal energy extraction that
+# doesn't require natural permeability or fracture stimulation.
+#
+# The AGS well setup for this example was provided by Alexander Rath (OMV)
+
+# Add required modules to namespace
+using Jutul, JutulDarcy, Fimbul # Core reservoir simulation framework
+using HYPRE # High-performance linear solvers
+using GLMakie # 3D visualization and plotting capabilities
+
+# Useful SI units
+meter, hour, day, watt = si_units(:meter, :hour, :day, :watt);
+
+# ## AGS setup
+# We consider an AGS system featuring a closed-loop configuration with a single
+# vertical well extending 2400 m deep, followed by two horizontal lateral
+# sections at depth for enhanced heat exchange with the surrounding rock. The
+# system includes a production well that returns heated water to the surface.
+
+# ### Create simulation case
+# We set up a scenario describing 50 years of operation with a water circulation
+# rate of 50 L/s at an injection temperature of 25°C. The simulation will
+# output results four times per year for detailed analysis.
+reports_per_year = 4 # Output frequency for results analysis
+case = Fimbul.ags(;
+    rate = 25meter^3/hour, # Water circulation rate
+    temperature_inj = convert_to_si(25.0, :Celsius), # Injection temperature
+    num_years = 50, # Years of operation
+    report_interval = si_unit(:year)/reports_per_year,
+    porosity = 0.01, # Low porosity rock matrix
+    permeability = 1e-3*si_unit(:darcy), # Low permeability formation
+    rock_thermal_conductivity = 2.5*watt/(meter*si_unit(:Kelvin)), # Rock thermal conductivity
+    rock_heat_capacity = 900.0*si_unit(:joule)/(si_unit(:kilogram)*si_unit(:Kelvin)) # Rock heat capacity
+);
+
+# ### Inspect model
+# Visualize the computational mesh and well configuration. The mesh is refined
+# around wells to accurately capture thermal and hydraulic processes in the
+# closed-loop system.
+msh = physical_representation(reservoir_model(case.model).data_domain)
+geo = tpfv_geometry(msh)
+fig = Figure(size = (800, 800))
+ax = Axis3(fig[1, 1]; zreversed = true, aspect = :data, perspectiveness = 0.5,
+    title = "AGS System: Closed-loop well and mesh")
+Jutul.plot_mesh_edges!( # Show computational mesh with transparency
+    ax, msh, alpha = 0.2)
+wells = get_model_wells(case.model)
+function plot_ags_wells( # Utility to plot wells in AGS system
+    ax; colors = [:black, :black])
+    for (i, (name, well)) in enumerate(wells)
+        color = colors[i]
+        label = string(name)
+        if haskey(well.perforations, :reservoir)
+            cells = well.perforations.reservoir
+            if length(cells) > 0
+                xy = geo.cell_centroids[1:2, cells]
+                xy = hcat(xy[:,1], xy)'
+                plot_mswell_values!(ax, case.model, name, xy;
+                    geo = geo, linewidth = 3, color = color, label = label)
+            end
+        end
+    end
+end
+plot_ags_wells(ax)
+fig
+
+# ## Simulate system
+# We set up the the simulator
+sim, cfg = setup_reservoir_simulator(case;
+    output_substates = true, # Store results from timesteps between
+    info_level = 0, # 0=progress bar, 1=basic, 2=detailed
+    initial_dt = 5.0, # Initial timestep [s]
+    presolve_wells = true, # Solve wells with fixed reservoir state at the beginning of each timestep
+    relaxation = true); # Enable relaxation in Newton solver
+
+# We add a specialized timestep selector to control solution quality during
+# thermal transients. This selector monitors temperature changes and adjusts
+# timesteps aiming at a maximum change of 5°C per timestep in both the reservoir
+# and the well.
+sel = VariableChangeTimestepSelector(:Temperature, 5.0; 
+    relative = false, model = :Reservoir)
+push!(cfg[:timestep_selectors], sel);
+sel = VariableChangeTimestepSelector(:Temperature, 5.0; 
+    relative = false, model = :AGS_supply)
+push!(cfg[:timestep_selectors], sel);
+
+# NOTE: depending in your system, the simulation may take a few minutes to
+# complete.
+results = simulate_reservoir(case; simulator = sim, config = cfg)
+
+# ## Interactive Visualization
+# Next, we analyze and visualize the simulation results interactively to
+# understand the AGS performance, thermal depletion patterns, and energy
+# production characteristics throughout the 50-year operational period.
+
+# ### Reservoir state evolution
+# It is often most informative to visualize the deviation from the initial
+# conditions to highlight the  extent of the thermal depletion zones around the
+# AGS system. We compute the change in reservoir variables to the initial state
+# for all timesteps.
+Δstates = JutulDarcy.delta_state(results.states, case.state0[:Reservoir])
+plot_res_args = (
+    resolution = (600, 800), aspect = :data, 
+    colormap = :seaborn_icefire_gradient, key = :Temperature,
+    well_arg = (markersize = 0.0, ),
+)
+plot_reservoir(case.model, Δstates; plot_res_args...)
+
+# ### Final temperature change in the reservoir
+# We visualize the final temperature change in the reservoir after 50 years of
+# operation, with a seubset of cells cut out for better visibility.
+
+# Define cells to cut out
+cut_out = geo.cell_centroids[1, :] .< 750.0
+cut_out = cut_out .|| geo.cell_centroids[2, :] .< 0.0
+cut_out = cut_out .&& geo.cell_centroids[3, :] .< 2400
+
+fig = Figure(size = (800, 800))
+tot_time = round(sum(case.dt)/si_unit(:year), digits = 1)
+ax = Axis3(fig[1, 1]; title = "Temperature after $(tot_time) years",
+zreversed = true, elevation = pi/8, 
+aspect = :data, perspectiveness = 0.5)
+plt = plot_cell_data!(ax, msh, Δstates[end][:Temperature]; 
+    cells = .!cut_out, colormap = :seaborn_icefire_gradient)
+[plot_well!(ax, msh, well; 
+    color = :black, markersize = 0.0, fontsize = 0.0, linewidth = 1) 
+    for well in values(wells)]
+
+Colorbar(fig[2,1], plt;
+    label = "ΔT (°C)", vertical = false)
+fig
+
+# ## Well Performance Analysis
+# Examine the well responses including flow rates, pressures, and temperatures.
+# The AGS system shows the circulation flow through the closed loop and the
+# thermal response as the system extracts heat from the surrounding rock matrix.
+plot_well_results(results.wells)
+
+# ## Lateral Section Analysis
+# We analyze the performance of the lateral sections in the AGS system by
+# extracting temperature and power data along the lateral segments over time. This
+# analysis helps to understand how effectively the two laterals exchange heat with
+# the reservoir and contribute to overall energy production.
+section_data = Fimbul.get_section_data_ags(
+    case, results.result.states, :AGS_supply)
+
+# ### Set up plotting utilities
+colors = collect(cgrad(:BrBg, 8, categorical = true))[[2, end-1]]
+function plot_lateral_data!(ax, time, data; stacked = false)
+    num_laterals = size(data, 2)-2
+    y_prev = zeros(size(data, 1))
+    for lno = 1:num_laterals
+        y = data[:, lno+1]
+        if stacked
+            x = vcat(time, reverse(time))
+            println("Size x: ", size(x), ", size y: ", size(y))
+            y .+= y_prev
+            y = vcat(y_prev, reverse(y))
+            poly!(ax, x, y; color = colors[lno],
+            strokecolor = :black, strokewidth = 1, label = "Lateral $lno")
+            y_prev = data[:, lno+1]
+        else
+            lines!(ax, time, y; color = colors[lno], linewidth = 4, label = "Lateral $lno")
+        end
+    end
+end
+
+# Plot lateral temperature and power over time
+fig = Figure(size = (800, 800))
+time = results.time ./ si_unit(:year)
+
+ax_tmp = Axis( # Panel 1: Lateral temperature
+    fig[1, 1:3]; title = "Lateral Temperature", 
+    ylabel = "Temperature (°C)", xlabel = "Time (years)")
+T = convert_from_si.(section_data[:Temperature], :Celsius)
+plot_lateral_data!(ax_tmp, time, T, stacked = false)
+hidexdecorations!(ax_tmp, grid = false)
+
+ax_pwr = Axis( # Panel 2: Lateral power
+    fig[2, 1:3]; title = "Lateral Power", 
+    ylabel = "Power (W)", xlabel = "Time (years)",
+limits = (nothing, (-0.05, 0.5)))
+MW = si_unit(:mega)*si_unit(:watt)
+plot_lateral_data!(ax_pwr, time, section_data[:Power]./MW, stacked = true)
+
+ax_lat = Axis( # Panel 3: Lateral well trajectories for reference
+    fig[1:2, 4]; aspect = DataAspect(), limits = ((-150, 100), nothing), 
+    xticks = [-150, 100])
+well_coords, _ = Fimbul.get_ags_trajectory()
+for (k, wc) in enumerate(well_coords[2:3])
+    lines!(ax_lat, wc[:, 2], wc[:, 1]; linewidth = 4, color = colors[k])
+end
+
+axislegend( # Add legend
+    ax_tmp; position = :rt, fontsize = 20)
+fig

examples/production/doublet_demo.jl

@@ -13,15 +13,15 @@ using GLMakie
 # 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();
+case = geothermal_doublet();
 
 # ## Inspect model
 # We first plot the computational mesh and wells. The mesh is refined around
 # the wells in the horizontal plane and vertically in and near the target
 # aquifer.
 msh = physical_representation(reservoir_model(case.model).data_domain)
 fig = Figure(size = (1200, 800))
-ax = Axis3(fig[1, 1], zreversed = true, aspect = plot_args.aspect)
+ax = Axis3(fig[1, 1], zreversed = true, aspect = :data)
 Jutul.plot_mesh_edges!(ax, msh, alpha = 1.0)
 wells = get_model_wells(case.model)
 for (k, w) in wells
@@ -31,7 +31,7 @@ fig
 
 # ### Plot reservoir properties
 # Next, we visualize the reservoir interactively.
-plot_reservoir(case.model; plot_args...)
+plot_reservoir(case.model; aspect = :data)
 
 # ## Simulate geothermal energy production
 # We simulate the geothermal doublet for 200 years. The producer is set to
@@ -47,7 +47,7 @@ results = simulate_reservoir(case; info_level = 0)
 # 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...)
+colormap = :seaborn_icefire_gradient, key = :Temperature, aspect = :data)
 
 # ### Plot well output
 # Next, we plot the well output to examine the production rates and temperatures.
@@ -72,7 +72,7 @@ for (n, state) in enumerate(states)
     geo = geo, linewidth = 4, color = colors[n], alpha = 0.25)
 end
 # Highlight selected timesteps with solid lines and labels
-timesteps = [7, 21, 65, 200]
+timesteps = [12, 25, 50, 100]
 for n in timesteps
     T = convert_from_si.(states[n][:Producer][:Temperature], :Celsius)
     plot_mswell_values!(ax, case.model, :Producer, T;
@@ -99,7 +99,7 @@ for state in results.states
     push!(Δstates, Δstate)
 end
 plot_reservoir(case.model, Δstates;
-colormap = :seaborn_icefire_gradient, key = :Temperature, plot_args...)
+colormap = :seaborn_icefire_gradient, key = :Temperature, aspect = :data)
 
 # ### 3D visualization of temperature changes
 # Finally, we plot the change in temperature at the same timesteps highlighted in

examples/storage/ates_demo.jl

@@ -40,7 +40,7 @@ darcy = si_unit(:darcy);
 # realistic geological and operational parameters for a medium-scale ATES
 # installation with 400m well spacing.
 num_years = 5
-case, layers = Fimbul.ates(;
+case = Fimbul.ates(;
     well_distance = 400.0, # Distance between wells [m]
     temperature_charge = convert_to_si(85, :Celsius),   # Hot injection temperature
     temperature_discharge = convert_to_si(20, :Celsius), # Cold injection temperature
@@ -67,7 +67,7 @@ colors = [:red, :blue]  # Hot well = red, Cold well = blue
 for (i, (k, w)) in enumerate(wells)
     plot_well!(ax, msh, w, color = colors[i], linewidth = 6)
 end
-plot_cell_data!(ax, msh, layers,
+plot_cell_data!(ax, msh, case.input_data[:layers],
     colormap = :rainbow,
     alpha = 0.3)
 fig
@@ -247,6 +247,7 @@ fig_wells
 
 # Extract temperature along a horizontal line in the aquifer layer
 # This transect passes through the center of the aquifer between the two wells
+layers = case.input_data[:layers]
 ijk = [cell_ijk(msh, c) for c in 1:number_of_cells(msh)]
 j = div(maximum(getindex.(ijk, 2)) + minimum(getindex.(ijk, 2)), 2)
 k = div(maximum(getindex.(ijk[layers.==3], 3)) + minimum(getindex.(ijk[layers.==3], 3)), 2)
@@ -270,7 +271,10 @@ function plot_aquifer_temperature!(fig, T_line, stage, cycle)
     else
         steps = dch_start[cycle]:dch_stop[cycle]
     end
-    colors = cgrad(:seaborn_icefire_gradient, length(steps), categorical = true)
+    colors = cgrad(:RdBu, length(steps), categorical = true)
+    if stage == "Charging"
+        colors = reverse(colors)
+    end
     for (n, T_n) in enumerate(T_line[steps])
         T_n = convert_from_si.(T_n, :Celsius)
         lines!(ax, x, T_n, color = colors[n], linewidth = 3,