Merge pull request #56 from sintefmath/dev

Bugfix to initialization and fault plots

Project.toml

@@ -1,7 +1,7 @@
 name = "JutulDarcy"
 uuid = "82210473-ab04-4dce-b31b-11573c4f8e0a"
 authors = ["Olav Møyner <olav.moyner@gmail.com>"]
-version = "0.2.29"
+version = "0.2.30"
 
 [deps]
 Artifacts = "56f22d72-fd6d-98f1-02f0-08ddc0907c33"

src/ext.jl

@@ -164,7 +164,7 @@ function plot_faults!(ax, mesh::UnstructuredMesh; kwarg...)
             if length(v) == 0
                 continue
             end
-            plot_mesh!(ax, mesh; faces = v, color = i, colorrange = (1, n), kwarg...)
+            plot_mesh!(ax, mesh; faces = v, color = i, colorrange = (1, max(n, 2)), kwarg...)
             i += 1
         end
     end

src/forces/bc.jl

@@ -165,7 +165,6 @@ function Jutul.apply_forces_to_equation!(acc, storage, model::SimulationModel{D,
 end
 
 function compute_bc_mass_fluxes(bc, state, nph)
-
     # Get reservoir properties
     p   = state.Pressure
     mu  = state.PhaseViscosities
@@ -181,7 +180,8 @@ function compute_bc_mass_fluxes(bc, state, nph)
     Δp = p[c] - bc.pressure
     q_tot = T_f*Δp
 
-    q = Vector{typeof(q_tot)}(undef, nph)
+    num_t = Base.promote_type(typeof(q_tot), eltype(kr), eltype(mu), eltype(rho))
+    q = zeros(MVector{nph, num_t})
     if q_tot > 0
         # Pressure inside is higher than outside, flow out from domain
         for ph in 1:nph
@@ -222,12 +222,10 @@ function compute_bc_mass_fluxes(bc, state, nph)
         end
     end
 
-    return q
-
+    return SVector{nph, num_t}(q)
 end
 
 function compute_bc_heat_fluxes(bc, state, nph)
-
     q = compute_bc_mass_fluxes(bc, state, nph)
     c = bc.cell
 

src/init/init.jl

@@ -265,6 +265,7 @@ function parse_state0_equil(model, datafile; normalize = :sum)
     has_oil = haskey(datafile["RUNSPEC"], "OIL")
     has_gas = haskey(datafile["RUNSPEC"], "GAS")
 
+    is_co2 = haskey(datafile["RUNSPEC"], "JUTUL_CO2BRINE")
     is_single_phase = (has_water + has_oil + has_gas) == 1
 
     sys = model.system
@@ -435,7 +436,12 @@ function parse_state0_equil(model, datafile; normalize = :sum)
                     contacts = []
                     contacts_pc = []
                 elseif nph == 2
-                    if has_oil && has_gas
+                    if is_co2
+                        contacts = (woc, )
+                        # TODO: Check sign here. Usually these models are
+                        # initialized without CO2.
+                        contacts_pc = (woc_pc, )
+                    elseif has_oil && has_gas
                         contacts = (goc, )
                         contacts_pc = (goc_pc, )
                     else

src/utils.jl

@@ -1335,6 +1335,7 @@ function partitioner_input(model, parameters; conn = :trans)
     if conn == :unit
         T = ones(Int, length(trans))
     else
+        trans = max.(trans, 1e-20)
         if conn == :trans
             T = copy(trans)
             T = length(T)*T./sum(T)