Merge pull request #66 from sintefmath/dev

Rock hysteresis, fix scaled capillary pressure regression, bump HYPRE

Project.toml

@@ -1,7 +1,7 @@
 name = "JutulDarcy"
 uuid = "82210473-ab04-4dce-b31b-11573c4f8e0a"
 authors = ["Olav Møyner <[email protected]>"]
-version = "0.2.33"
+version = "0.2.34"
 
 [deps]
 AlgebraicMultigrid = "2169fc97-5a83-5252-b627-83903c6c433c"
@@ -55,16 +55,16 @@ DocStringExtensions = "0.9.3"
 ForwardDiff = "0.10.30"
 GLMakie = "0.10.0"
 GeoEnergyIO = "1.1.12"
-HYPRE = "1.4.0"
-Jutul = "0.2.37"
+HYPRE = "1.6.0"
+Jutul = "0.2.39"
 Krylov = "0.9.1"
 LazyArtifacts = "1"
 LinearAlgebra = "1"
 LoopVectorization = "0.12.115"
 MAT = "0.10.3"
 MultiComponentFlash = "1.1.15"
 OrderedCollections = "1.6.2"
-PartitionedArrays = "0.3.2"
+PartitionedArrays = "0.5"
 Polyester = "0.6.11, 0.7.3"
 Polynomials = "3, 4"
 PrecompileTools = "1.0.1"

docs/src/man/basics/systems.md

@@ -120,10 +120,15 @@ StandardBlackOilSystem
 
 The black-oil equations is an extension of the immiscible description to handle limited miscibility between the phases. Originally developed for certain types of oil and gas simulation, these equations are useful when the number of components is low and tabulated values for dissolution and vaporization are available.
 
-The assumptions of the black-oil model is that the "oil" and "gas" pseudo-components have uniform composition throughout the domain. JutulDarcy supports two- and three-phase black oil flow. The difference between two and three phases amounts to an additional immiscible aqueous phase that is identical to that of the previous section. For that reason, we focus on the miscible pseudo-components:
+The assumptions of the black-oil model is that the "oil" and "gas" pseudo-components have uniform composition throughout the domain. JutulDarcy supports two- and three-phase black oil flow. The difference between two and three phases amounts to an additional immiscible aqueous phase that is identical to that of the previous section. For that reason, we focus on the miscible pseudo-components, oil:
+
+```math
+r_o = \rho_o^s \left( \frac{\partial}{\partial t}( (b_o S_o + R_v b_g (1 - S_o)) \phi) + \nabla \cdot ( b_o \vec{v}_o + R_v b_o \vec{v}_g) - q_o^s \right ),
+```
+
+and gas,
 
 ```math
-r_o = \rho_o^s \left( \frac{\partial}{\partial t}( (b_o S_o + R_v b_g (1 - S_o)) \phi) + \nabla \cdot ( b_o \vec{v}_o + R_v b_o \vec{v}_g) - q_o^s \right ) \\
 r_g = \rho_g^s \left( \frac{\partial}{\partial t}( (b_g S_g + R_s b_o S_o) \phi) + \nabla \cdot ( b_g \vec{v}_g + R_s b_g \vec{v}_o) - q_g^s \right )
 ```
 

examples/intro_sensitivities.jl

@@ -138,7 +138,7 @@ controls[:Producer] = P_ctrl
 
 forces = setup_reservoir_forces(model, control = controls)
 case = JutulCase(model, dt, forces, parameters = parameters, state0 = state0)
-result = simulate_reservoir(case, output_substates = true);
+result = simulate_reservoir(case, output_substates = true, info_level = -1);
 ##
 ws, states = result
 ws(:Producer, :grat)

examples/two_phase_unstable_gravity.jl

@@ -62,7 +62,7 @@ state0 = setup_reservoir_state(model, Pressure = p0, Saturations = s0)
 @. μ[2, :] = 5e-3
 # Convert time-steps from days to seconds
 timesteps = repeat([10.0*3600*24], 20)
-_, states, = simulate_reservoir(state0, model, timesteps, parameters = parameters, info_level = 1);
+_, states, = simulate_reservoir(state0, model, timesteps, parameters = parameters, info_level = -1);
 # ## Plot results
 # Plot initial saturation
 tmp = reshape(state0[:Reservoir][:Saturations][1, :], nx, nz)

examples/validation_compositional.jl

@@ -1,4 +1,4 @@
-# # Validation of equation-of-state compositiona simulator
+# # Validation of equation-of-state compositional simulator
 # This example solves a 1D two-phase, three component miscible displacement
 # problem and compares against existing simulators (E300, AD-GPRS) to verify
 # correctness.
