Merge pull request #71 from sintefmath/dev

Upscaling/coarsening of models, CGNet example, relperm example, LET functions, new GUIs, rewritten well plotter

.gitignore

@@ -31,3 +31,4 @@ Manifest.toml
 /docs/node_modules
 docs/package-lock.json
 docs/package.json
+*.dat

Project.toml

@@ -1,7 +1,7 @@
 name = "JutulDarcy"
 uuid = "82210473-ab04-4dce-b31b-11573c4f8e0a"
 authors = ["Olav Møyner <[email protected]>"]
-version = "0.2.34"
+version = "0.2.35"
 
 [deps]
 AlgebraicMultigrid = "2169fc97-5a83-5252-b627-83903c6c433c"
@@ -53,10 +53,10 @@ Dates = "1"
 DelimitedFiles = "1.6"
 DocStringExtensions = "0.9.3"
 ForwardDiff = "0.10.30"
-GLMakie = "0.10.0"
+GLMakie = "0.10.13"
 GeoEnergyIO = "1.1.12"
-HYPRE = "1.6.0"
-Jutul = "0.2.39"
+HYPRE = "1.6.0, 1.7"
+Jutul = "0.2.40"
 Krylov = "0.9.1"
 LazyArtifacts = "1"
 LinearAlgebra = "1"

docs/Project.toml

@@ -9,6 +9,7 @@ GraphMakie = "1ecd5474-83a3-4783-bb4f-06765db800d2"
 HYPRE = "b5ffcf37-a2bd-41ab-a3da-4bd9bc8ad771"
 Jutul = "2b460a1a-8a2b-45b2-b125-b5c536396eb9"
 JutulDarcy = "82210473-ab04-4dce-b31b-11573c4f8e0a"
+LBFGSB = "5be7bae1-8223-5378-bac3-9e7378a2f6e6"
 LayeredLayouts = "f4a74d36-062a-4d48-97cd-1356bad1de4e"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 Literate = "98b081ad-f1c9-55d3-8b20-4c87d4299306"

docs/make.jl

