Geometry updates

examples/production/egs_demo.jl

@@ -30,7 +30,7 @@ meter, day, watt = si_units(:meter, :day, :watt);
 well_depth = 2500.0meter # Vertical depth to horizontal well section [m]
 well_spacing_x = 100.0meter # Horizontal separation between production wells [m]
 well_spacing_z = sqrt(3/4*well_spacing_x^2) # Vertical offset between injector and producers [m]
-well_lateral = 500.0meter; # Length of horizontal well section [m]
+well_lateral = 1000.0meter; # Length of horizontal well section [m]
 
 # Well coordinates are defined as a vector of nx3 arrays, each containing the
 # (x,y,z) coordinates of a well trajectory. The first well is the injector,
@@ -52,12 +52,12 @@ fracture_spacing = well_lateral/8; # Spacing between discrete fractures [m]
 # We set up a scenario describing 10 years of operation with a water injection
 # rate of 9250 m³/day (approximately 107 liters/second) at a temperature of
 # 25°C. The simulation will output results four times per year for analysis.
-reports_per_year = 4 # Output frequency for results analysis
+num_years = 10 # Total simulation period [years]
 case = Fimbul.egs(well_coords, fracture_radius, fracture_spacing;
     rate = 9250meter^3/day, # Water injection rate
     temperature_inj = convert_to_si(25.0, :Celsius), # Injection temperature
-    num_years = 10, # Years of operation
-    schedule_args = (report_interval = si_unit(:year)/reports_per_year,)
+    num_years = num_years, # Years of operation
+    schedule_args = (report_interval = si_unit(:year)/4,)
 );
 
 # ### Inspect model
@@ -101,12 +101,11 @@ fig
 # we use increment-based criteria that directly monitor changes in primary
 # variables (temperature and pressure) between Newton iterations.
 sim, cfg = setup_reservoir_simulator(case;
-    info_level = 0, # 0=progress bar, 1=basic, 2=detailed
-    tol_cnv = Inf, # Disable CNV tolerance - use other criteria
-    inc_tol_dT = 1e-2, # Temperature increment tolerance [K]
-    inc_tol_dp_rel = 1e-3, # Relative pressure increment tolerance [-]
+    info_level = 2, # 0=progress bar, 1=basic, 2=detailed
+    output_substates = true, # Output results at each Newton iteration
     initial_dt = 5.0, # Initial timestep [s]
-    relaxation = true); # Enable relaxation in Newton solver
+    relaxation = true # Enable relaxation in Newton solver
+);
 
 # We add a specialized timestep selector to control solution quality during
 # thermal transients. These selectors monitor temperature changes and adjust
@@ -153,8 +152,7 @@ plot_well_results(results.wells)
 # network at selected timesteps to understand thermal depletion patterns.
 
 # Extract fracture temperature changes and prepare plotting
-dt = case.dt
-time = convert_from_si.(cumsum(dt), :year)
+time = convert_from_si.(results.time, :year)
 ΔT = [Δstate[:Temperature][is_fracture] for Δstate in Δstates]
 colorrange = extrema(vcat(ΔT...))
 x = geo.cell_centroids[:, is_fracture]
@@ -163,7 +161,9 @@ limits = Tuple([xl .+ (xl[2]-xl[1]).*(-0.1, 0.1) for xl in xlim])
 
 # Visualize fracture temperature at three representative timesteps
 fig = Figure(size = (650, 800))
-steps = Int.(round.(range(1, length(ΔT), 3)))
+n_steps = length(ΔT)
+steps = Int.(round.([0.125, 0.25, 1.0].* n_steps))
+
 for (n, ΔT_n) in enumerate(ΔT[steps])
     ax_n = Axis3(fig[n, 1],
     perspectiveness = 0.5,
@@ -191,20 +191,22 @@ fig
 # ### Fracture metrics
 # Finally, we analyze key performance metrics for each individual fracture over
 # the entire simulation period, including temperature evolution, thermal power
-# production, annual energy production. Notice how the first three fractures
-# contribute more than 50 % of the total energy production throughout the
-# simulation period. The last three fractures contribute only 21 % in the first
-# year, but their relative contribution increases to 27 % in the final year.
+# production, annual energy production. The temperature is almost identical out
+# of each fracture. However, we notice the first and last fractures contribute
+# slightly more to the total energy production.
 
 # Extract fracture data from simulation results
-states = results.result.states;
-fdata_inj = Fimbul.get_egs_fracture_data(states, case.model, :Injector; geo=geo);
-fdata_prod = Fimbul.get_egs_fracture_data(states, case.model, :Producer; geo=geo);
+states, dt, report_step = Jutul.expand_to_ministeps(results.result)
+time = cumsum(dt)./si_unit(:year)
+
+fdata_inj = Fimbul.get_egs_fracture_data(states, dt,case.model, :Injector; geo=geo);
+fdata_prod = Fimbul.get_egs_fracture_data(states, dt, case.model, :Producer; geo=geo);
 nsteps = length(dt)
 
-colors = reverse(cgrad(:seaborn_icefire_gradient, size(fdata_inj[:Temperature], 2), categorical = true))
+colors = cgrad(:BrBg, size(fdata_inj[:Temperature], 2), categorical = true)
 function plot_fracture_data( # Utility for plotting fracture data
-    ax, data, stacked = false)
+    ax, time, data, stacked = false)
+    nsteps = length(time)
     df_prev = zeros(nsteps)
     for (fno, df) in enumerate(eachcol(data))
         if stacked