@@ -25,7 +25,7 @@ using GLMakie
 dpth = JutulDarcy.GeoEnergyIO.test_input_file_path("SIMPLE_COMP")
 data_path = joinpath(dpth, "SIMPLE_COMP.DATA")
 case = setup_case_from_data_file(data_path)
-result = simulate_reservoir(case, info_level = 1)
+result = simulate_reservoir(case, info_level = -1)
 ws, states = result;
 # ## Plot solutions and compare
 # The 1D displacement can be plotted as a line plot. We pick a step midway

examples/validation_norne_nohyst.jl

@@ -8,7 +8,7 @@ using Jutul, JutulDarcy, GLMakie, DelimitedFiles, HYPRE, GeoEnergyIO
 norne_dir = GeoEnergyIO.test_input_file_path("NORNE_NOHYST")
 data_pth = joinpath(norne_dir, "NORNE_NOHYST.DATA")
 data = parse_data_file(data_pth)
-case = setup_case_from_data_file(data)
+case = setup_case_from_data_file(data);
 # ## Unpack the case to see basic data structures
 model = case.model
 parameters = case.parameters

src/CO2Properties/generation.jl

@@ -11,14 +11,14 @@ P <= 35 MPa
 
 # Arguments:
 
-    - T: Scalar with temperature value in Kelvin
-    - P: Scalar with pressure value in bar
-    - S: Salt mass fraction
+  - T: Scalar with temperature value in Kelvin
+  - P: Scalar with pressure value in bar
+  - S: Salt mass fraction
 
 # Outputs
 
-    - rho_brine:    Scalar with density value in kg/m3
-    - c_brine:      Scalar with compressibility value in 1/kPa
+  - rho_brine:    Scalar with density value in kg/m3
+  - c_brine:      Scalar with compressibility value in 1/kPa
 
 """
 function pvt_brine_RoweChou1970(T, P, S; check::Bool = true)
@@ -67,12 +67,12 @@ Hassanzadeh et al., IJGGC (2008).
 
 # Arguments:
 
-    - T: Scalar of temperature value in Kelvin
-    - rho: Scalar of density value in kg/m^3
+  - T: Scalar of temperature value in Kelvin
+  - rho: Scalar of density value in kg/m^3
 
 # Outputs
 
-    - mu: Dynamic viscosity in Pa*s
+  - mu: Dynamic viscosity in Pa*s
 """
 function viscosity_co2_Fenghour1998(T, rho; check = true)
 
@@ -119,12 +119,13 @@ These authors used the data of Rowe and Chou (1970), Zarembo & Fedorov
 P valid from 5 to 100 MPa, T from 20 to 350 C (Adams & Bachu, 2002)
 
 # Arguments
-    - T: Temperature value in Kelvin
-    - P: Pressure value in bar
-    - w_nacl: Salt (NaCl) mass fraction
+
+  - T: Temperature value in Kelvin
+  - P: Pressure value in bar
+  - w_nacl: Salt (NaCl) mass fraction
 
 # Outputs 
-    - rho_b: Scalar with brine density in kg/m3
+  - rho_b: Scalar with brine density in kg/m3
 """
 function pvt_brine_BatzleWang1992(T, P, w_nacl; check::Bool = true)
 
@@ -175,19 +176,19 @@ to 600 bar, for (1) CO2 compressibility factor, (2) CO2 fugacity coefficient and
 
 # Arguments
 
-    - T: Scalar with temperature value in Kelvin
-    - P: Scalar with pressure value in bar
+  - T: Scalar with temperature value in Kelvin
+  - P: Scalar with pressure value in bar
 
 # Optional arguments:
 
-    - a_m: Intermolecular attraction constant (of the mixture) in bar*cm^6*K^0.5/mol^2
-    - b_m: Intermolecular repulsion constant (of the mixture) in cm^3/mol
+  - a_m: Intermolecular attraction constant (of the mixture) in bar*cm^6*K^0.5/mol^2
+  - b_m: Intermolecular repulsion constant (of the mixture) in cm^3/mol
 
 # Outputs
 
-    - V_m: Scalar with molar volume in [cm^3/mol]
-    - rhox: Scalar with density in [mol/m^3]
-    - rho: Scalar with density in [kg/m^3]
+  - V_m: Scalar with molar volume in [cm^3/mol]
+  - rhox: Scalar with density in [mol/m^3]
+  - rho: Scalar with density in [kg/m^3]
 """
 function pvt_co2_RedlichKwong1949(T, P, a_m=7.54 * 10^7 - 4.13 * 10^4 * T, b_m=27.8; check = true)
     # Check if T, P conditions are within range
