Merge pull request #85 from sintefmath/dev

Compositional improvements

docs/src/extras/faq.md

@@ -108,7 +108,11 @@ reservoir_model(model2).output_variables
  :PhaseMassMobilities
 ```
 
-Can also pass `extra_outputs = [:PhaseMobilities]`  to output specific variables.
+You can also pass `extra_outputs = [:PhaseMobilities]` as a keyword argument to `setup_reservoir_model` to make the resulting model output specific variables.
+
+### What is the unit and sign convention for well rates?
+
+Well results are given in strict SI, which means that rates are generally given in ``m^3/s``. Rates are positive for injection (mass entering the reservoir domain) and negative for production (leaving the reservoir domain).
 
 ## Plotting
 

ext/JutulDarcyPartitionedArraysExt/JutulDarcyPartitionedArraysExt.jl

@@ -24,7 +24,7 @@ module JutulDarcyPartitionedArraysExt
     function setup_reservoir_simulator_parray(
             case::JutulCase,
             backend::PArrayBackend;
-            conn = :unit,
+            conn = :logtrans,
             np = missing,
             kwarg...
         )

src/JutulDarcy.jl

@@ -1,5 +1,3 @@
-__precompile__(false)
-
 """
 $(README)
 """

src/formulations/pressure/functions.jl

@@ -2,11 +2,9 @@ function store_total_fluxes!(vT, model, state)
     sys = model.system
     nph = number_of_phases(sys)
     N = model.domain.representation.neighborship
+    @assert length(vT) == size(N, 2) "vT must have the same length (was $(length(vT))) as the number of faces $(size(N, 2))"
     μ = state.PhaseViscosities
     kr = state.RelativePermeabilities
-    function mob(ix)
-        return 
-    end
     for face in eachindex(vT)
         l = N[1, face]
         r = N[2, face]
@@ -23,4 +21,45 @@ function store_total_fluxes!(vT, model, state)
         end
         vT[face] = v
     end
+    return vT
+end
+
+function store_total_fluxes(model, state)
+    nf = number_of_faces(model.data_domain)
+    vT = zeros(nf)
+    return store_total_fluxes!(vT, model, state)
+end
+
+function store_phase_fluxes!(v_phases, model, state)
+    sys = model.system
+    nph = number_of_phases(sys)
+    N = model.domain.representation.neighborship
+
+    @assert size(v_phases, 1) == nph
+    @assert size(v_phases, 2) == size(N, 2)
+
+    μ = state.PhaseViscosities
+    kr = state.RelativePermeabilities
+    for face in axes(v_phases, 2)
+        l = N[1, face]
+        r = N[2, face]
+        tpfa = TPFA(l, r, 1)
+        upw = SPU(l, r)
+        f_t = Jutul.DefaultFlux()
+        common = flux_primitives(face, state, model, f_t, tpfa, upw)
+        for ph in 1:nph
+            mob = ix -> kr[ph, ix]/μ[ph, ix]
+            q = darcy_phase_kgrad_potential(face, ph, state, model, f_t, tpfa, upw, common)
+            λ_f = upwind(upw, mob, q)
+            v_phases[ph, face] = λ_f * q
+        end
+    end
+    return v_phases
+end
+
+function store_phase_fluxes(model, state)
+    nf = number_of_faces(model.data_domain)
+    nph = number_of_phases(model.system)
+    v = zeros(nph, nf)
+    return store_phase_fluxes!(v, model, state)
 end

src/input_simulation/data_input.jl

@@ -1904,7 +1904,7 @@ function select_injector_mixture_spec(sys::CompositionalSystem, name, streams, t
     mix = zeros(Float64, ncomp)
     if uppercase(type) == "WATER" || uppercase(type) == "WAT"
         has_water || throw(ArgumentError("Cannot have WATER injector without water phase."))
-        mix[1] = 1.0
+        mix[end] = 1.0
         rho = rho_s[1]
         phases_mix = ((1, 1.0), (2, 0.0), (3, 0.0))
     else
@@ -1925,7 +1925,7 @@ function select_injector_mixture_spec(sys::CompositionalSystem, name, streams, t
         )
         z_mass /= sum(z_mass)
         for i in eachindex(z_mass)
-            mix[i+offset] = z_mass[i]
+            mix[i] = z_mass[i]
         end
         @assert sum(mix) ≈ 1.0 "Sum of mixture was $(sum(mix)) != 1 for mole mixture $(z) as mass $z_mass"
 

src/multicomponent/flux.jl

@@ -61,8 +61,8 @@ function add_diffusive_component_flux(q, S, ρ, X, Y, D, face, grad, lv)
 
     @inbounds for i in eachindex(q)
         q_i = q[i]
-        dX = gradient(X, l, grad)
-        dY = gradient(Y, l, grad)
+        dX = gradient(X, i, grad)
+        dY = gradient(Y, i, grad)
 
         q_i += diff_mass_l*dX + diff_mass_v*dY
         q = setindex(q, q_i, i)

src/multicomponent/multicomponent.jl

@@ -192,8 +192,8 @@ function compositional_criterion(state, dt, active, r, nc, w, sl, liquid_density
         valw = dt*r[end, ix]/(water_density[i]*vol[i])
         if valw > e[end]
             e[end] = abs(valw)
-            mb[end] += valw
         end
+        mb[end] += valw
     end
     @. mb = dt*abs(mb)/total_mass
     return (e, mb)

src/multicomponent/variables/flash.jl

@@ -233,7 +233,7 @@ function update_flash_result(S, m, eos, phase_state, K, cond_prev, stability, x,
 
     p_val = value(P)
     T_val = value(T)
-    @. z = max(value(Z), 1e-8)
+    @. z = max(value(Z), MultiComponentFlash.MINIMUM_COMPOSITION)
 
     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)

src/multicomponent/variables/saturations.jl

@@ -16,16 +16,19 @@ end
     @inbounds for i in ix
         S_other = ImmiscibleSaturation[i]
         fr = FlashResults[i]
-        rem = one(T) - S_other + 0*MINIMUM_COMPOSITIONAL_SATURATION
+        S_eos = one(T) - S_other
         if fr.state == MultiComponentFlash.two_phase_lv
-            S_v_pure = MultiComponentFlash.two_phase_vapor_saturation(eos, Pressure[i], Temperature[i], fr)
-            S_v = rem*S_v_pure
+            S_l_pure, S_v_pure = phase_saturations(eos, Pressure[i], Temperature[i], fr)
         elseif fr.state == MultiComponentFlash.single_phase_v
-            S_v = rem
+            S_l_pure = zero(T)
+            S_v_pure = one(T)
         else
-            S_v = zero(T)
+            S_l_pure = one(T)
+            S_v_pure = zero(T)
         end
-        S_l = rem - S_v
+        S_v = S_eos*S_v_pure
+        S_l = S_eos*S_l_pure
+
         Sat[l, i] = S_l
         Sat[v, i] = S_v
         Sat[a, i] = S_other

src/multiphase.jl

@@ -510,8 +510,17 @@ function capillary_pressure(model, s)
     return(pc, ref_index)
 end
 
-get_reference_phase_index(::SinglePhaseSystem) = 1
+function get_reference_phase_index(::SinglePhaseSystem)
+    return 1
+end
+
+"""
+    get_reference_phase_index(system::JutulSystem)
 
+Get the index of the reference phase in the system. The reference phase is the
+pressure for which the pressure is given in the system. For single-phase systems
+and models without capillary pressure this is unambigious.
+"""
 function get_reference_phase_index(system::JutulSystem)
     mphases = get_phases(system)
     return get_reference_phase_index(mphases)
@@ -521,11 +530,18 @@ function get_reference_phase_index(sys::MultiPhaseSystem)
     return sys.reference_phase_index
 end
 
+"""
+    get_reference_phase_index(mphases)
+
+Get the index of the reference phase for a set of phases. The reference phase is
+selected as the liquid phase if present, otherwise the aqueous phase. If neither
+is present the first phase is selected.
+"""
 function get_reference_phase_index(mphases)
     function find_phase(k)
         ix = 0
         for (i, ph) in enumerate(mphases)
-            if ph isa LiquidPhase
+            if ph isa k
                 ix = i
                 break
             end