diff --git a/Project.toml b/Project.toml index f7c21c41..82916974 100644 --- a/Project.toml +++ b/Project.toml @@ -56,7 +56,7 @@ Krylov = "0.9.1" LinearAlgebra = "1" LoopVectorization = "0.12.115" MAT = "0.10.3" -MultiComponentFlash = "1.1.10" +MultiComponentFlash = "1.1.13" OrderedCollections = "1.6.2" PartitionedArrays = "0.3.2" Polyester = "0.6.11, 0.7.3" diff --git a/src/CO2Properties/setup.jl b/src/CO2Properties/setup.jl index 06ed7654..1e9a881b 100644 --- a/src/CO2Properties/setup.jl +++ b/src/CO2Properties/setup.jl @@ -1,10 +1,10 @@ -function setup_reservoir_model_co2_brine(reservoir::DataDomain, - T = missing; - thermal = false, - composite = thermal, - kwarg... +function setup_reservoir_model_co2_brine(reservoir::DataDomain; + temperature = missing, + thermal = false, + composite = thermal, + kwarg... ) - tables = co2_brine_property_tables(T) + tables = co2_brine_property_tables(temperature) rho = JutulDarcy.BrineCO2MixingDensities(tables[:density]) mu = JutulDarcy.PTViscosities(tables[:viscosity]) if thermal diff --git a/src/NLDD/simulator.jl b/src/NLDD/simulator.jl index babb0a02..b1d767a1 100644 --- a/src/NLDD/simulator.jl +++ b/src/NLDD/simulator.jl @@ -335,7 +335,11 @@ function get_solve_status(rep, ok) if length(rep) > 1 status = local_solved_in_multiple_steps else - if length(only(rep)[:steps]) == 0 + rep = only(rep) + nstep = length(rep[:steps]) + if nstep == 0 + status = local_already_converged + elseif nstep == 1 && !haskey(rep[:steps][1], :linear_solver) status = local_already_converged else status = local_solved_in_single_step @@ -673,6 +677,10 @@ function solve_subdomain(sim, i, dt, forces, cfg) ok, report = Jutul.solve_timestep!(sim, dt, f, max_iter, cfg, finalize = true, prepare = true, update_explicit = false) catch e + if e isa InterruptException + # Don't leave the user trapped... + rethrow(e) + end @warn "Subdomain $i threw exception: $e" end return (ok, report) diff --git a/src/NLDD/utils.jl b/src/NLDD/utils.jl index 6277449f..e21aa13b 100644 --- a/src/NLDD/utils.jl +++ b/src/NLDD/utils.jl @@ -5,6 +5,7 @@ import JutulDarcy: set_default_cnv_mb! function simulator_config(sim::NLDDSimulator; method = :nldd, + subdomain_info_level = -1, inner_tol_mul = 1.0, inner_tol_final = 10.0, inner_max_timestep_cuts = 2, @@ -19,6 +20,9 @@ function simulator_config(sim::NLDDSimulator; inc_tol_dz = Inf, kwarg... ) + if subdomain_info_level isa Real + subdomain_info_level = i -> i == subdomain_info_level + end check_before_solve = isinf(inc_tol_dp_abs) && isinf(inc_tol_dp_rel) && isinf(inc_tol_dz) cfg = simulator_config(sim.simulator) set_default_cnv_mb!(cfg, sim.simulator.model, @@ -68,7 +72,7 @@ function simulator_config(sim::NLDDSimulator; min_nonlinear_iterations = inner_min_nonlinear_iterations, tol_factor_final_iteration = inner_tol_final, check_before_solve = check_before_solve, - info_level = -1 + info_level = Int(subdomain_info_level(i)) ) # Small hack for reservoir stuff set_default_cnv_mb!(subconfigs[i], subsims[i].model, tol_cnv = inner_tol_mul*tol_cnv, diff --git a/src/cpr.jl b/src/cpr.jl index e9b0240d..67657344 100644 --- a/src/cpr.jl +++ b/src/cpr.jl @@ -367,7 +367,7 @@ function true_impes!(w, acc, r, n, bz, arg...) elseif bz == 8 true_impes_8!(w, acc, r, n, bz, arg...) else - true_impes_gen!(w, acc, r, n, bz, arg...) + true_impes_gen!(w, acc, r, n, Val(bz), arg...) end end @@ -447,9 +447,10 @@ function true_impes_8!(w, acc, r, n, bz, s, scaling) end end -function true_impes_gen!(w, acc, r, n, bz, p_scale, scaling) +function true_impes_gen!(w, acc, r, n, ::Val{bz}, p_scale, scaling) where bz r_p = SVector{bz}(r) - A = MMatrix{bz, bz, eltype(r)}(zeros(bz, bz)) + r_T = eltype(r) + A = MMatrix{bz, bz, r_T}(undef) for cell in 1:n @inbounds for i = 1:bz v = acc[i, cell] @@ -458,7 +459,7 @@ function true_impes_gen!(w, acc, r, n, bz, p_scale, scaling) A[j, i] = v.partials[j] end end - invert_w!(w, A, r_p, cell, bz, scaling) + invert_w!(w, SMatrix{bz, bz, r_T}(A), r_p, cell, bz, scaling) end end diff --git a/src/facility/well_presolve.jl b/src/facility/well_presolve.jl index 48d11b43..c0a88bd8 100644 --- a/src/facility/well_presolve.jl +++ b/src/facility/well_presolve.jl @@ -38,7 +38,7 @@ end function Jutul.simulator_config!(cfg, sim, ::PrepareStepWellSolver, storage) add_option!(cfg, :well_iterations, 25, "Well iterations to be performed before step.", types = Int, values = 0:10000) - add_option!(cfg, :well_outer_iterations_limit, Inf, "Number of outer iterations for each solve where wells are to be solved.", types = Union{Int, Float64}) + add_option!(cfg, :well_outer_iterations_limit, 1, "Number of outer iterations for each solve where wells are to be solved.", types = Union{Int, Float64}) add_option!(cfg, :well_acceptance_factor, 10.0, "Accept well pre-solve results at this relaxed factor.", types = Float64) add_option!(cfg, :well_info_level, -1, "Info level for well solver.", types = Int) end diff --git a/src/flux.jl b/src/flux.jl index f6749c62..ef95da40 100644 --- a/src/flux.jl +++ b/src/flux.jl @@ -89,7 +89,8 @@ function face_average_density(model::Union{CompositionalModel, ThermalCompositio elseif s_r_tiny ρ_avg = ρ_l else - ρ_avg = (s_l*ρ_r + s_r*ρ_l)/(s_l + s_r) + # alt def: (s_l*ρ_r + s_r*ρ_l)/(s_l + s_r) + ρ_avg = 0.5*(ρ_r + ρ_l) end end return ρ_avg diff --git a/src/init/init.jl b/src/init/init.jl index 319ad265..049c2dcc 100644 --- a/src/init/init.jl +++ b/src/init/init.jl @@ -5,6 +5,7 @@ function equilibriate_state(model, contacts, cells = missing, rs = missing, rv = missing, + composition = missing, kwarg... ) model = reservoir_model(model) @@ -25,7 +26,7 @@ function equilibriate_state(model, contacts, init = Dict{Symbol, Any}() init = equilibriate_state!(init, pts, model, sys, contacts, datum_depth, datum_pressure; - cells = cells, rv = rv, rs = rs, kwarg... + cells = cells, rv = rv, rs = rs, composition = composition, kwarg... ) is_blackoil = sys isa StandardBlackOilSystem @@ -51,6 +52,21 @@ function equilibriate_state(model, contacts, init[:Rv] = Rv end end + is_compositional = sys isa CompositionalSystem + if is_compositional + @assert !ismissing(composition) + has_wat = has_other_phase(sys) + nc = length(init[:Pressure]) + ncomp = number_of_components(sys) - has_wat + zmf = zeros(ncomp, nc) + for i in 1:nc + zmf[:, i] = composition(pts[i]) + end + init[:OverallMoleFractions] = zmf + if has_wat + init[:ImmiscibleSaturation] = init[:Saturations][1, :] + end + end return init end @@ -58,6 +74,7 @@ function equilibriate_state!(init, depths, model, sys, contacts, depth, datum_pr cells = 1:length(depths), rs = missing, rv = missing, + composition = missing, T_z = missing, s_min = missing, contacts_pc = missing, @@ -95,24 +112,43 @@ function equilibriate_state!(init, depths, model, sys, contacts, depth, datum_pr relperm = last(relperm) end - pvt = rho.pvt reg = Int[pvtnum] function density_f(p, z, ph) - pvt_i = pvt[ph] - if phases[ph] == LiquidPhase() && disgas - rs_max = table_by_region(sys.rs_max, pvtnum) - Rs = min(rs(z), rs_max(p)) - b = JutulDarcy.shrinkage(pvt_i, reg, p, Rs, 1) - rho = b*(rhoOS + Rs*rhoGS) - elseif phases[ph] == VaporPhase() && vapoil - rv_max = table_by_region(sys.rv_max, pvtnum) - Rv = min(rv(z), rv_max(p)) - b = JutulDarcy.shrinkage(pvt_i, reg, p, Rv, 1) - rho = b*(rhoGS + Rv*rhoOS) + if rho isa ThreePhaseCompositionalDensitiesLV || rho isa TwoPhaseCompositionalDensities + if phases[ph] == AqueousPhase() + phase_density = rho_s[ph]*JutulDarcy.shrinkage(rho.immiscible_pvt, p) + else + @assert !ismissing(composition) "Composition must be present for equilibrium." + @assert !ismissing(T_z) "Temperature must be present for equilibrium." + z_i = Vector{Float64}(composition(z)) + T = T_z(z) + eos = model.system.equation_of_state + f = MultiComponentFlash.flashed_mixture_2ph(eos, (p = p, T = T, z = z_i)) + rho_l, rho_v = mass_densities(eos, p, T, f) + if phases[ph] == VaporPhase() || rho_l == 0 + phase_density = rho_v + else + phase_density = rho_l + end + end else - rho = rho_s[ph]*JutulDarcy.shrinkage(pvt_i, reg, p, 1) + pvt = rho.pvt + pvt_i = pvt[ph] + if phases[ph] == LiquidPhase() && disgas + rs_max = table_by_region(sys.rs_max, pvtnum) + Rs = min(rs(z), rs_max(p)) + b = JutulDarcy.shrinkage(pvt_i, reg, p, Rs, 1) + phase_density = b*(rhoOS + Rs*rhoGS) + elseif phases[ph] == VaporPhase() && vapoil + rv_max = table_by_region(sys.rv_max, pvtnum) + Rv = min(rv(z), rv_max(p)) + b = JutulDarcy.shrinkage(pvt_i, reg, p, Rv, 1) + phase_density = b*(rhoGS + Rv*rhoOS) + else + phase_density = rho_s[ph]*JutulDarcy.shrinkage(pvt_i, reg, p, 1) + end end - return rho + return phase_density end pressures = determine_hydrostatic_pressures(depths, depth, zmin, zmax, contacts, datum_pressure, density_f, contacts_pc) if nph > 1 @@ -207,6 +243,8 @@ function parse_state0_equil(model, datafile) init = Dict{Symbol, Any}() sol = datafile["SOLUTION"] + props = datafile["PROPS"] + G = physical_representation(model.data_domain) nc = number_of_cells(G) nph = number_of_phases(model.system) @@ -220,12 +258,17 @@ function parse_state0_equil(model, datafile) npvt = GeoEnergyIO.InputParser.number_of_tables(datafile, :pvtnum) nsat = GeoEnergyIO.InputParser.number_of_tables(datafile, :satnum) - if haskey(sol, "RTEMP") - Ti = convert_to_si(only(sol["RTEMP"]), :Celsius) + if haskey(sol, "RTEMP") || haskey(props, "RTEMP") + if haskey(props, "RTEMP") + rtmp = only(props["RTEMP"]) + else + rtmp = only(sol["RTEMP"]) + end + Ti = convert_to_si(rtmp, :Celsius) T_z = z -> Ti - elseif haskey(sol, "TEMPVD") - z = vec(sol["TEMPVD"][:, 1]) - Tvd = vec(sol["TEMPVD"][:, 2] + 273.15) + elseif haskey(props, "TEMPVD") + z = vec(props["TEMPVD"][:, 1]) + Tvd = vec(props["TEMPVD"][:, 2] + 273.15) T_z = get_1d_interpolator(z, Tvd) else T_z = missing @@ -359,9 +402,34 @@ function parse_state0_equil(model, datafile) else rv = missing end + if haskey(props, "ZMFVD") + ztab = props["ZMFVD"][ereg] + N = size(ztab, 2) - 1 + Sz = SVector{N, Float64} + ztab_z = Float64[] + ztab_static = Sz[] + for row in size(ztab, 1) + mf_i = Sz(ztab[row, 2:end]) + z_i = ztab[row, 1] + push!(ztab_z, z_i) + push!(ztab_static, mf_i) + if row == 1 + push!(ztab_z, z_i-1.0) + push!(ztab_static, mf_i) + end + if row == size(ztab, 1) + push!(ztab_z, z_i+1.0) + push!(ztab_static, mf_i) + end + end + composition = get_1d_interpolator(ztab_z, ztab_static) + else + composition = missing + end subinit = equilibriate_state( model, contacts, datum_depth, datum_pressure, + composition = composition, cells = cells, pvtnum = preg, contacts_pc = contacts_pc, diff --git a/src/input_simulation/data_input.jl b/src/input_simulation/data_input.jl index d2dcf14c..93fa218c 100644 --- a/src/input_simulation/data_input.jl +++ b/src/input_simulation/data_input.jl @@ -426,8 +426,8 @@ function filter_inactive_completions!(completions_vector, g) if !haskey(tuple_active, k) ix = cell_index(g, tupl, throw = false) if isnothing(ix) - jutul_message("Removing well", - "Well $well completion $tupl was declared using COMPDAT, but that cell is not active in model. Skipped.", + jutul_message("$well completion", + "Removed COMPDAT as $tupl is not active in processed mesh.", color = :yellow ) tuple_active[k] = false @@ -824,7 +824,7 @@ function parse_physics_types(datafile; pvt_region = missing) push!(phases, VaporPhase()) push!(rhoS, rhoGS) - cnames = props["CNAMES"] + cnames = copy(props["CNAMES"]) acf = props["ACF"] mw = props["MW"] p_c = props["PCRIT"] @@ -837,6 +837,7 @@ function parse_physics_types(datafile; pvt_region = missing) A_ij = nothing end mp = MolecularProperty.(mw, p_c, T_c, V_c, acf) + @assert length(cnames) == length(mp) mixture = MultiComponentMixture(mp, A_ij = A_ij, names = cnames) if haskey(props, "EOS") eos_str = uppercase(props["EOS"]) @@ -853,7 +854,12 @@ function parse_physics_types(datafile; pvt_region = missing) else eos_type = PengRobinson() end - eos = GenericCubicEOS(mixture, eos_type) + if haskey(props, "SSHIFT") + vshift = Tuple(props["SSHIFT"]) + else + vshift = nothing + end + eos = GenericCubicEOS(mixture, eos_type, volume_shift = vshift) sys = MultiPhaseCompositionalSystemLV(eos, phases, reference_densities = rhoS) else if has_oil diff --git a/src/multicomponent/variables/flash.jl b/src/multicomponent/variables/flash.jl index 7699f827..2af92804 100644 --- a/src/multicomponent/variables/flash.jl +++ b/src/multicomponent/variables/flash.jl @@ -7,12 +7,24 @@ mutable struct InPlaceFlashBuffer end end -struct FlashResults <: ScalarVariable - storage - method - update_buffer +@enum PERFORMED_FLASH_TYPE FLASH_SINGLE_PHASE_TESTED FLASH_SINGLE_PHASE_BYPASSED FLASH_FULL FLASH_RESTARTED FLASH_PURE_WATER +struct FlashResults{F_t, S_t, B_t} <: ScalarVariable + storage::S_t + method::F_t + update_buffer::B_t use_threads::Bool - function FlashResults(system; method = SSIFlash(), threads = Threads.nthreads() > 1) + tolerance::Float64 + tolerance_bypass::Float64 + stability_bypass::Bool + reuse_guess::Bool + function FlashResults(system; + method::MultiComponentFlash.AbstractFlash = SSIFlash(), + threads = Threads.nthreads() > 1, + tolerance = 1e-8, + tolerance_bypass = 10, + reuse_guess = false, + stability_bypass = false + ) eos = system.equation_of_state nc = MultiComponentFlash.number_of_components(eos) if has_other_phase(system) @@ -28,14 +40,30 @@ struct FlashResults <: ScalarVariable N = 1 end for i = 1:N - s = flash_storage(eos, method = method, inc_jac = true, diff_externals = true, npartials = n, static_size = true) + s = flash_storage(eos, + method = method, + inc_jac = true, + inc_bypass = stability_bypass, + diff_externals = true, + npartials = n, + static_size = true + ) push!(storage, s) b = InPlaceFlashBuffer(nc) push!(buffers, b) end storage = tuple(storage...) buffers = tuple(buffers...) - new(storage, method, buffers, threads) + new{typeof(method), typeof(storage), typeof(buffers)}( + storage, + method, + buffers, + threads, + tolerance, + tolerance_bypass, + stability_bypass, + reuse_guess + ) end end @@ -91,8 +119,28 @@ function flash_entity_loop!(flash_results, fr, model, eos, Pressure, Temperature update_flash_buffer!(buf, eos, Pressure, Temperature, OverallMoleFractions) buf_z = buf.z buf_forces = buf.forces + flash_entity_loop_impl!(flash_results, ix, S, fr, m, eos, buf_z, buf_forces, Pressure, Temperature, OverallMoleFractions, sw) +end + +function flash_entity_loop_impl!(flash_results, ix, S, fr, m, eos, buf_z, buf_forces, Pressure, Temperature, OverallMoleFractions, sw) + extra_print = false + if extra_print + code_log = Dict( + FLASH_SINGLE_PHASE_TESTED => 0, + FLASH_SINGLE_PHASE_BYPASSED => 0, + FLASH_FULL => 0, + FLASH_RESTARTED => 0, + FLASH_PURE_WATER => 0 + ) + end @inbounds for i in ix - flash_results[i] = internal_flash!(flash_results[i], S, m, eos, buf_z, buf_forces, Pressure, Temperature, OverallMoleFractions, sw, i) + flash_results[i], code = internal_flash!(flash_results[i], S, fr, m, eos, buf_z, buf_forces, Pressure, Temperature, OverallMoleFractions, sw, i) + if extra_print + code_log[code] += 1 + end + end + if extra_print + @info "Flash is done:" code_log end end @@ -141,7 +189,7 @@ function update_flash_buffer!(buf, eos::KValuesEOS, Pressure, Temperature, Overa return nothing end -function internal_flash!(f, S, m, eos, buf_z, buf_forces, Pressure, Temperature, OverallMoleFractions, sw, i) +function internal_flash!(f, S, var, m, eos, buf_z, buf_forces, Pressure, Temperature, OverallMoleFractions, sw, i) @inline function immiscible_sat(::Nothing, i) return 0.0 end @@ -159,32 +207,187 @@ function internal_flash!(f, S, m, eos, buf_z, buf_forces, Pressure, Temperature, K = f.K x = f.liquid.mole_fractions y = f.vapor.mole_fractions - - return update_flash_result(S, m, eos, K, x, y, buf_z, buf_forces, P, T, Z, Sw) + V = f.V + b = f.critical_distance + + return update_flash_result(S, m, eos, f.state, K, f.flash_cond, f.flash_stability, x, y, buf_z, V, buf_forces, P, T, Z, Sw, + critical_distance = b, + tolerance = var.tolerance, + stability_bypass = var.stability_bypass, + reuse_guess = var.reuse_guess, + tolerance_bypass = var.tolerance_bypass + ) end end - -function update_flash_result(S, m, eos, K, x, y, z, forces, P, T, Z, Sw = 0.0) +import MultiComponentFlash: michelsen_critical_point_measure! +function update_flash_result(S, m, eos, phase_state, K, cond_prev, stability, x, y, z, V, forces, P, T, Z, Sw = 0.0; + critical_distance::Float64 = NaN, + tolerance::Float64 = 1e-8, + tolerance_bypass::Float64 = 10.0, + stability_bypass::Bool = false, + reuse_guess::Bool = false, + ) + + # We can check for bypass if feature is enabled, we have a critical distance + # that is finite and we were in single phase previously. + # do_full_flash(c) = flash_2ph!(S, K, eos, c, NaN, method = m, extra_out = false, z_min = nothing) + + p_val = value(P) + T_val = value(T) @. z = max(value(Z), 1e-8) - # Conditions - c = (p = value(P), T = value(T), z = z) - # Perform flash - if is_pure_single_phase(Sw) - vapor_frac = NaN + + new_cond = (p = p_val, T = T_val, z = z) + do_flash, critical_distance, single_phase_code = single_phase_bypass_check(eos, new_cond, cond_prev, phase_state, Sw, critical_distance, stability_bypass, tolerance_bypass) + + if do_flash + vapor_frac, stability, return_code = two_phase_flash_implementation!(K, S, m, eos, phase_state, x, y, new_cond, V, reuse_guess) + is_single_phase = isnan(vapor_frac) + if is_single_phase && stability_bypass + # If we did a full flash and single-phase prevailed outside the + # shadow region it is time to update the critical distance for + # future reference. + if stability.liquid.trivial && stability.vapor.trivial + # Outside shadow region, flash bypass can be enabled + critical_distance = michelsen_critical_point_measure!(S.bypass, eos, new_cond.p, new_cond.T, new_cond.z) + else + # Inside shadow region, we cannot use flash bypass directly + critical_distance = NaN + end + end + # Update the condition only if we actually did a flash, otherwise we + # keep the value where we last flashed around for stability bypass + # testing. + @. cond_prev.z = z + cond_prev = (p = p_val, T = T_val, cond_prev.z) else - vapor_frac = flash_2ph!(S, K, eos, c, NaN, method = m, extra_out = false, z_min = nothing) + return_code = single_phase_code + is_single_phase = true end force_coefficients!(forces, eos, (p = P, T = T, z = Z)) - if isnan(vapor_frac) + update_condition = false + if is_single_phase # Single phase condition. Life is easy. - Z_L, Z_V, V, phase_state = single_phase_update!(P, T, Z, x, y, forces, eos, c) + Z_L, Z_V, V, phase_state = single_phase_update!(P, T, Z, x, y, forces, eos, new_cond) else # Two-phase condition: We have some work to do. - Z_L, Z_V, V, phase_state = two_phase_update!(S, P, T, Z, x, y, K, vapor_frac, forces, eos, c) + Z_L, Z_V, V, phase_state = two_phase_update!(S, P, T, Z, x, y, K, vapor_frac, forces, eos, new_cond) + # Reset critical distance to NaN and set condition to last flash since we are now in two-phase region. + critical_distance = NaN end - out = FlashedMixture2Phase(phase_state, K, V, x, y, Z_L, Z_V) - return out + out = FlashedMixture2Phase(phase_state, K, V, x, y, Z_L, Z_V, critical_distance, cond_prev, stability) + return (out, return_code) +end + +function single_phase_bypass_check(eos, new_cond, old_cond, phase_state, Sw, critical_distance, stability_bypass, ϵ) + is_pure_water = is_pure_single_phase(Sw) + was_single_phase = phase_state == MultiComponentFlash.single_phase_l || phase_state == MultiComponentFlash.single_phase_v + + if is_pure_water + # Set this to make sure that a proper stability test happens if we leave + # the pure water condition. + critical_distance = NaN + need_two_phase_flash = false + code = FLASH_PURE_WATER + elseif stability_bypass && was_single_phase && isfinite(critical_distance) + z_diff = p_diff = T_diff = 0.0 + # This is put early while we have the old z values in buffer + for i in eachindex(old_cond.z, new_cond.z) + z_diff = max(z_diff, abs(old_cond.z[i] - value(new_cond.z[i]))) + end + p_diff = abs(new_cond.p - old_cond.p) + T_diff = abs(new_cond.T - old_cond.T) + + b_old = critical_distance + z_crit = z_diff ≥ b_old/ϵ + p_crit = p_diff ≥ b_old*new_cond.p/ϵ + T_crit = T_diff ≥ b_old*ϵ + need_two_phase_flash = z_crit || p_crit || T_crit + if need_two_phase_flash + code = FLASH_SINGLE_PHASE_TESTED + else + code = FLASH_SINGLE_PHASE_BYPASSED + end + if false && need_two_phase_flash + # Hacked in expensive test to check if the stability bypass is ok. + # Can be manually activated. + V = flash_2ph(eos, new_cond) + if isnan(V) + println("Bypass OK.") + else + error("Flash bypass was wrong for conditions $new_cond giving $V") + end + end + else + need_two_phase_flash = true + code = FLASH_SINGLE_PHASE_TESTED + end + return (need_two_phase_flash, critical_distance, code) +end + +function two_phase_flash_implementation!(K, S, m, eos, old_phase_state, x, y, flash_cond, V, reuse_guess) + # Have to do some kind of flash, could be single or two-phase. + need_full_flash = true + if reuse_guess && old_phase_state == MultiComponentFlash.two_phase_lv + # If the model was previously two-phase and this option is enabled, we + # can try reusing the previous solution as an initial guess. This + # follows Rasmussen et al where a theta estimating phase partition is + # used, adapted to the K-value initial guess used by our flash. + V = value(V) + z = flash_cond.z + K = estimate_K_values_from_previous_flash!(K, V, x, y, z) + try + vapor_frac, K, stats = MultiComponentFlash.flash_2ph_impl!(S, K, eos, flash_cond, V, + method = SSINewtonFlash(swap_iter = 2), + maxiter = 20, + z_min = nothing, + check = false + ) + need_full_flash = !stats.converged || V <= 0.0 || V >= 1.0 + catch e + jutul_message("Flash", "Exception ocurred in flash: $(typeof(e)), falling back to SSI with stability test.", color = :red) + end + end + + if need_full_flash + vapor_frac, K, stats = MultiComponentFlash.flash_2ph_impl!(S, K, eos, flash_cond, NaN, + method = m, + z_min = nothing + ) + code = FLASH_FULL + else + code = FLASH_RESTARTED + end + return (vapor_frac, stats.stability, code) +end + +function estimate_K_values_from_previous_flash!(K, V, x, y, z) + ϵ = MultiComponentFlash.MINIMUM_COMPOSITION + x_t = 0.0 + y_t = 0.0 + for i in eachindex(z, x, y, K) + z_i = z[i] + x_i = value(x[i]) + y_i = value(y[i]) + + # x_i = max(x_i, ϵ) + # y_i = max(y_i, ϵ) + θ_i = V*y_i/(V*y_i + (1.0 - V)*x_i) + # Use K as working buffer + K[i] = θ_i + # Sum up molar amounts + x_t += (1.0 - θ_i)*z[i] + y_t += θ_i*z[i] + end + + for i in eachindex(K) + # Normalize estimates to mole fractions and get K-values + θ_i = K[i] + liquid = (1.0 - θ_i)/x_t + vapor = θ_i/y_t + K[i] = vapor/liquid # y_i / x_i + end + return K end function get_compressibility_factor(forces, eos, P, T, Z, phase = :unknown) diff --git a/src/multicomponent/wells.jl b/src/multicomponent/wells.jl index 4b62c284..9aa2b420 100644 --- a/src/multicomponent/wells.jl +++ b/src/multicomponent/wells.jl @@ -172,7 +172,8 @@ function separator_flash!(flash, eos, cond, z) Temperature = cond.T update_flash_buffer!(buf, eos, Pressure, Temperature, z) forces = buf.forces - return update_flash_result(fstorage, SSIFlash(), eos, f.K, x, y, buf.z, forces, Pressure, Temperature, z, 0.0) + result, code = update_flash_result(fstorage, SSIFlash(), eos, f.state, f.K, f.flash_cond, f.flash_stability, x, y, buf.z, NaN, forces, Pressure, Temperature, z, 0.0) + return result end function flash_stream!(moles::SVector{N, T}, flash, eos, cond) where {N, T} @@ -255,5 +256,5 @@ function compositional_surface_densities(state, system, S_l::S_T, S_v::S_T, rho_ volume[l] = S_l*rem volume[v] = S_v*rem @assert sum(volume) ≈ 1.0 "Volume should sum to 1, was $(sum(volume))" - return (rhoS, volume) + return (SVector(rhoS), SVector(volume)) end diff --git a/src/utils.jl b/src/utils.jl index eab5db52..762d693a 100644 --- a/src/utils.jl +++ b/src/utils.jl @@ -138,8 +138,8 @@ function reservoir_system(; flow = missing, thermal = missing, kwarg...) end """ - setup_reservoir_model(reservoir, system; wells = [], ) - setup_reservoir_model(reservoir, system; wells = [], context = DefaultContext(), reservoir_context = nothing, backend = :csc, ) + model, parameters = setup_reservoir_model(reservoir, system; wells = [], ) + model, parameters = setup_reservoir_model(reservoir, system; wells = [w1, w2], backend = :csr, ) Set up a reservoir `MultiModel` for a given reservoir `DataDomain` typically set up from [`reservoir_domain`](@ref) and an optional vector of wells that are @@ -147,6 +147,68 @@ created using [`setup_vertical_well`](@ref) and [`setup_well`](@ref). The routine automatically sets up a facility and couples the wells with the reservoir and that facility. + +# Keyword arguments + +- `wells=[]`: Vector of wells (e.g. from [`setup_well`](@ref)) that are to be + used in the model. Each well must have a unique name. +- `extra_out=true`: Return both the model and the parameters instead of just the + model. + +## Advanced model setup +- `split_wells=false`: Add a facility model for each well instead of one + facility model that controls all wells. This must be set to `true` if you want + to use MPI or nonlinear domain decomposition. +- `backend=:csc`: Backend to use. Can be `:csc` for serial compressed sparse + column CSC matrix, `:csr` for parallel compressed sparse row matrix. `:csr`` + is recommended when using MPI. +- `context=DefaultContext()`: Context used for entire model. Not recommended to + set up manually, use `backend` instead. +- `assemble_wells_together=true`: Assemble wells in a single big matrix rather + than many small matrices. +- `extra_outputs=Symbol[]`: Extra output variables for reservoir model. Defaults + to "typical" values seen in reservoir simulation. Valid values: Vector of + symbols to be output, `true` for all variables and `false` for the minimal set + required to restart simulations (typically only the primary variables) + +## Increment and variable options +These options govern the range of values and the maximum allowable change of +properties over a single Newton iteration. Changing values for maximum change +will not change the equations themselves, but the values will change the rate of +nonlinear solver convergence. Typically, smaller values are more conservative +and reduce numerical difficulties, but can significantly increase the number of +iterations and the reduce the length of the average time-step. Setting very +small values can make it infeasible to solve the problems in a reasonable time. + +Note that relative values are usually given relative to the cell value. If your +expected output values are close to zero (e.g. for near-atmospheric pressure) +low values can lead to slow convergence. + +- `dp_max_abs=nothing`: Maximum allowable pressure change in SI units (Pa) +- `dp_max_rel=0.2`: Maximum allowable relative pressure change (default is 20%) +- `dp_max_abs_well=convert_to_si(50, :bar)`: Maximum allowable pressure change + for wells in SI units (Pa) +- `dp_max_rel_well=nothing`: Maximum allowable relative pressure change in well +- `ds_max=0.2`: Maximum change in saturations +- `dz_max=`: Maximum change in composition (for compositional models only) +- `p_min=`: Minimum pressure in model (hard limit) +- `p_max=`: Maximum pressure in model (hard limit) +- `dr_max=`: Maximum change in Rs/Rv for blackoil models over a Newton + iteration. Taken relative to the saturated value of the cell. +- `dT_max_rel=`: Maximum relative change in temperature (JutulDarcy uses Kelvin, + so comments about changing limits near zero above does not apply to typical + reservoir temperatures) +- `dT_max_abs=`: Maximum absolute change in temperature (in °K/°C) +- `fast_flash=false`: Shorthand to enable `flash_reuse_guess` and + `flash_stability_bypass`. These options can together speed up the time spent + in flash solver for compositional models. Options are based on "Increasing the + Computational Speed of Flash Calculations With Applications for Compositional, + Transient Simulations" by Rasmussen et al (2006). +- `flash_reuse_guess=fast_flash`: Reuse previous flash guess when a cell remains + in two-phase. +- `flash_stability_bypass=fast_flash`: Bypass stability testing for cells + outside the two-phase and shadow region. + """ function setup_reservoir_model(reservoir::DataDomain, system::JutulSystem; wells = [], @@ -169,6 +231,9 @@ function setup_reservoir_model(reservoir::DataDomain, system::JutulSystem; dr_max = Inf, dT_max_rel = nothing, dT_max_abs = nothing, + fast_flash = false, + flash_reuse_guess = fast_flash, + flash_stability_bypass = fast_flash, parameters = Dict{Symbol, Any}(), kwarg... ) @@ -199,14 +264,25 @@ function setup_reservoir_model(reservoir::DataDomain, system::JutulSystem; ds_max = ds_max, dz_max = dz_max, dT_max_rel = dT_max_rel, - dT_max_abs = dT_max_abs + dT_max_abs = dT_max_abs, + flash_reuse_guess = flash_reuse_guess, + flash_stability_bypass = flash_stability_bypass ) - for k in extra_outputs - if haskey(rmodel.secondary_variables, k) - push!(rmodel.output_variables, k) + if extra_outputs isa Bool + if extra_outputs + for k in keys(rmodel.secondary_variables) + push!(rmodel.output_variables, k) + end + unique!(rmodel.output_variables) + end + else + for k in extra_outputs + if haskey(rmodel.secondary_variables, k) + push!(rmodel.output_variables, k) + end + unique!(rmodel.output_variables) end - unique!(rmodel.output_variables) end if haskey(reservoir, :diffusion) || haskey(reservoir, :diffusivity) rmodel.parameters[:Diffusivities] = Diffusivities() @@ -235,7 +311,9 @@ function setup_reservoir_model(reservoir::DataDomain, system::JutulSystem; ds_max = ds_max, dz_max = dz_max, dT_max_rel = dT_max_rel, - dT_max_abs = dT_max_abs + dT_max_abs = dT_max_abs, + flash_reuse_guess = flash_reuse_guess, + flash_stability_bypass = flash_stability_bypass ) models[wname] = wmodel if split_wells @@ -267,7 +345,19 @@ function setup_reservoir_model(reservoir::DataDomain, label::Symbol; kwarg...) return setup_reservoir_model(reservoir, Val(label); kwarg...) end -function set_reservoir_variable_defaults!(model; p_min, p_max, dp_max_abs, dp_max_rel, ds_max, dz_max, dr_max, dT_max_rel = nothing, dT_max_abs = nothing) +function set_reservoir_variable_defaults!(model; + p_min, + p_max, + dp_max_abs, + dp_max_rel, + ds_max, + dz_max, + dr_max, + dT_max_rel = nothing, + dT_max_abs = nothing, + flash_reuse_guess = false, + flash_stability_bypass = flash_reuse_guess + ) # Replace various variables - if they are available replace_variables!(model, OverallMoleFractions = OverallMoleFractions(dz_max = dz_max), throw = false) replace_variables!(model, Saturations = Saturations(ds_max = ds_max), throw = false) @@ -275,6 +365,12 @@ function set_reservoir_variable_defaults!(model; p_min, p_max, dp_max_abs, dp_ma replace_variables!(model, ImmiscibleSaturation = ImmiscibleSaturation(ds_max = ds_max), throw = false) replace_variables!(model, BlackOilUnknown = BlackOilUnknown(ds_max = ds_max, dr_max = dr_max), throw = false) + rmodel = reservoir_model(model) + if rmodel isa CompositionalModel + sys = rmodel.system + flash = FlashResults(sys, stability_bypass = flash_stability_bypass, reuse_guess = flash_reuse_guess) + replace_variables!(model, FlashResults = flash, throw = false) + end p_def = Pressure( max_abs = dp_max_abs, max_rel = dp_max_rel, @@ -393,8 +489,9 @@ list a few of the most relevant entries here for convenience: - `info_level = 0`: Output level. Set to 0 for minimal output, -1 for no output and 1 or more for increasing verbosity. - `output_path`: Path to write output to. -- `relaxation=Jutul.NoRelaxation`: Dampening used for solves. Can be set to - `Jutul.SimpleRelaxation()` for difficult models. +- `relaxation=Jutul.NoRelaxation()`: Dampening used for solves. Can be set to + `Jutul.SimpleRelaxation()` for difficult models. Equivialent option is to set + `true` for relaxation and `false` for no relaxation. - `failure_cuts_timestep=true`: Cut timestep instead of throwing an error when numerical issues are encountered (e.g. linear solver divergence). - `max_timestep_cuts=25`: Maximum number of timestep cuts before a solver gives @@ -427,6 +524,7 @@ function setup_reservoir_simulator(case::JutulCase; inc_tol_dz = Inf, set_linear_solver = true, timesteps = :auto, + relaxation = false, presolve_wells = false, parray_arg = Dict{Symbol, Any}(), linear_solver_arg = Dict{Symbol, Any}(), @@ -511,10 +609,18 @@ function setup_reservoir_simulator(case::JutulCase; else extra_kwarg[:linear_solver] = linear_solver end + if relaxation isa Bool + if relaxation + relaxation = SimpleRelaxation() + else + relaxation = NoRelaxation() + end + end cfg = simulator_config(sim; extra_kwarg..., timestep_selectors = sel, info_level = info_level, + relaxation = relaxation, max_timestep = max_timestep, failure_cuts_timestep = failure_cuts_timestep, max_timestep_cuts = max_timestep_cuts, diff --git a/test/multimodel.jl b/test/multimodel.jl index 9a4393c3..a916b346 100644 --- a/test/multimodel.jl +++ b/test/multimodel.jl @@ -131,6 +131,7 @@ end @testset "Block=$b, backend=$b, defaulted linsolve=$default_linsolve" begin arg = (general_ad = gen_ad, backend = backend, use_blocks = b, default_linsolve = default_linsolve) test_compositional_with_wells(; arg...) + test_compositional_with_wells(; fast_flash = true, arg...) test_immiscible_with_wells(; arg...) test_blackoil_with_wells(; arg...) end