@@ -266,17 +267,17 @@ Calculate a CO2 pseudo activity coefficient based on a virial expansion
 of excess Gibbs energy.
 
 # Arguments
-    - T: Scalar with temperature value in Kelvin
-    - P: Scalar with pressure value in bar
-    - m_io: Vector where each entry corresponds to the
+  - T: Scalar with temperature value in Kelvin
+  - P: Scalar with pressure value in bar
+  - m_io: Vector where each entry corresponds to the
     molality of a particular ion in the initial brine solution. The
     order is as follows:
     [ Na(+),   K(+),  Ca(2+), Mg(2+), Cl(-), SO4(2-)]
 
 # Outputs 
-    - V_m: Scalar with molar volume in [cm^3/mol]
-    - rhox: Scalar with density in [mol/m^3]
-    - rho: Scalar with density in [kg/m^3]
+  - V_m: Scalar with molar volume in [cm^3/mol]
+  - rhox: Scalar with density in [mol/m^3]
+  - rho: Scalar with density in [kg/m^3]
 """
 function activity_co2_DS2003(T, P, m_io; check = true)
     if check
@@ -351,7 +352,6 @@ function compute_salinity(inp_mole_fractions=Float64[], names=String[])
     # mass_salt = sum(mass_fractions)
     # 1000*mass_salt/(1.0 - mass_salt)
     w_salt = mass_salt ./ (h2o_mass + mass_salt)
-
     # TODO: Check this part.
     # m_salt is molality of species k in solution (mol k / kg solvent)
     # % w_k          mass fraction of species k in aqueous phase
@@ -381,14 +381,14 @@ who provided experimental data up to P = 200 bar, T extrapolated to
 
 # Arguments
 
-    - T: Scalar with temperature value in Kelvin
-    - P: Scalar with pressure value in bar
-    - m_nacl: Salt molality (NaCl) in mol/kg solvent
-    - w_co2: Mass fraction of CO2 in the aqueous solution (i.e. brine)
+  - T: Scalar with temperature value in Kelvin
+  - P: Scalar with pressure value in bar
+  - m_nacl: Salt molality (NaCl) in mol/kg solvent
+  - w_co2: Mass fraction of CO2 in the aqueous solution (i.e. brine)
 
 # Outputs
 
-    - mu_b_co2: Scalar with dynamic viscosity in Pa*s
+  - mu_b_co2: Scalar with dynamic viscosity in Pa*s
 
 """
 function viscosity_brine_co2_mixture_IC2012(T, P, m_nacl, w_co2; check=true)
@@ -474,16 +474,16 @@ viscosities are valid.
 
 Arguments:
 
-    - x: Mole fraction of each component
-    - M: Molar mass of each component
-    - mu: Viscosity of each component in user chosen units
+  - x: Mole fraction of each component
+  - M: Molar mass of each component
+  - mu: Viscosity of each component in user chosen units
 
 Each input should be a Float64 Vector of length n where n is the total
 number of components
 
 # Outputs
 
