Merge pull request #54 from sintefmath/dev

Nonlinear finite-volume schemes, constant pressure aquifers, converter for CO2STORE, simplify thermal

Artifacts.toml

@@ -0,0 +1,7 @@
+[CO2Tables_CSP11]
+git-tree-sha1 = "fe3911cdbf4fe04bdc36628f687b50b79f309e6a"
+lazy = true
+
+    [[CO2Tables_CSP11.download]]
+    sha256 = "b20fd6fbffdc9e2ab6f6e147e799ecc57fd2b8b5888923ccad3b4d437eab3479"
+    url = "https://github.com/sintefmath/JutulDarcy.jl/releases/download/v0.2.27/csp11.tar.gz"

Project.toml

@@ -1,9 +1,10 @@
 name = "JutulDarcy"
 uuid = "82210473-ab04-4dce-b31b-11573c4f8e0a"
 authors = ["Olav Møyner <olav.moyner@gmail.com>"]
-version = "0.2.27"
+version = "0.2.28"
 
 [deps]
+Artifacts = "56f22d72-fd6d-98f1-02f0-08ddc0907c33"
 AlgebraicMultigrid = "2169fc97-5a83-5252-b627-83903c6c433c"
 DataStructures = "864edb3b-99cc-5e75-8d2d-829cb0a9cfe8"
 Dates = "ade2ca70-3891-5945-98fb-dc099432e06a"
@@ -13,6 +14,7 @@ ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
 GeoEnergyIO = "3b1dd628-313a-45bb-9d8d-8f3c48dcb5d4"
 Jutul = "2b460a1a-8a2b-45b2-b125-b5c536396eb9"
 Krylov = "ba0b0d4f-ebba-5204-a429-3ac8c609bfb7"
+LazyArtifacts = "4af54fe1-eca0-43a8-85a7-787d91b784e3"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 LoopVectorization = "bdcacae8-1622-11e9-2a5c-532679323890"
 MAT = "23992714-dd62-5051-b70f-ba57cb901cac"
@@ -41,17 +43,19 @@ JutulDarcyGLMakieExt = "GLMakie"
 JutulDarcyPartitionedArraysExt = ["PartitionedArrays", "MPI", "HYPRE"]
 
 [compat]
+Artifacts = "1"
 AlgebraicMultigrid = "0.5.1, 0.6.0"
 DataStructures = "0.18.13"
 Dates = "1"
 DelimitedFiles = "1.9.1"
 DocStringExtensions = "0.9.3"
 ForwardDiff = "0.10.30"
 GLMakie = "0.10.0"
-GeoEnergyIO = "1.1.3"
+GeoEnergyIO = "1.1.5"
 HYPRE = "1.4.0"
-Jutul = "0.2.35"
+Jutul = "0.2.37"
 Krylov = "0.9.1"
+LazyArtifacts = "1"
 LinearAlgebra = "1"
 LoopVectorization = "0.12.115"
 MAT = "0.10.3"

docs/make.jl

@@ -71,6 +71,7 @@ function build_jutul_darcy_docs(build_format = nothing; build_examples = true, b
             prettyurls=get(ENV, "CI", "false") == "true",
             canonical="https://sintefmath.github.io/JutulDarcy.jl",
             edit_link="main",
+            size_threshold_ignore = ["ref/jutul.md", "docstrings.md"],
             assets=String["assets/citations.css"],
         )
     end

docs/src/man/basics/systems.md

@@ -159,6 +159,5 @@ For additional details, please see Chapter 8 - Compositional Simulation with the
 
 ```@docs
 reservoir_system
-ThermalSystem
-Jutul.CompositeSystem
+JutulDarcy.add_thermal_to_model!
 ```

examples/co2_sloped.jl

@@ -24,30 +24,120 @@ for (i, pt) in enumerate(points)
     x, y, z = pt
     x_u = 2*π*x/1000.0
     w = 0.2
