Tweak defaults (CPRW, CSP11C) and make compositional fluxes use poten…

src/CSP11/setup.jl

@@ -47,7 +47,7 @@ function setup_spe11_case(domain_or_dims;
         nstep_injection2 = 50,
         nstep_migration = 100,
         use_reporting_steps = false,
-        divide_c_wells = true,
+        divide_c_wells = false,
         kgrad = :tpfa,
         pc_limit = 20*si_unit(:bar),
         kwarg...

src/flux.jl

@@ -14,7 +14,6 @@ end
 end
 
 @inline function component_mass_fluxes!(q, face, state, model, flux_type, kgrad, upw)
-    T = eltype(q)
     phase_fluxes = darcy_phase_mass_fluxes(face, state, model, flux_type, kgrad, upw)
     for ph in eachindex(phase_fluxes)
         q_ph = phase_fluxes[ph]
@@ -25,6 +24,10 @@ end
 
 @inline function darcy_phase_mass_fluxes(face, state, model, flux_type, kgrad, upw, phases = eachphase(model.system))
     dpot = darcy_permeability_potential_differences(face, state, model, flux_type, kgrad, upw, phases)
+    return darcy_phase_mass_fluxes(face, state, model, flux_type, kgrad, upw, phases, dpot)
+end
+
+@inline function darcy_phase_mass_fluxes(face, state, model, flux_type, kgrad, upw, phases, dpot)
     F(phase) = darcy_phase_mass_flux(face, phase, state, model, flux_type, kgrad, upw, dpot[phase])
     return map(
         F,
@@ -50,6 +53,10 @@ end
 
 @inline function darcy_phase_volume_fluxes(face, state, model, flux_type, kgrad, upw, phases = eachphase(model.system))
     dpot = darcy_permeability_potential_differences(face, state, model, flux_type, kgrad, upw, phases)
+    return darcy_phase_volume_fluxes(face, state, model, flux_type, kgrad, upw, phases, dpot)
+end
+
+@inline function darcy_phase_volume_fluxes(face, state, model, flux_type, kgrad, upw, phases, dpot)
     F(phase) = darcy_phase_volume_flux(face, phase, state, model, flux_type, kgrad, upw, dpot[phase])
     return map(
         F,
@@ -148,7 +155,7 @@ end
     if s_l + s_r > ϵ
         ρ_avg =  (s_l*ρ_l + s_r*ρ_r)/(s_l + s_r)
     else
-        ρ_avg = 0.5*(ρ_r + ρ_l)
+        ρ_avg = (ρ_r + ρ_l)/2
     end
     return ρ_avg
 end

src/multicomponent/flux.jl

@@ -6,33 +6,35 @@
 
     X = state.LiquidMassFractions
     Y = state.VaporMassFractions
-    mass_fluxes = darcy_phase_mass_fluxes(face, state, model, flux_type, kgrad, upw)
-    q = compositional_fluxes!(q, X, Y, mass_fluxes, upw, aqua, ph_ix, component_count)
+    phases = eachphase(model.system)
+    dpot = darcy_permeability_potential_differences(face, state, model, flux_type, kgrad, upw, phases)
+    mass_fluxes = darcy_phase_mass_fluxes(face, state, model, flux_type, kgrad, upw, phases, dpot)
+    q = compositional_fluxes!(q, X, Y, mass_fluxes, dpot, upw, aqua, ph_ix, component_count)
 
     Dl = phase_diffusivities(state, LiquidPhase())
     Dv = phase_diffusivities(state, VaporPhase())
     return add_diffusive_component_flux(q, Dl, Dv, face, state, model, kgrad, component_count)
 end
 
-@inline function compositional_fluxes!(q, X, Y, mass_fluxes, upw, ::Val{false},
+@inline function compositional_fluxes!(q, X, Y, mass_fluxes, dpot, upw, ::Val{false},
         phase_ix, component_count)
     l, v = phase_ix
-    return inner_compositional!(q, X, Y, mass_fluxes[l], mass_fluxes[v], upw, component_count)
+    return inner_compositional!(q, X, Y, mass_fluxes[l], mass_fluxes[v], dpot[l], dpot[v], upw, component_count)
 end
 
-@inline function compositional_fluxes!(q, X, Y, mass_fluxes, upw, ::Val{true},
+@inline function compositional_fluxes!(q, X, Y, mass_fluxes, dpot, upw, ::Val{true},
         phase_ix, component_count::Val{N}) where N
     a, l, v = phase_ix
-    q = inner_compositional!(q, X, Y, mass_fluxes[l], mass_fluxes[v], upw, component_count)
+    q = inner_compositional!(q, X, Y, mass_fluxes[l], mass_fluxes[v], dpot[l], dpot[v], upw, component_count)
     return setindex(q, mass_fluxes[a], N + 1)
 end
 
-@inline function inner_compositional!(q::StaticVector, X, Y, q_l, q_v, upw, ::Val{N}) where N
+@inline function inner_compositional!(q::StaticVector, X, Y, q_l, q_v, dpot_l, dpot_v, upw, ::Val{N}) where N
     # Static component indices avoid SVector copies and large GPU local-memory buffers.
     return map(q, eachindex(q)) do q_i, i
         if i <= N
-            X_f = upwind(upw, cell -> @inbounds(X[i, cell]), q_l)
-            Y_f = upwind(upw, cell -> @inbounds(Y[i, cell]), q_v)
+            X_f = upwind(upw, cell -> @inbounds(X[i, cell]), dpot_l)
+            Y_f = upwind(upw, cell -> @inbounds(Y[i, cell]), dpot_v)
             oftype(q_i, q_l*X_f + q_v*Y_f)
         else
             q_i

src/utils.jl

@@ -1173,14 +1173,15 @@ these settings. Later TUNING changes in the schedule are ignored.
 
 ## Linear solver options
 
-- `linear_solver=:bicgstab`: iterative solver to use (provided model supports
-  it). Typical options are `:bicgstab` or `:gmres` Can alternatively pass a
-  linear solver instance.
-- `precond=:cpr`: preconditioner for iterative solver. For larger problems, CPR
-  variants are recommended. In order of strength and cost:
-   - `:cpr` standard Constrained-Pressure-Residual with ILU(0) second stage
+- `linear_solver=:gmres`: iterative solver to use (provided model supports it).
+  Typical options are `:bicgstab`, `:fgmres` or `:gmres` Can alternatively pass
+  a linear solver instance.
+- `precond=:cprw`: preconditioner for iterative solver. For larger problems, CPR
+  variants are recommended. `:cprw` is the default for most configurations. In
+  order of strength and cost:
+     - `:cprw` CPRW with ILU(0) second stage. Faster for problems with wells
      (strong preconditioner)
-   - `:cprw` CPRW with ILU(0) second stage. Faster for problems with wells
+   - `:cpr` standard Constrained-Pressure-Residual with ILU(0) second stage
      (strong preconditioner)
    - `:ilu0` block-incomplete-LU (intermediate strength preconditioner)
    - `:spai0`: Sparse Approximate Inverse of lowest order (weak preconditioner)
@@ -1315,15 +1316,11 @@ function setup_reservoir_simulator(case::JutulCase;
         mode = :default,
         ka_backend = missing,
         method = :newton,
-        precond = :cpr,
-        linear_solver = :bicgstab,
+        precond = missing,
+        linear_solver = ifelse(mode == :parray, :bicgstab, :gmres),
         linear_solver_backend = :auto,
         max_dt = missing,
-        max_timestep = if ismissing(max_dt)
-            reservoir_simulation_input_defaults(case, :max_timestep)
-        else
-            max_dt
-        end,
+        max_timestep = ifelse(ismissing(max_dt), reservoir_simulation_input_defaults(case, :max_timestep), max_dt),
         min_timestep = reservoir_simulation_input_defaults(case, :min_timestep),
         rtol = missing,
         initial_dt = reservoir_simulation_input_defaults(case, :initial_dt),
@@ -1374,6 +1371,24 @@ function setup_reservoir_simulator(case::JutulCase;
         reduce_memory = true,
         kwarg...
     )
+    if ismissing(precond)
+        rmodel = reservoir_model(case)
+        is_csr = rmodel.context isa ParallelCSRContext
+        is_single_node = mode != :mpi
+        is_newton = method == :newton
+        is_multimodel = case.model isa MultiModel
+        is_tracer = haskey(rmodel.equations, :tracers)
+        if is_multimodel
+            has_groups = case.model.groups isa AbstractVector && maximum(case.model.groups) > 1
+        else
+            has_groups = false
+        end
+        if is_multimodel && is_newton && is_csr && is_single_node && has_groups && !is_tracer
+            precond = :cprw
+        else
+            precond = :cpr
+        end
+    end
     warn_unsupported_tuning(case.input_data)
     ka_mode = mode isa Symbol && (mode == :ka || startswith(String(mode), "ka_"))
     sim_type_overridden = float_type != Float64 || index_type != Int

test/compositional_flux.jl

@@ -21,9 +21,9 @@ using Jutul, JutulDarcy, StaticArrays, Test
         q = zeros(SVector{n, Float64})
         q_traced = zeros(MVector{n, typeof(first(traced_fluxes))})
         numeric = JutulDarcy.compositional_fluxes!(q, liquid, vapor,
-            mass_fluxes, upw, Val(aqueous), phase_indices, Val(3))
+            mass_fluxes, mass_fluxes, upw, Val(aqueous), phase_indices, Val(3))
         traced = JutulDarcy.compositional_fluxes!(q_traced,
-            liquid_traced, vapor_traced, traced_fluxes, upw,
+            liquid_traced, vapor_traced, traced_fluxes, mass_fluxes, upw,
             Val(aqueous), phase_indices, Val(3))
 
         expected = 2.0*liquid[:, 1] - 3.0*vapor[:, 2]

test/dfm.jl

@@ -516,7 +516,7 @@ end
     @testset "immiscible simulation runs" begin
         setup = setup_fractured_cube_test(fluid_model = :two_phase)
 
-        sim_args = (info_level = -1, max_nonlinear_iterations = 15, max_timestep_cuts = 25)
+        sim_args = (info_level = -1, max_nonlinear_iterations = 15, max_timestep_cuts = 25, linear_solver = :bicgstab)
         res_fd  = simulate_reservoir(setup.case_fd;  sim_args...)
         res_dfm = simulate_reservoir(setup.case_dfm; sim_args...)
 
@@ -560,7 +560,7 @@ end
         # differences in adaptive time-stepping between the two models do not
         # dominate the temperature comparison.
         sim_args = (info_level = -1, initial_dt = 5.0,
-            tol_cnv = 1e-5, tol_cnve = 1e-6,
+            tol_cnv = 1e-5, tol_cnve = 1e-6, linear_solver = :bicgstab,
             max_timestep = setup.case_fd.dt[1]/4)
         _, cfg = setup_reservoir_simulator(setup.case_dfm; sim_args...)
         @test cfg[:tolerances][:Reservoir][:energy_conservation].CNV == sim_args.tol_cnve