@@ -33,6 +33,9 @@ function build_jutul_darcy_docs(build_format = nothing;
         "Quarter-five-spot with variation" => "five_spot_ensemble",
         "Gravity circulation with CPR preconditioner" => "two_phase_unstable_gravity",
         "CO2 property calculations" => "co2_props",
+        "Relative permeability" => "relperms",
+        "Model coarsening" => "model_coarsening",
+        "History matching a coarse model - CGNet" => "cgnet_egg",
         "CO2 injection in saline aquifer" => "co2_sloped",
         "Compositional with five components" => "compositional_5components",
         "Parameter matching of Buckley-Leverett" => "optimize_simple_bl",

docs/src/man/highlevel.md

@@ -22,33 +22,44 @@ Jutul.UnstructuredMesh
 
 ### Reservoir
 
+Once a mesh has been set up, we can turn it into a reservoir with static properties:
+
 ```@docs
 reservoir_domain
 ```
 
 ### Wells
 
+Wells are most easily created using utilities that act directly on a reservoir domain:
+
 ```@docs
 setup_well
 setup_vertical_well
 ```
 
 ### Model
 
+A single, option-heavy function is used to set up the reservoir model and default parameters:
+
 ```@docs
 setup_reservoir_model
 ```
 
 ### Initial state
 
+The initial state can be set up by explicitly setting all primary variables. JutulDarcy also contains functionality for initial hydrostatic equilibriation of the state, but this is at the moment most easily set up using input files with the `EQUIL` keyword.
+
 ```@docs
 setup_reservoir_state
+JutulDarcy.equilibriate_state
 ```
 
-TODO: Write about hydrostatic equilbriation.
-
 ## Simulation
 
+Simulating is done by either setting up a reservoir simulator and then simulating, or by using the convenience function that automatically sets up a simulator for you.
+
+There are a number of different options available to tweak the tolerances, timestepping and behavior of the simulation. It is advised to read through the documentation in detail before running very large simulations.
+
 ```@docs
 simulate_reservoir
 setup_reservoir_simulator

docs/src/refs.bib

@@ -125,3 +125,14 @@ @article{salo_co2
   pages={34--48},
   year={2024}
 }
+
+@article{cgnet1,
+  title={Data-driven modelling with coarse-grid network models},
+  author={Lie, Knut-Andreas and Krogstad, Stein},
+  journal={Computational Geosciences},
+  volume={28},
+  number={2},
+  pages={273--287},
+  year={2024},
+  publisher={Springer}
+}

examples/cgnet_egg.jl

@@ -0,0 +1,252 @@
+# # History matching a coarse model - CGNet
+# This example demonstrates how to calibrate a coarse model against results from
+# a fine model. We do this by optimizing the parameters of the coarse model to
+# match the well curves. This is a implementation of the method described in
+# [cgnet1](@cite). This also serves as a demonstration of how to use the
+# simulator for history matching, as the fine model results can stand in for
+# real field observations.
+#
+# ## Load and simulate Egg base case
+using Jutul, JutulDarcy, HYPRE, GeoEnergyIO, GLMakie
+import LBFGSB as lb
+
+egg_dir = JutulDarcy.GeoEnergyIO.test_input_file_path("EGG")
+data_pth = joinpath(egg_dir, "EGG.DATA")
+
+fine_case = setup_case_from_data_file(data_pth)
+simulated_fine = simulate_reservoir(fine_case)
+plot_reservoir(fine_case, simulated_fine.states, key = :Saturations, step = 100)
+# ## Create initial coarse model and simulate
+coarse_case = JutulDarcy.coarsen_reservoir_case(fine_case, (25, 25, 5), method = :ijk)
+simulated_coarse = simulate_reservoir(coarse_case, info_level = -1)
+plot_reservoir(coarse_case, simulated_coarse.states, key = :Saturations, step = 100)
+# ## Setup optimization
+# We set up the optimization problem by defining the objective function as a sum
+# of squared mismatches for all well observations, for all time-steps. We also
+# define limits for the parameters, and set up the optimization problem.
+function setup_optimization_cgnet(case_c, case_f, result_f)
+    states_f = result_f.result.states
+    wells_results, = result_f
+    model_c = case_c.model
+    state0_c = setup_state(model_c, case_c.state0);
+    param_c = setup_parameters(model_c)
+    forces_c = case_c.forces
+    dt = case_c.dt
+    model_f = case_f.model
+
+    bhp = JutulDarcy.BottomHolePressureTarget(1.0)
+    wells = collect(keys(JutulDarcy.get_model_wells(case_f)))
+
+    day = si_unit(:day)
+    wrat_scale = (1/150)*day
+    orat_scale = (1/80)*day
+    grat_scale = (1/1000)*day
+
+    w = []
+    matches = []
+    signs = []
+    sys = reservoir_model(model_f).system
+    wrat = SurfaceWaterRateTarget(-1.0)
+    orat = SurfaceOilRateTarget(-1.0)
+    grat = SurfaceGasRateTarget(-1.0)
+
+    push!(matches, bhp)
+    push!(w, 1.0/si_unit(:bar))
+    push!(signs, 1)
+
+    for phase in JutulDarcy.get_phases(sys)
+        if phase == LiquidPhase()
+            push!(matches, orat)
+            push!(w, orat_scale)
+            push!(signs, -1)
+
+        elseif phase == VaporPhase()
+            push!(matches, grat)
+            push!(w, grat_scale)
+            push!(signs, -1)
+        else
+            @assert phase == AqueousPhase()
+            push!(matches, wrat)
+            push!(w, wrat_scale)
+            push!(signs, -1)
+        end
+    end
+    signs = zeros(Int, length(signs))
+    o_scale = 1.0/(sum(dt)*length(wells))
+    G = (model_c, state_c, dt, step_no, forces) -> well_mismatch(
+        matches,
+        wells,
+        model_f,
+        states_f,
+        model_c,
+        state_c,
+        dt,
+        step_no,
+        forces,
+        weights = w,
+        scale = o_scale,
+        signs = signs
+    )
+
+    @assert Jutul.evaluate_objective(G, model_f, states_f, dt, case_f.forces) == 0.0
+    ##
+    cfg = optimization_config(model_c, param_c,
+        use_scaling = true,
+        rel_min = 0.001,
+        rel_max = 1000
+    )
+    for (k, v) in cfg
+        for (ki, vi) in v
+            if ki == :FluidVolume
+                vi[:active] = k == :Reservoir
+            end
+            if ki == :ConnateWater
+                vi[:active] = false
+            end
+            if ki in [:TwoPointGravityDifference, :PhaseViscosities, :PerforationGravityDifference]
+                vi[:active] = false
+            end
+            if ki in [:WellIndices, :Transmissibilities]
+                vi[:active] = true
+                vi[:abs_min] = 0.0
+                vi[:abs_max] = 1e-6
+            end
+
+        end
+    end
+    opt_setup = setup_parameter_optimization(model_c, state0_c, param_c, dt, forces_c, G, cfg);
+    x0 = opt_setup.x0
+    F0 = opt_setup.F!(x0)
+    dF0 = opt_setup.dF!(similar(x0), x0)
+    println("Initial objective: $F0, gradient norm $(sum(abs, dF0))")
+    return opt_setup
+end
+
+# ## Define the optimization loop
+# JutulDarcy can use any optimization package that can work with gradients and
+# limits, here we use the `LBFGSB` package.
+function optimize_cgnet(opt_setup)
+    lower = opt_setup.limits.min
+    upper = opt_setup.limits.max
+    x0 = opt_setup.x0
+    n = length(x0)
+    setup = Dict(:lower => lower, :upper => upper, :x0 => copy(x0))
+
+    prt = 1
+    f! = (x) -> opt_setup.F_and_dF!(NaN, nothing, x)
+    g! = (dFdx, x) -> opt_setup.F_and_dF!(NaN, dFdx, x)
+    results, final_x = lb.lbfgsb(f!, g!, x0, lb=lower, ub=upper,
+        iprint = prt,
+        factr = 1e12,
+        maxfun = 50,
+        maxiter = 50,
+        m = 20
+    )
+    return (final_x, results, setup)
+end
+# ## Run the optimization
+opt_setup = setup_optimization_cgnet(coarse_case, fine_case, simulated_fine);
+final_x, results, setup = optimize_cgnet(opt_setup);
+# ### Transfer the results back
+# The optimization is generic and works on a single long vector that represents
+# all our active parameters. We can devectorize this vector back into the nested
+# representation used by the model itself and simulate.
+tuned_case = deepcopy(opt_setup.data[:case])
+model_c = coarse_case.model
+model_f = fine_case.model
+param_c = tuned_case.parameters
+data = opt_setup.data
+devectorize_variables!(param_c, model_c, final_x, data[:mapper], config = data[:config])
+
+
+simulated_tuned = simulate_reservoir(tuned_case, info_level = -1);
+# ### Plot the results interactively
+using GLMakie
+
+wells_f, = simulated_fine
+wells_c, = simulated_coarse
+wells_t, states_t, time = simulated_tuned
+
+plot_well_results([wells_f, wells_c, wells_t], time, names = ["Fine", "CGNet-initial", "CGNet-tuned"])
+# ### Create a function to compare individual wells
+# We next compare individual wells to see how the optimization has affected the
+# match between the coarse scale and fine scale. As we can see, we have
+# reasonably good match between the original model with about 18,000 cells and
+# the coarse model with about 3000 cells. Even better match could be possible by
+# adding more coarse blocks, or also optimizing for example the relative
+# permeability parameters for the coarse model.
+#
+# We plot the water cut and total rate for the production wells, and the bottom
+# hole pressure for the injection wells.
+function plot_tuned_well(k, kw; lposition = :lt)
+    fig = Figure()
+    ax = Axis(fig[1, 1], title = "$k", xlabel = "days", ylabel = "$kw")
+    t = wells_f.time./si_unit(:day)
+    if kw == :wcut
+        F = x -> x[k, :wrat]./x[k, :lrat]
+    else
+        F = x -> abs.(x[k, kw])
+    end
+
+    lines!(ax, t, F(wells_f), label = "Fine-scale")
+    lines!(ax, t, F(wells_c), label = "Initial coarse")
+    lines!(ax, t, F(wells_t), label = "Tuned coarse")
+    axislegend(position = lposition)
+    fig
+end
+# ### Plot PROD1 water cut
+plot_tuned_well(:PROD1, :wcut)
+# ### Plot PROD2 water cut
+plot_tuned_well(:PROD2, :wcut)
+# ### Plot PROD4 water cut
+plot_tuned_well(:PROD4, :wcut)
+# ### Plot PROD1 total rate
+plot_tuned_well(:PROD1, :rate)
+# ### Plot PROD2 total rate
+plot_tuned_well(:PROD2, :rate)
+# ### Plot PROD4 total rate
+plot_tuned_well(:PROD4, :rate)
+# ### Plot INJECT1 bhp
+plot_tuned_well(:INJECT1, :bhp, lposition = :rt)
+# ### Plot INJECT4 bhp
+plot_tuned_well(:INJECT4, :bhp, lposition = :rt)
+
+# ## Plot the objective evaluations during optimization
+fig = Figure()
+ys = log10
+is = x -> x
+ax1 = Axis(fig[1, 1], yscale = ys, title = "Objective evaluations", xlabel = "Iterations", ylabel = "Objective")
+plot!(ax1, opt_setup[:data][:obj_hist][2:end])
+fig
+# ## Plot the evoluation of scaled parameters
+# We show the difference between the initial and final values of the scaled
+# parameters, as well as the lower bound.
+#
+# JutulDarcy maps the parameters to a single vector for optimization with values
+# that are approximately in the box limit range (0, 1). This is convenient for
+# optimizers, but can also be useful when plotting the parameters, even if the
+# units are not preserved in this visualization, only the magnitude.
+fig = Figure(size = (800, 600))
+ax1 = Axis(fig[1, 1], title = "Scaled parameters", ylabel = "Scaled value")
+scatter!(ax1, setup[:x0], label = "Initial X")
+scatter!(ax1, final_x, label = "Final X", markersize = 5)
+lines!(ax1, setup[:lower], label = "Lower bound")
+axislegend()
+
+trans = data[:mapper][:Reservoir][:Transmissibilities]
+
+function plot_bracket(v, k)
+    start = v.offset+1
+    stop = v.offset+v.n
+    y0 = setup[:lower][start]
+    y1 = setup[:lower][stop]
+    bracket!(ax1, start, y0, stop, y1,
+    text = "$k", offset = 1, orientation = :down)
+end
+
+for (k, v) in pairs(data[:mapper][:Reservoir])
+    plot_bracket(v, k)
+end
+ylims!(ax1, (-0.2*maximum(final_x), nothing))
+fig
+

examples/co2_props.jl

@@ -130,7 +130,7 @@ plot_property(:viscosity, tab_water, :CO2, "Pure phase viscosity, no salts")
 # ``y_i = K_i(p, T) x_i``
 #
 plot_property(:K, tab_water, :CO2, "K value of CO2 component in aqueous phase, no salts")
-# ## Compare viscosity with and without salts
-plot_brine_comparison(:density, 30.0, tab_water, tab_brine, :H2O, "Pure phase viscosity with salts")
 # ## Compare K value with and without salts
 plot_brine_comparison(:K, 30.0, tab_water, tab_brine, :CO2, "K value of CO2 component in aqueous phase")
+# ## Compare density with and without salts
+plot_brine_comparison(:density, 30.0, tab_water, tab_brine, :H2O, "Pure phase density")