-    dz = 0.05*x + w*(30*sin(2.0*x_u) + 20*sin(5.0*x_u) + 10*sin(10.0*x_u) + 5*sin(25.0*x_u))
+    dz = 0.05*x + 0.05*abs(x - 500.0)+ w*(30*sin(2.0*x_u) + 20*sin(5.0*x_u))
     points[i] = pt + [0, 0, dz]
 end
+# ## Find and plot cells intersected by a deviated injector well
+# We place a single injector well. This well was unfortunately not drilled
+# completely straight, so we cannot directly use `add_vertical_well` based on
+# logical indices. We instead define a matrix with three columns x, y, z that
+# lie on the well trajectory and use utilities from `Jutul` to find the cells
+# intersected by the trajectory.
+import Jutul: find_enclosing_cells, plot_mesh_edges
+trajectory = [
+    745.0 0.5 45;    # First point
+    760.0 0.5 70;    # Second point
+    810.0 0.5 100.0  # Third point
+]
+
+wc = find_enclosing_cells(mesh, trajectory)
+
+fig, ax, plt = plot_mesh_edges(mesh, z_is_depth = true)
+plot_mesh!(ax, mesh, cells = wc, transparency = true, alpha = 0.3)
+lines!(ax, trajectory', color = :red)
+fig
+
 # ## Set up simulation model
 # We set up a domain and a single injector. We pass the special :co2brine
 # argument in place of the system to the reservoir model setup routine. This
 # will automatically set up a compositional two-component CO2-H2O model with the
 # appropriate functions for density, viscosity and miscibility.
 #
-# Note that this model can be run with a thermal mode by setting 
-domain = reservoir_domain(mesh, permeability = 1.0Darcy, porosity = 0.3, temperature = convert_to_si(30.0, :Celsius))
-Injector = setup_well(domain, (65, 1, 1), name = :Injector)
-model, parameters = setup_reservoir_model(domain, :co2brine, wells = Injector);
-# ## Find the boundary and set increased volume
-# We find the left and right boundary of the model and increase the volume of
-# those cells. This mimicks a constant pressure boundary condition.
+# Note that this model by default is isothermal, but we still need to specify a
+# temperature when setting up the model. This is because the properties of CO2
+# strongly depend on temperature, even when thermal transport is not solved.
+domain = reservoir_domain(mesh, permeability = 0.3Darcy, porosity = 0.3, temperature = convert_to_si(30.0, :Celsius))
+Injector = setup_well(domain, wc, name = :Injector, simple_well = true)
+
+model = setup_reservoir_model(domain, :co2brine, wells = Injector, extra_out = false);
+# ## Customize model by adding relative permeability with hysteresis
+# We define three relative permeability functions: kro(so) for the brine/liquid
+# phase and krg(g) for both drainage and imbibition. Here we limit the
+# hysteresis to only the non-wetting gas phase, but either combination of
+# wetting or non-wetting hysteresis is supported.
+#
+# Note that we import a few utilities from JutulDarcy that are not exported by
+# default since hysteresis falls under advanced functionality.
+import JutulDarcy: table_to_relperm, add_relperm_parameters!, brooks_corey_relperm
+so = range(0, 1, 10)
+krog_t = so.^2
+krog = PhaseRelativePermeability(so, krog_t, label = :og)
+
+# Higher resolution for second table
+sg = range(0, 1, 50)
+
+# Evaluate Brooks-Corey to generate tables
+tab_krg_drain = brooks_corey_relperm.(sg, n = 2, residual = 0.1)
+tab_krg_imb = brooks_corey_relperm.(sg, n = 3, residual = 0.25)
+
+krg_drain  = PhaseRelativePermeability(sg, tab_krg_drain, label = :g)
+krg_imb  = PhaseRelativePermeability(sg, tab_krg_imb, label = :g)
+
+fig, ax, plt = lines(sg, tab_krg_drain, label = "krg drainage")
+lines!(ax, sg, tab_krg_imb, label = "krg imbibition")
+lines!(ax, 1 .- so, krog_t, label = "kro")
+axislegend()
+fig
+## Define a relative permeability variable
+# JutulDarcy uses type instances to define how different variables inside the
+# simulation are evaluated. The `ReservoirRelativePermeabilities` type has
+# support for up to three phases with w, ow, og and g relative permeabilities
+# specified as a function of their respective phases. It also supports
+# saturation regions.
+#
+# Note: If regions are used, all drainage curves come first followed by equal
+# number of imbibition curves. Since we only have a single (implicit) saturation
+# region, the krg input should have two entries: One for drainage, and one for
+# imbibition.
+#
+# We also call `add_relperm_parameters` to the model. This makes sure that when
+# hysteresis is enabled, we track maximum saturation for hysteresis in each
+# reservoir cell.
+import JutulDarcy: KilloughHysteresis, ReservoirRelativePermeabilities
+krg = (krg_drain, krg_imb) 
+H_g = KilloughHysteresis() # Other options: CarlsonHysteresis, JargonHysteresis
+relperm = ReservoirRelativePermeabilities(g = krg, og = krog, hysteresis_g = H_g)
+replace_variables!(model, RelativePermeabilities = relperm)
+add_relperm_parameters!(model)
+# ## Define approximate hydrostatic pressure
+# The initial pressure of the water-filled domain is assumed to be at
+# hydrostatic equilibrium.
+nc = number_of_cells(mesh)
+p0 = zeros(nc)
+depth = domain[:cell_centroids][3, :]
+g = Jutul.gravity_constant
+@. p0 = 200bar + depth*g*1000.0
+# ## Set up initial state and parameters
+state0 = setup_reservoir_state(model,
+    Pressure = p0,
+    OverallMoleFractions = [1.0, 0.0],
+)
+parameters = setup_parameters(model)
+
+# ## Find the boundary and apply a constant pressureboundary condition
+# We find cells on the left and right boundary of the model and set a constant
+# pressure boundary condition to represent a bounding aquifer that retains the
+# initial pressure far away from injection.
+
 boundary = Int[]
-for cell in 1:number_of_cells(mesh)
+for cell in 1:nc
     I, J, K = cell_ijk(mesh, cell)
     if I == 1 || I == nx
         push!(boundary, cell)
     end
 end
-parameters[:Reservoir][:FluidVolume][boundary] *= 1000;
+bc = flow_boundary_condition(boundary, domain, p0[boundary], fractional_flow = [1.0, 0.0])
+
 # ## Plot the model
 plot_reservoir(model)
 # ## Set up schedule
@@ -58,37 +148,41 @@ nstep = 25
 nstep_shut = 25
 dt_inject = fill(365.0day, nstep)
 pv = pore_volume(model, parameters)
-inj_rate = 0.0075*sum(pv)/sum(dt_inject)
+inj_rate = 0.05*sum(pv)/sum(dt_inject)
 
 rate_target = TotalRateTarget(inj_rate)
 I_ctrl = InjectorControl(rate_target, [0.0, 1.0],
     density = 900.0,
 )
 # Set up forces for use in injection
 controls = Dict(:Injector => I_ctrl)
-forces_inject = setup_reservoir_forces(model, control = controls)
+forces_inject = setup_reservoir_forces(model, control = controls, bc = bc)
 # Forces with shut wells
-forces_shut = setup_reservoir_forces(model)
+forces_shut = setup_reservoir_forces(model, bc = bc)
 dt_shut = fill(365.0day, nstep_shut);
 # Combine the report steps and forces into vectors of equal length
 dt = vcat(dt_inject, dt_shut)
-forces = vcat(fill(forces_inject, nstep), fill(forces_shut, nstep_shut));
-# ## Set up initial state
-state0 = setup_reservoir_state(model,
-    Pressure = 200bar,
-    OverallMoleFractions = [1.0, 0.0],
-)
+forces = vcat(
+    fill(forces_inject, nstep),
+    fill(forces_shut, nstep_shut)
+);
+# ## Add some more outputs for plotting
+rmodel = reservoir_model(model)
+push!(rmodel.output_variables, :RelativePermeabilities)
+push!(rmodel.output_variables, :PhaseViscosities)
 # ## Simulate the schedule
 # We set a maximum internal time-step of 30 days to ensure smooth convergence
 # and reduce numerical diffusion.
 wd, states, t = simulate_reservoir(state0, model, dt,
     parameters = parameters,
     forces = forces,
-    max_timestep = 30day
+    max_timestep = 90day
 )
-# ## Plot the density of brine
-# The density of brine depends on the CO2 concentration and gives a good
-# visualization of where the mass of CO2 exists.
+# ## 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.
 using GLMakie
 function plot_co2!(fig, ix, x, title = "")
     ax = Axis3(fig[ix, 1],
@@ -102,8 +196,33 @@ function plot_co2!(fig, ix, x, title = "")
 end
 fig = Figure(size = (900, 1200))
 for (i, step) in enumerate([1, 5, nstep, nstep+nstep_shut])
-    plot_co2!(fig, i, states[step][:PhaseMassDensities][1, :], "Brine density report step $step/$(nstep+nstep_shut)")
+    plot_co2!(fig, i, log10.(states[step][:OverallMoleFractions][2, :]), "log10 of CO2 mole fraction at report step $step/$(nstep+nstep_shut)")
+end
+fig
+# ## Plot all relative permeabilities for all time-steps
+# We can plot all relative permeability evaluations. This both verifies that the
+# hysteresis model is active, but also gives an indication to how many cells are
+# exhibiting imbibition during the simulation.
+kro_val = Float64[]
+krg_val = Float64[]
+sg_val = Float64[]
+for state in states
+    kr_state = state[:RelativePermeabilities]
+    s_state = state[:Saturations]
+    for c in 1:nc
+        push!(kro_val, kr_state[1, c])
+        push!(krg_val, kr_state[2, c])
+        push!(sg_val, s_state[2, c])
+    end
 end
+
+fig = Figure()
+ax = Axis(fig[1, 1], title = "Relative permeability during simulation")
+fig, ax, plt = scatter(sg_val, kro_val, label = "kro", alpha = 0.3)
+scatter!(ax, sg_val, krg_val, label = "krg", alpha = 0.3)
+axislegend()
 fig
 # ## Plot result in interactive viewer
-plot_reservoir(model, states)
+# If you have interactive plotting available, you can explore the results
+# yourself.
+plot_reservoir(model, states)

src/CO2Properties/CO2Properties.jl

@@ -1,5 +1,5 @@
 module CO2Properties
-    using Jutul, JutulDarcy, StaticArrays, MultiComponentFlash, DelimitedFiles
+    using Jutul, JutulDarcy, StaticArrays, MultiComponentFlash, DelimitedFiles, Artifacts, LazyArtifacts
     include("kvalues.jl")
     include("props.jl")
     include("reading.jl")

src/CO2Properties/reading.jl

@@ -1,12 +1,13 @@
 function read_solubility_table(name; fix = true)
     x, header = DelimitedFiles.readdlm(name, ',', header = true, skipstart = 7)
+    header = map(strip, header)
     if fix
         new_header = Symbol[]
         lookup = Dict(
             "# temperature [°C]" => :T,
-            " phase pressure [Pa]" => :p,
-            "         x_CO2 [-]" => :x_co2,
-            "         y_H2O [-]" => :y_h2o,
+            "phase pressure [Pa]" => :p,
+            "x_CO2 [mol/mol]" => :x_co2,
+            "y_H2O [mol/mol]" => :y_h2o,
         )
         for label in header
             push!(new_header, lookup[label])