diff --git a/Project.toml b/Project.toml index b9303555..3774d4c4 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "JutulDarcy" uuid = "82210473-ab04-4dce-b31b-11573c4f8e0a" authors = ["Olav Møyner "] -version = "0.2.33" +version = "0.2.34" [deps] AlgebraicMultigrid = "2169fc97-5a83-5252-b627-83903c6c433c" @@ -55,8 +55,8 @@ 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" @@ -64,7 +64,7 @@ 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" diff --git a/docs/src/man/basics/systems.md b/docs/src/man/basics/systems.md index a769f18b..cdfb4f3c 100644 --- a/docs/src/man/basics/systems.md +++ b/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 ) ``` diff --git a/examples/intro_sensitivities.jl b/examples/intro_sensitivities.jl index 3ee362dc..a6adf095 100644 --- a/examples/intro_sensitivities.jl +++ b/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) diff --git a/examples/two_phase_unstable_gravity.jl b/examples/two_phase_unstable_gravity.jl index 4390cd14..be9704a3 100644 --- a/examples/two_phase_unstable_gravity.jl +++ b/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) diff --git a/examples/validation_compositional.jl b/examples/validation_compositional.jl index b3e6c55d..014d5dbe 100644 --- a/examples/validation_compositional.jl +++ b/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 diff --git a/examples/validation_norne_nohyst.jl b/examples/validation_norne_nohyst.jl index c47688fe..b94624af 100644 --- a/examples/validation_norne_nohyst.jl +++ b/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 diff --git a/src/CO2Properties/generation.jl b/src/CO2Properties/generation.jl index 0fc1a6c7..771758c7 100644 --- a/src/CO2Properties/generation.jl +++ b/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) diff --git a/src/deck_support.jl b/src/deck_support.jl index b9e6992d..51fef939 100644 --- a/src/deck_support.jl +++ b/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 diff --git a/src/deck_types.jl b/src/deck_types.jl index 128a8844..5d536bb7 100644 --- a/src/deck_types.jl +++ b/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 diff --git a/src/input_simulation/mrst_input.jl b/src/input_simulation/mrst_input.jl index c5035eb6..1e0597f2 100644 --- a/src/input_simulation/mrst_input.jl +++ b/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 diff --git a/src/variables/capillary.jl b/src/variables/capillary.jl index e21c71b8..5fffa0b8 100644 --- a/src/variables/capillary.jl +++ b/src/variables/capillary.jl @@ -136,7 +136,13 @@ end cap = pc.pc npc = size(Δp, 1) scale = pc.scaling + reference_ph = get_reference_phase_index(model.system) if npc == 1 + if reference_ph == 1 + w = 2 + else + w = 1 + end pcow = only(cap) @inbounds for c in ix reg = region(pc.regions, c) @@ -145,12 +151,20 @@ end Δp[1, c] = scale[1, c]*pcow_c(sw) end elseif npc == 2 + if reference_ph == 1 + w, g = 2, 3 + elseif reference_ph == 2 + w, g = 1, 3 + else + @assert reference_ph == 3 + w, g = 1, 2 + end pcow, pcog = cap if isnothing(pcow) @inbounds for c in ix reg = region(pc.regions, c) pcog_c = table_by_region(pcog, reg) - sg = Saturations[3, c] + sg = Saturations[g, c] Δp[1, c] = 0 Δp[2, c] = scale[2, c]*pcog_c(sg) end @@ -158,8 +172,8 @@ end @inbounds for c in ix reg = region(pc.regions, c) pcow_c = table_by_region(pcow, reg) - sw = Saturations[1, c] - Δp[1, c] = -scale[1, c]*pcow_c(sw) + sw = Saturations[w, c] + Δp[1, c] = scale[1, c]*pcow_c(sw) Δp[2, c] = 0 end else @@ -167,9 +181,9 @@ end reg = region(pc.regions, c) pcow_c = table_by_region(pcow, reg) pcog_c = table_by_region(pcog, reg) - sw = Saturations[1, c] - sg = Saturations[3, c] - Δp[1, c] = -scale[1, c]*pcow_c(sw) + sw = Saturations[w, c] + sg = Saturations[g, c] + Δp[1, c] = scale[1, c]*pcow_c(sw) Δp[2, c] = scale[2, c]*pcog_c(sg) end end diff --git a/src/variables/relperm/hysteresis.jl b/src/variables/relperm/hysteresis.jl index 3a46a8a8..ff57151d 100644 --- a/src/variables/relperm/hysteresis.jl +++ b/src/variables/relperm/hysteresis.jl @@ -51,16 +51,24 @@ struct MaxSaturations <: PhaseVariables end function Jutul.update_parameter_before_step!(s_max, ::MaxSaturations, storage, model, dt, forces) s = storage.state.Saturations - update_max_saturations!(s_max, s) + update_max_hysteresis_value!(s_max, s) return s_max end -function update_max_saturations!(s_max, s) - for i in eachindex(s_max, s) - s_prev = s_max[i] - s_now = value(s[i]) - if s_now > s_prev - s_max[i] = replace_value(s_prev, s_now) +struct MaxPressure <: ScalarVariable end + +function Jutul.update_parameter_before_step!(p_max, ::MaxPressure, storage, model, dt, forces) + p = storage.state.Pressure + update_max_hysteresis_value!(p_max, p) + return p_max +end + +function update_max_hysteresis_value!(v_max, v) + for i in eachindex(v_max, v) + v_prev = v_max[i] + v_now = value(v[i]) + if v_now > v_prev + v_max[i] = replace_value(v_prev, v_now) end end end