@@ -233,26 +235,33 @@ function make_axis( # Utility for making axes
     return ax
 end
 
+first_step = findfirst(time .> 1/104)
+
 ax = make_axis( # Panel 1: fracture temperature
     "Temperature", "T (°C)", 1)
-temperature = fdata_prod[:Temperature]
-plot_fracture_data(ax, convert_from_si.(temperature, :Celsius), false)
+temperature = fdata_prod[:Temperature][first_step:end, :]
+plot_fracture_data(ax, time[first_step:end], convert_from_si.(temperature, :Celsius), false)
 hidexdecorations!(ax, grid = false)
 
 ax_pwr = make_axis( # Panel 2: thermal power production
     "Thermal power", "Power (MW)", 2)
-effect = .-(fdata_prod[:EnergyFlux] .+ fdata_inj[:EnergyFlux])
-plot_fracture_data(ax_pwr, effect./1e6, true)
+
+Er = fdata_prod[:TotalThermalEnergy]
+ΔErΔt⁻¹ = diff(vcat(Er[1,:]', Er), dims=1)./dt
+power = .-(fdata_prod[:EnergyFlux] .+ fdata_inj[:EnergyFlux]) .+ ΔErΔt⁻¹
+
+# power = fdata_prod[:EnergyFlux] # Power production from producers only --- IGNORE ---
+plot_fracture_data(ax_pwr, time[first_step:end], power[first_step:end, :]./1e6, true)
 hidexdecorations!(ax_pwr, grid = false)
 
 ax = make_axis( # Panel 3: annual energy production per fracture
     "Annual energy production per fracture", "Energy (GWh)", 3)
 energy_per_year, cat, dodge = [], Int[], Int[]
-n = reports_per_year
 GWh = si_unit(:giga)*si_unit(:watt)*si_unit(:hour)
-for (fno, eff_f) in enumerate(eachcol(effect))
+ix = vcat(0, [findfirst(isapprox.(time, y; atol=1e-2)) for y in 1:num_years]) .+1
+for (fno, pwr_f) in enumerate(eachcol(power))
     ## Integrate power over annual periods to get energy [GWh]
-    energy_f = [sum(eff_f[k:k+n-1].*dt[k:k+n-1])./GWh for k = 1:n:length(dt)-n+1]
+    energy_f = [sum(pwr_f[ix[k]:ix[k+1]-1].*dt[ix[k]:ix[k+1]-1])./GWh for k = 1:length(ix)-1]
     push!(energy_per_year, energy_f)
     push!(cat, 1:length(energy_f)...)
     push!(dodge, fill(fno, length(energy_f))...)
@@ -270,4 +279,4 @@ dodge = dodge, color = colors[dodge], strokecolor = :black, strokewidth = 1)
 
 Legend( # Add legend indicating fracture numbers
     fig[2:3,2], ax_pwr)
-fig
+fig

src/cases/ags.jl

@@ -152,14 +152,25 @@ function get_ags_constraints(well_coords; hxy_min)
         push!(well_coords_2x, wc_right)
     end
 
-    cell_constraints = []
+    cell_constraints = Vector{Matrix{Float64}}()
     for wc in well_coords_2x
-        xw = [Tuple(x) for x in eachrow(unique(wc[:, 1:2], dims=1))]
-        for cc in cell_constraints
-            cor = [norm(collect(x) .- collect(y), 2) for x in xw, y in cc]
-            xw = [xw[i] for i in eachindex(xw) if all(cor[i, :] .>= Δ)]
+        # Get unique coordinates as 2×N matrix
+        cc_new = unique(wc[:, 1:2], dims=1)'  # Transpose to get 2×N
+        # Filter based on distances to existing constraints
+        if !isempty(cell_constraints)
+            # Convert matrix to vector of points for distance comparison
+            for cc in cell_constraints
+                remove = falses(size(cc_new, 2))
+                for (k, x) in enumerate(eachcol(cc_new))
+                    if any(norm(x - y, 2) < Δ for y in eachcol(cc))
+                        remove[k] = true
+                    end
+                end
+                cc_new = cc_new[:, .!remove]
+            end
         end
-        push!(cell_constraints, xw)
+
+        push!(cell_constraints, cc_new)
     end
 
     return cell_constraints

src/cases/btes.jl

@@ -55,7 +55,7 @@ function btes(;
 
     # ## Create mesh
     x = fibonacci_pattern_2d(num_wells; spacing = well_spacing)
-    well_coordinates = map(x -> [x], x)
+    well_coordinates = [x[:, i:i] for i in 1:size(x, 2)]  # Convert 2×N to vector of 2×1 matrices
     hz = diff(depths)./n_z
     hxy = well_spacing/n_xy
 
@@ -82,7 +82,7 @@ function btes(;
     for (wno, xw) in enumerate(well_coordinates)
         name = Symbol("B$wno")
         println("Adding well $name ($wno/$num_wells)")
-        xw = xw[1]
+        xw = xw[:, 1]  # Extract 2D coordinates from 2×1 matrix
         d = max(norm(xw, 2))
         v = (d > 0) ? xw./d : (1.0, 0.0)
         # Shift coordiates a bit to avoid being exactly on the node

src/cases/doublet.jl

@@ -56,9 +56,9 @@ function geothermal_doublet(;
 
     well_coords = [trajectory_inj, trajectory_prod]
 
-    xw_inj = [Tuple(x) for x in eachrow(unique(trajectory_inj[:, 1:2], dims=1))]
-    xw_prod = [Tuple(x) for x in eachrow(unique(trajectory_prod[:, 1:2], dims=1))]
-    xw = vcat([xw_inj], [xw_prod])
+    xw_inj = unique(trajectory_inj[:, 1:2], dims=1)'  # Convert to 2×N matrix
+    xw_prod = unique(trajectory_prod[:, 1:2], dims=1)'  # Convert to 2×N matrix 
+    xw = [xw_inj, xw_prod]
 
     # Create the mesh
     thickness = diff(depths)