-    - viscMixture: Scalar viscosity of the mixture in same units as mu
+  - viscMixture: Scalar viscosity of the mixture in same units as mu
 """
 function viscosity_gas_mixture_Davidson1993(x, M, mu)
 

src/deck_support.jl

@@ -120,6 +120,15 @@ end
     end
 end
 
+@jutul_secondary function update_pore_volume!(pv, Φ::HystereticTableCompressiblePoreVolume, model, MaxPressure, Pressure, StaticFluidVolume, ix)
+    @inbounds for i in ix
+        reg = region(Φ.regions, i)
+        F = table_by_region(Φ.tab, reg)
+        p = max(Pressure[i], MaxPressure[i])
+        pv[i] = StaticFluidVolume[i]*F(p)
+    end
+end
+
 function Jutul.line_plot_data(model::SimulationModel, k::DeckPhaseVariables)
     deck_pvt_type(::DeckShrinkageFactors) = :shrinkage
     deck_pvt_type(::DeckPhaseMassDensities) = :density

src/deck_types.jl

@@ -629,6 +629,25 @@ function Jutul.subvariable(p::TableCompressiblePoreVolume, map::FiniteVolumeGlob
     )
 end
 
+struct HystereticTableCompressiblePoreVolume{V, R} <: ScalarVariable
+    tab::V
+    regions::R
+    function HystereticTableCompressiblePoreVolume(tab; regions = nothing)
+        check_regions(regions, length(tab))
+        tab = region_wrap(tab, regions)
+        new{typeof(tab), typeof(regions)}(tab, regions)
+    end
+end
+
+function Jutul.subvariable(p::HystereticTableCompressiblePoreVolume, map::FiniteVolumeGlobalMap)
+    c = map.cells
+    regions = Jutul.partition_variable_slice(p.regions, c)
+    return HystereticTableCompressiblePoreVolume(
+        p.tab,
+        regions = regions
+    )
+end
+
 struct ScalarPressureTable{V, R} <: ScalarVariable
     tab::V
     regions::R
@@ -657,6 +676,35 @@ end
     end
 end
 
+struct HystereticScalarPressureTable{V, R} <: ScalarVariable
+    tab::V
+    regions::R
+    function HystereticScalarPressureTable(tab; regions = nothing)
+        check_regions(regions, length(tab))
+        tab = region_wrap(tab, regions)
+        new{typeof(tab), typeof(regions)}(tab, regions)
+    end
+end
+
+function Jutul.subvariable(p::HystereticScalarPressureTable, map::FiniteVolumeGlobalMap)
+    c = map.cells
+    regions = Jutul.partition_variable_slice(p.regions, c)
+    return HystereticScalarPressureTable(
+        p.tab,
+        regions = regions
+    )
+end
+
+
+@jutul_secondary function update_variable!(pv, Φ::HystereticScalarPressureTable, model, MaxPressure, Pressure, ix)
+    @inbounds for i in ix
+        reg = region(Φ.regions, i)
+        F = table_by_region(Φ.tab, reg)
+        p = max(Pressure[i], MaxPressure[i])
+        pv[i] = F(p)
+    end
+end
+
 function add_lower_pvto(data, pos, rs)
     ref_p = 101325.0
     first_offset = pos[2]-1

src/input_simulation/mrst_input.jl

@@ -823,7 +823,7 @@ function set_deck_specialization!(model, runspec, props, satnum, oil, water, gas
         set_deck_relperm!(svar, param, sys, runspec, props; oil = oil, water = water, gas = gas, satnum = satnum)
         set_deck_pc!(svar, param, sys, props; oil = oil, water = water, gas = gas, satnum = satnum, is_co2 = is_co2)
     end
-    set_deck_pvmult!(svar, param, sys, props, model.data_domain)
+    set_deck_pvmult!(svar, runspec, param, sys, props, model.data_domain)
 end
 
 function set_thermal_deck_specialization!(model, props, pvtnum, oil, water, gas)
@@ -888,7 +888,7 @@ function set_deck_relperm!(vars, param, sys, runspec, props; kwarg...)
     add_relperm_parameters!(param, kr)
 end
 
-function set_deck_pvmult!(vars, param, sys, props, reservoir)
+function set_deck_pvmult!(vars, runspec, param, sys, props, reservoir)
     # Rock compressibility (if present)
     if haskey(reservoir, :rocknum)
         regions = reservoir[:rocknum]
@@ -902,9 +902,20 @@ function set_deck_pvmult!(vars, param, sys, props, reservoir)
     if haskey(props, "ROCKTAB")
         rt = vec(props["ROCKTAB"])
         tab = map(x -> get_1d_interpolator(x[:, 1], x[:, 2]), rt)
-        ϕ = TableCompressiblePoreVolume(tab, regions = regions)
         tab_perm = map(x -> get_1d_interpolator(x[:, 1], x[:, 3]), rt)
-        vars[:PermeabilityMultiplier] = ScalarPressureTable(tab_perm, regions = regions)
+        rockcomp = get(runspec, "ROCKCOMP", ["REVERS", 1, "NO", "CZ", 0.0])
+        if rockcomp[1] == "REVERS"
+            ϕ = TableCompressiblePoreVolume(tab, regions = regions)
+            Kfn = ScalarPressureTable(tab_perm, regions = regions)
+        else
+            if rockcomp[1] != "IRREVERS"
+                jutul_message("ROCKCOMP", "Only IRREVERS and REVERS are supported, using IRREVERS fallback for $(rockcomp[1])")
+            end
+            ϕ = HystereticTableCompressiblePoreVolume(tab, regions = regions)
+            Kfn = HystereticScalarPressureTable(tab_perm, regions = regions)
+            param[:MaxPressure] = MaxPressure()
+        end
+        vars[:PermeabilityMultiplier] = Kfn
     elseif haskey(props, "ROCK")
         rock = props["ROCK"]
         if rock isa Matrix
@@ -933,7 +944,7 @@ function set_deck_pvmult!(vars, param, sys, props, reservoir)
         static = param[:FluidVolume]
         delete!(param, :FluidVolume)
         param[:StaticFluidVolume] = static
-        vars[:FluidVolume] = wrap_reservoir_variable(sys, ϕ, :flow)
+        vars[:FluidVolume] = ϕ
     end
 end