Merge pull request #58 from sintefmath/numerical-aquifers

Support for numerical aquifers and cleanup of capillary pressure
sintefmath · Sep 22, 2024 · faf16d9 · faf16d9 · moyner · Sep 22, 2024
2 parents 8ecc3cb + 72b7206
commit faf16d9
Show file tree

Hide file tree

Showing 15 changed files with 457 additions and 103 deletions.
diff --git a/Project.toml b/Project.toml
@@ -1,7 +1,7 @@
 name = "JutulDarcy"
 uuid = "82210473-ab04-4dce-b31b-11573c4f8e0a"
 authors = ["Olav Møyner <[email protected]>"]
-version = "0.2.30"
+version = "0.2.31"
 
 [deps]
 Artifacts = "56f22d72-fd6d-98f1-02f0-08ddc0907c33"
@@ -53,7 +53,7 @@ DelimitedFiles = "1.9.1"
 DocStringExtensions = "0.9.3"
 ForwardDiff = "0.10.30"
 GLMakie = "0.10.0"
-GeoEnergyIO = "1.1.6"
+GeoEnergyIO = "1.1.9"
 HYPRE = "1.4.0"
 Jutul = "0.2.37"
 Krylov = "0.9.1"

diff --git a/docs/src/index.md b/docs/src/index.md
@@ -11,7 +11,22 @@ DocTestSetup = quote
 end
 ```
 
-Documentation for [JutulDarcy.jl](https://github.com/sintefmath/JutulDarcy.jl). The documentation is currently limited to docstrings and a series of examples. The examples are sorted by complexity. We suggest you start with [Gravity segregation example](Gravity segregation example).
+Documentation for [JutulDarcy.jl](https://github.com/sintefmath/JutulDarcy.jl). The documentation consists of a mixture of technical documentation, docstrings and examples. The examples are sorted by complexity. We suggest you start with [Gravity segregation example](Gravity segregation example).
+
+To generate a folder that contains the examples locally, you can run the following code to create a folder `jutuldarcy_examples` in your current working directory:
+
+```julia
+using JutulDarcy
+generate_jutuldarcy_examples()
+```
+
+Alternatively, a folder can be specified if you want the examples to be placed outside your present working directory:
+
+
+```julia
+using JutulDarcy
+generate_jutuldarcy_examples("/home/username/")
+```
 
 JutulDarcy builds upon the general features found in [Jutul.jl](https://github.com/sintefmath/Jutul.jl). You may also find it useful to look at the [Jutul.jl documentation](https://sintefmath.github.io/Jutul.jl/dev/).
 

diff --git a/examples/co2_brine_2d_vertical.jl b/examples/co2_brine_2d_vertical.jl
@@ -19,7 +19,6 @@ bic = zeros(2, 2)
 
 mixture = MultiComponentMixture([h2o, co2], A_ij = bic, names = ["H2O", "CO2"])
 eos = GenericCubicEOS(mixture, PengRobinson())
-#-
 # ## Set up domain and wells
 using Jutul, JutulDarcy, GLMakie
 nx = 50
@@ -35,7 +34,6 @@ res = reservoir_domain(g, porosity = 0.3, permeability = K)
 prod = setup_well(g, K, [(nx, ny, 1)], name = :Producer)
 # Set up an injector in the opposite corner, perforated in bottom layer
 inj = setup_well(g, K, [(1, 1, nz)], name = :Injector)
-#-
 # ## Define system and realize on grid
 rhoLS = 844.23*kg/meter^3
 rhoVS = 126.97*kg/meter^3
@@ -46,15 +44,15 @@ model, parameters = setup_reservoir_model(res, sys, wells = [inj, prod]);
 push!(model[:Reservoir].output_variables, :Saturations)
 kr = BrooksCoreyRelativePermeabilities(sys, 2.0, 0.0, 1.0)
 model = replace_variables!(model, RelativePermeabilities = kr)
-T0 = repeat([303.15*Kelvin], 1, nc)
+T0 = fill(303.15*Kelvin, nc)
 parameters[:Reservoir][:Temperature] = T0
 state0 = setup_reservoir_state(model, Pressure = 50*bar, OverallMoleFractions = [1.0, 0.0]);
 
 # ## Define schedule
 # 5 year (5*365.24 days) simulation period
-dt0 = repeat([1]*day, 26)
-dt1 = repeat([10.0]*day, 180)
-dt = append!(dt0, dt1)
+dt0 = fill(1*day, 26)
+dt1 = fill(10.0*day, 180)
+dt = cat(dt0, dt1, dims = 1)
 rate_target = TotalRateTarget(9.5066e-06*meter^3/sec)
 I_ctrl = InjectorControl(rate_target, [0, 1], density = rhoVS)
 bhp_target = BottomHolePressureTarget(50*bar)
@@ -78,18 +76,17 @@ function plot_vertical(x, t)
     Colorbar(fig[1, 2], plot)
     fig
 end
-#-
+
 # ### Plot final CO2 mole fraction
 plot_vertical(z, "CO2")
-#-
+
 # ### Plot final vapor saturation 
 sg = states[end][:Saturations][2, :]
 plot_vertical(sg, "Vapor saturation")
-#-
+
 # ### Plot final pressure
 p = states[end][:Pressure]
 plot_vertical(p./bar, "Pressure [bar]")
 
-#-
 # ### Plot in interactive viewer
 plot_reservoir(model, states, step = length(dt), key = :Saturations)
diff --git a/examples/co2_sloped.jl b/examples/co2_sloped.jl
@@ -23,7 +23,8 @@ Darcy, bar, kg, meter, day, yr = si_units(:darcy, :bar, :kilogram, :meter, :day,
 # surface has a large impact on where the CO2 migrates.
 cart_dims = (nx, 1, nz)
 physical_dims = (1000.0, 1.0, 50.0)
-mesh = UnstructuredMesh(CartesianMesh(cart_dims, physical_dims))
+cart_mesh = CartesianMesh(cart_dims, physical_dims)
+mesh = UnstructuredMesh(cart_mesh)
 
 points = mesh.node_points
 for (i, pt) in enumerate(points)
@@ -200,9 +201,9 @@ println("Boundary condition added to $(length(bc)) cells.")
 # ## Plot the model
 plot_reservoir(model)
 # ## Set up schedule
-# We set up 25 years of injection and 1000 years of migration where the well is
-# shut. The density of the injector is set to 900 kg/m^3, which is roughly
-# the density of CO2 at in-situ conditions.
+# We set up 25 years of injection and 475 years of migration where the well is
+# shut. The density of the injector is set to 900 kg/m^3, which is roughly the
+# density of CO2 at these high-pressure in-situ conditions.
 nstep = 25
 nstep_shut = 475
 dt_inject = fill(365.0day, nstep)
@@ -239,10 +240,10 @@ wd, states, t = simulate_reservoir(state0, model, dt,
     max_timestep = 90day
 )
 # ## Plot the CO2 mole fraction
-# We plot log10 of the CO2 mole fraction. We use log10 to account for the fact
-# that the mole fraction in cells made up of only the aqueous phase is much
-# smaller than that of cells with only the gaseous phase, where there is almost
-# just CO2.
+# We plot the overall CO2 mole fraction. We scale the color range to account for
+# the fact that the mole fraction in cells made up of only the aqueous phase is
+# much smaller than that of cells with only the gaseous phase, where there is
+# almost just CO2.
 using GLMakie
 function plot_co2!(fig, ix, x, title = "")
     ax = Axis3(fig[ix, 1],
@@ -251,15 +252,15 @@ function plot_co2!(fig, ix, x, title = "")
         elevation = 0.05,
         aspect = (1.0, 1.0, 0.3),
         title = title)
-    plt = plot_cell_data!(ax, mesh, x, colormap = :seaborn_icefire_gradient)
+    plt = plot_cell_data!(ax, mesh, x, colormap = :seaborn_icefire_gradient, colorrange = (0.0, 0.1))
     Colorbar(fig[ix, 2], plt)
 end
 fig = Figure(size = (900, 1200))
 for (i, step) in enumerate([1, 5, nstep, nstep+nstep_shut])
     if use_immiscible
         plot_co2!(fig, i, states[step][:Saturations][2, :], "CO2 plume saturation at report step $step/$(nstep+nstep_shut)")
     else
-        plot_co2!(fig, i, log10.(states[step][:OverallMoleFractions][2, :]), "log10 of CO2 mole fraction at report step $step/$(nstep+nstep_shut)")
+        plot_co2!(fig, i, states[step][:OverallMoleFractions][2, :], "CO2 mole fraction at report step $step/$(nstep+nstep_shut)")
     end
 end
 fig

diff --git a/src/deck_types.jl b/src/deck_types.jl
@@ -620,6 +620,15 @@ struct TableCompressiblePoreVolume{V, R} <: ScalarVariable
     end
 end
 
+function Jutul.subvariable(p::TableCompressiblePoreVolume, map::FiniteVolumeGlobalMap)
+    c = map.cells
+    regions = Jutul.partition_variable_slice(p.regions, c)
+    return TableCompressiblePoreVolume(
+        p.tab,
+        regions = regions
+    )
+end
+
 struct ScalarPressureTable{V, R} <: ScalarVariable
     tab::V
     regions::R
@@ -630,6 +639,15 @@ struct ScalarPressureTable{V, R} <: ScalarVariable
     end
 end
 
+function Jutul.subvariable(p::ScalarPressureTable, map::FiniteVolumeGlobalMap)
+    c = map.cells
+    regions = Jutul.partition_variable_slice(p.regions, c)
+    return ScalarPressureTable(
+        p.tab,
+        regions = regions
+    )
+end
+
 @jutul_secondary function update_variable!(pv, Φ::ScalarPressureTable, model, Pressure, ix)
     @inbounds for i in ix
         reg = region(Φ.regions, i)

diff --git a/src/facility/wells/wells.jl b/src/facility/wells/wells.jl
@@ -63,10 +63,21 @@ function common_well_setup(nr; dz = nothing, WI = nothing, gravity = gravity_con
 end
 
 """
-    setup_well(D::DataDomain, reservoir_cells; skin = 0.0, Kh = nothing, radius = 0.1, dir = :z)
+    setup_well(D::DataDomain, reservoir_cells; skin = 0.0, Kh = nothing, radius = 0.1, dir = :z, name = :Well)
+    w = setup_well(D, 1, name = :MyWell)         # Cell 1 in the grid
+    w = setup_well(D, (2, 5, 1), name = :MyWell) # Position (2, 5, 1) in logically structured mesh
+    w2 = setup_well(D, [1, 2, 3], name = :MyOtherWell)
+
 
 Set up a well in `reservoir_cells` with given skin factor and radius. The order
 of cells matter as it is treated as a trajectory.
+
+The `name` keyword argument can be left defaulted if your model will only have a
+single well (named `:Well`). It is highly recommended to provide this whenever a
+well is set up.
+
+`reservoir_cells` can be one of the following: A Vector of cells, a single cell,
+a Vector of `(I, J, K)` Tuples or a single Tuple of the same type.
 """
 function setup_well(D::DataDomain, reservoir_cells; cell_centers = D[:cell_centroids], kwarg...)
     K = D[:permeability]
@@ -205,10 +216,12 @@ function Jutul.plot_primitives(mesh::MultiSegmentWell, plot_type; kwarg...)
 end
 
 """
-    setup_vertical_well(D::DataDomain, i, j; <kwarg>)
+    setup_vertical_well(D::DataDomain, i, j; name = :MyWell, <kwarg>)
 
-Set up a vertical well with a `DataDomain` input that represents the porous
-medium / reservoir where the wells it to be placed.
+Set up a vertical well with a [`DataDomain`](@ref) input that represents the porous
+medium / reservoir where the wells it to be placed. See [`SimpleWell`](@ref),
+[`MultiSegmentWell`](@ref) and [`setup_well`](@ref) for more details about possible keyword
+arguments.
 """
 function setup_vertical_well(D::DataDomain, i, j; cell_centers = D[:cell_centroids], kwarg...)
     K = D[:permeability]

diff --git a/src/init/init.jl b/src/init/init.jl
@@ -152,7 +152,7 @@ function equilibriate_state!(init, depths, model, sys, contacts, depth, datum_pr
         return phase_density
     end
     # Find the reference phase. It is either liquid
-    ref_phase = find_reference_phase(model.system)
+    ref_phase = get_reference_phase_index(model.system)
     pressures = determine_hydrostatic_pressures(depths, depth, zmin, zmax, contacts, datum_pressure, density_f, contacts_pc, ref_phase)
     if nph > 1
         relperm = model.secondary_variables[:RelativePermeabilities]
@@ -234,42 +234,21 @@ function equilibriate_state!(init, depths, model, sys, contacts, depth, datum_pr
     return init
 end
 
-function find_reference_phase(system)
-    mphases = get_phases(system)
-    function find_phase(k)
-        ix = 0
-        for (i, ph) in enumerate(mphases)
-            if ph isa LiquidPhase
-                ix = i
-                break
-            end
-        end
-        return ix
-    end
-    l = find_phase(LiquidPhase)
-    a = find_phase(AqueousPhase)
-    if l > 0
-        phase = l
-    elseif a > 0
-        phase = a
-    else
-        # Last phase is water or something custom, send out 1.
-        phase = 1
-    end
-    return phase
-end
-
-
 function parse_state0_equil(model, datafile; normalize = :sum)
+    sys = model.system
+    d = model.data_domain
+
     has_water = haskey(datafile["RUNSPEC"], "WATER")
-    has_oil = haskey(datafile["RUNSPEC"], "OIL")
-    has_gas = haskey(datafile["RUNSPEC"], "GAS")
+    if sys isa CompositionalSystem
+        has_oil = has_gas = true
+    else
+        has_oil = haskey(datafile["RUNSPEC"], "OIL")
+        has_gas = haskey(datafile["RUNSPEC"], "GAS")
+    end
 
     is_co2 = haskey(datafile["RUNSPEC"], "JUTUL_CO2BRINE")
     is_single_phase = (has_water + has_oil + has_gas) == 1
 
-    sys = model.system
-    d = model.data_domain
     has_sat_reg = haskey(d, :satnum)
     ncells = number_of_cells(d)
     if has_sat_reg
@@ -392,9 +371,6 @@ function parse_state0_equil(model, datafile; normalize = :sum)
                             cap = cap[2:end]
                         end
                         ix = unique(i -> cap[i], 1:length(cap))
-                        if i == 1 && get_phases(model.system)[1] isa AqueousPhase
-                            @. cap *= -1
-                        end
                         s = s[ix]
                         cap = cap[ix]
                         if length(s) == 1
@@ -446,11 +422,11 @@ function parse_state0_equil(model, datafile; normalize = :sum)
                         contacts_pc = (goc_pc, )
                     else
                         contacts = (woc, )
-                        contacts_pc = (-woc_pc, )
+                        contacts_pc = (woc_pc, )
                     end
                 else
                     contacts = (woc, goc)
-                    contacts_pc = (-woc_pc, goc_pc)
+                    contacts_pc = (woc_pc, goc_pc)
                 end
 
                 if disgas || vapoil