Skip to content

Quarter-five-spot example

The quarter-five-spot is a standard test problem that simulates 1/4 of the five spot well pattern by assuming axial symmetry. The problem contains an injector in one corner and the producer in the opposing corner, with a significant volume of fluids injected into the domain.

julia
using JutulDarcy, Jutul
nx = 50;

Setup

We define a function that, for a given porosity field, computes a solution with an estimated permeability field. For assumptions and derivation of the specific form of the Kozeny-Carman relation used in this example, see Lie, Knut-Andreas. An introduction to reservoir simulation using MATLAB/GNU Octave: User guide for the MATLAB Reservoir Simulation Toolbox (MRST). Cambridge University Press, 2019, Section 2.5.2

julia
function perm_kozeny_carman(Φ)
    return ((Φ^3)*(1e-5)^2)/(0.81*72*(1-Φ)^2);
end

function simulate_qfs(porosity = 0.2)
    Dx = 1000.0
    Dz = 10.0
    Darcy = 9.869232667160130e-13
    Darcy, bar, kg, meter, Kelvin, day, sec = si_units(:darcy, :bar, :kilogram, :meter, :Kelvin, :day, :second)

    mesh = CartesianMesh((nx, nx, 1), (Dx, Dx, Dz))
    K = perm_kozeny_carman.(porosity)
    domain = reservoir_domain(mesh, permeability = K, porosity = porosity)
    Inj = setup_vertical_well(domain, 1, 1, name = :Injector);
    Prod = setup_vertical_well(domain, nx, nx, name = :Producer);
    phases = (LiquidPhase(), VaporPhase())
    rhoLS = 1000.0*kg/meter^3
    rhoGS = 700.0*kg/meter^3
    rhoS = [rhoLS, rhoGS]
    sys = ImmiscibleSystem(phases, reference_densities = rhoS)
    model, parameters = setup_reservoir_model(domain, sys, wells = [Inj, Prod])
    c = [1e-6/bar, 1e-6/bar]
    ρ = ConstantCompressibilityDensities(p_ref = 150*bar, density_ref = rhoS, compressibility = c)
    kr = BrooksCoreyRelativePermeabilities(sys, [2.0, 2.0])
    replace_variables!(model, PhaseMassDensities = ρ, RelativePermeabilities = kr);

    state0 = setup_reservoir_state(model, Pressure = 150*bar, Saturations = [1.0, 0.0])
    dt = repeat([30.0]*day, 12*10)
    dt = vcat([0.1, 1.0, 10.0], dt)
    inj_rate = Dx*Dx*Dz*0.2/sum(dt) # 1 PVI if average porosity is 0.2

    rate_target = TotalRateTarget(inj_rate)
    I_ctrl = InjectorControl(rate_target, [0.0, 1.0], density = rhoGS)
    bhp_target = BottomHolePressureTarget(50*bar)
    P_ctrl = ProducerControl(bhp_target)
    controls = Dict()
    controls[:Injector] = I_ctrl
    controls[:Producer] = P_ctrl
    forces = setup_reservoir_forces(model, control = controls)
    return simulate_reservoir(state0, model, dt, parameters = parameters, forces = forces)
end
simulate_qfs (generic function with 2 methods)

Simulate base case

This will give the solution with uniform porosity of 0.2.

julia
ws, states, report_time = simulate_qfs();
Jutul: Simulating 9 years, 44.69 weeks as 123 report steps
╭────────────────┬───────────┬───────────────┬──────────╮
 Iteration type   Avg/step   Avg/ministep     Total 
 123 steps  141 ministeps  (wasted) 
├────────────────┼───────────┼───────────────┼──────────┤
 Newton         │    3.3252 │       2.90071 │  409 (0) │
 Linearization  │   4.47154 │       3.90071 │  550 (0) │
 Linear solver  │    10.252 │       8.94326 │ 1261 (0) │
 Precond apply  │   20.5041 │       17.8865 │ 2522 (0) │
╰────────────────┴───────────┴───────────────┴──────────╯
╭───────────────┬────────┬────────────┬────────╮
 Timing type      Each    Relative   Total 
     ms  Percentage       s 
├───────────────┼────────┼────────────┼────────┤
 Properties    │ 0.2486 │     5.57 % │ 0.1017 │
 Equations     │ 0.2148 │     6.47 % │ 0.1182 │
 Assembly      │ 0.2563 │     7.72 % │ 0.1410 │
 Linear solve  │ 0.3007 │     6.73 % │ 0.1230 │
 Linear setup  │ 1.8231 │    40.82 % │ 0.7457 │
 Precond apply │ 0.1967 │    27.15 % │ 0.4960 │
 Update        │ 0.0773 │     1.73 % │ 0.0316 │
 Convergence   │ 0.0684 │     2.06 % │ 0.0376 │
 Input/Output  │ 0.0299 │     0.23 % │ 0.0042 │
 Other         │ 0.0681 │     1.53 % │ 0.0279 │
├───────────────┼────────┼────────────┼────────┤
 Total         │ 4.4666 │   100.00 % │ 1.8268 │
╰───────────────┴────────┴────────────┴────────╯

Plot the solution of the base case

We observe a radial flow pattern initially, before coning occurs near the producer well once the fluid has reached the opposite corner. The uniform permeability and porosity gives axial symmetry at x=y.

julia
using GLMakie
to_2d(x) = reshape(vec(x), nx, nx)
get_sat(state) = to_2d(state[:Saturations][2, :])
nt = length(report_time)
fig = Figure()
h = nothing
ax = Axis(fig[1, 1])
h = contourf!(ax, get_sat(states[nt÷3]))
ax = Axis(fig[1, 2])
h = contourf!(ax, get_sat(states[nt]))
Colorbar(fig[1, end+1], h)
fig

Create 10 realizations

We create a small set of realizations of the same model, with porosity that is uniformly varying between 0.05 and 0.3. This is not especially sophisticated geostatistics - for a more realistic approach, take a look at GeoStats.jl. The main idea is to get significantly different flow patterns as the porosity and permeability changes.

julia
N = 10
saturations = []
wells = []
report_step = nt
for i = 1:N
    poro = 0.05 .+ 0.25*rand(Float64, (nx*nx))
    ws_i, states_i, rt = simulate_qfs(poro)
    push!(wells, ws_i)
    push!(saturations, get_sat(states_i[report_step]))
end
Jutul: Simulating 9 years, 44.69 weeks as 123 report steps
╭────────────────┬───────────┬───────────────┬────────────╮
 Iteration type   Avg/step   Avg/ministep       Total 
 123 steps  195 ministeps    (wasted) 
├────────────────┼───────────┼───────────────┼────────────┤
 Newton         │   7.08943 │       4.47179 │   872 (75) │
 Linearization  │    8.6748 │       5.47179 │  1067 (80) │
 Linear solver  │   17.0976 │       10.7846 │ 2103 (114) │
 Precond apply  │   34.1951 │       21.5692 │ 4206 (228) │
╰────────────────┴───────────┴───────────────┴────────────╯
╭───────────────┬────────┬────────────┬────────╮
 Timing type      Each    Relative   Total 
     ms  Percentage       s 
├───────────────┼────────┼────────────┼────────┤
 Properties    │ 0.2536 │     5.10 % │ 0.2212 │
 Equations     │ 0.2242 │     5.51 % │ 0.2392 │
 Assembly      │ 0.2643 │     6.50 % │ 0.2820 │
 Linear solve  │ 0.2708 │     5.44 % │ 0.2362 │
 Linear setup  │ 2.1008 │    42.21 % │ 1.8319 │
 Precond apply │ 0.1935 │    18.75 % │ 0.8140 │
 Update        │ 0.0833 │     1.67 % │ 0.0726 │
 Convergence   │ 0.5006 │    12.31 % │ 0.5341 │
 Input/Output  │ 0.0287 │     0.13 % │ 0.0056 │
 Other         │ 0.1185 │     2.38 % │ 0.1034 │
├───────────────┼────────┼────────────┼────────┤
 Total         │ 4.9772 │   100.00 % │ 4.3401 │
╰───────────────┴────────┴────────────┴────────╯
Jutul: Simulating 9 years, 44.69 weeks as 123 report steps
╭────────────────┬───────────┬───────────────┬──────────╮
 Iteration type   Avg/step   Avg/ministep     Total 
 123 steps  145 ministeps  (wasted) 
├────────────────┼───────────┼───────────────┼──────────┤
 Newton         │   3.47967 │       2.95172 │  428 (0) │
 Linearization  │   4.65854 │       3.95172 │  573 (0) │
 Linear solver  │   14.8943 │       12.6345 │ 1832 (0) │
 Precond apply  │   29.7886 │        25.269 │ 3664 (0) │
╰────────────────┴───────────┴───────────────┴──────────╯
╭───────────────┬────────┬────────────┬────────╮
 Timing type      Each    Relative   Total 
     ms  Percentage       s 
├───────────────┼────────┼────────────┼────────┤
 Properties    │ 0.2491 │     4.41 % │ 0.1066 │
 Equations     │ 0.2186 │     5.18 % │ 0.1253 │
 Assembly      │ 0.4702 │    11.13 % │ 0.2694 │
 Linear solve  │ 0.3827 │     6.77 % │ 0.1638 │
 Linear setup  │ 2.2334 │    39.50 % │ 0.9559 │
 Precond apply │ 0.1890 │    28.61 % │ 0.6924 │
 Update        │ 0.0790 │     1.40 % │ 0.0338 │
 Convergence   │ 0.0685 │     1.62 % │ 0.0393 │
 Input/Output  │ 0.0308 │     0.18 % │ 0.0045 │
 Other         │ 0.0685 │     1.21 % │ 0.0293 │
├───────────────┼────────┼────────────┼────────┤
 Total         │ 5.6549 │   100.00 % │ 2.4203 │
╰───────────────┴────────┴────────────┴────────╯
Jutul: Simulating 9 years, 44.69 weeks as 123 report steps
╭────────────────┬───────────┬───────────────┬──────────╮
 Iteration type   Avg/step   Avg/ministep     Total 
 123 steps  156 ministeps  (wasted) 
├────────────────┼───────────┼───────────────┼──────────┤
 Newton         │   4.21138 │       3.32051 │  518 (0) │
 Linearization  │   5.47967 │       4.32051 │  674 (0) │
 Linear solver  │   20.3496 │       16.0449 │ 2503 (0) │
 Precond apply  │   40.6992 │       32.0897 │ 5006 (0) │
╰────────────────┴───────────┴───────────────┴──────────╯
╭───────────────┬────────┬────────────┬────────╮
 Timing type      Each    Relative   Total 
     ms  Percentage       s 
├───────────────┼────────┼────────────┼────────┤
 Properties    │ 0.2481 │     4.48 % │ 0.1285 │
 Equations     │ 0.2151 │     5.05 % │ 0.1449 │
 Assembly      │ 0.2570 │     6.04 % │ 0.1732 │
 Linear solve  │ 0.4156 │     7.50 % │ 0.2153 │
 Linear setup  │ 2.2172 │    40.03 % │ 1.1485 │
 Precond apply │ 0.1869 │    32.62 % │ 0.9357 │
 Update        │ 0.0771 │     1.39 % │ 0.0399 │
 Convergence   │ 0.0681 │     1.60 % │ 0.0459 │
 Input/Output  │ 0.0304 │     0.17 % │ 0.0048 │
 Other         │ 0.0620 │     1.12 % │ 0.0321 │
├───────────────┼────────┼────────────┼────────┤
 Total         │ 5.5383 │   100.00 % │ 2.8688 │
╰───────────────┴────────┴────────────┴────────╯
Jutul: Simulating 9 years, 44.69 weeks as 123 report steps
╭────────────────┬───────────┬───────────────┬────────────╮
 Iteration type   Avg/step   Avg/ministep       Total 
 123 steps  205 ministeps    (wasted) 
├────────────────┼───────────┼───────────────┼────────────┤
 Newton         │   7.88618 │       4.73171 │   970 (30) │
 Linearization  │   9.55285 │       5.73171 │  1175 (32) │
 Linear solver  │   29.7561 │       17.8537 │ 3660 (120) │
 Precond apply  │   59.5122 │       35.7073 │ 7320 (240) │
╰────────────────┴───────────┴───────────────┴────────────╯
╭───────────────┬────────┬────────────┬────────╮
 Timing type      Each    Relative   Total 
     ms  Percentage       s 
├───────────────┼────────┼────────────┼────────┤
 Properties    │ 0.2505 │     4.89 % │ 0.2430 │
 Equations     │ 0.2841 │     6.72 % │ 0.3339 │
 Assembly      │ 0.2613 │     6.18 % │ 0.3070 │
 Linear solve  │ 0.3480 │     6.80 % │ 0.3376 │
 Linear setup  │ 2.2264 │    43.48 % │ 2.1596 │
 Precond apply │ 0.1876 │    27.65 % │ 1.3733 │
 Update        │ 0.0793 │     1.55 % │ 0.0769 │
 Convergence   │ 0.0706 │     1.67 % │ 0.0830 │
 Input/Output  │ 0.0270 │     0.11 % │ 0.0055 │
 Other         │ 0.0480 │     0.94 % │ 0.0466 │
├───────────────┼────────┼────────────┼────────┤
 Total         │ 5.1199 │   100.00 % │ 4.9663 │
╰───────────────┴────────┴────────────┴────────╯
Jutul: Simulating 9 years, 44.69 weeks as 123 report steps
╭────────────────┬───────────┬───────────────┬────────────╮
 Iteration type   Avg/step   Avg/ministep       Total 
 123 steps  180 ministeps    (wasted) 
├────────────────┼───────────┼───────────────┼────────────┤
 Newton         │   5.86992 │       4.01111 │   722 (15) │
 Linearization  │   7.33333 │       5.01111 │   902 (16) │
 Linear solver  │   29.8293 │       20.3833 │ 3669 (106) │
 Precond apply  │   59.6585 │       40.7667 │ 7338 (212) │
╰────────────────┴───────────┴───────────────┴────────────╯
╭───────────────┬────────┬────────────┬────────╮
 Timing type      Each    Relative   Total 
     ms  Percentage       s 
├───────────────┼────────┼────────────┼────────┤
 Properties    │ 0.2494 │     4.38 % │ 0.1801 │
 Equations     │ 0.2159 │     4.74 % │ 0.1947 │
 Assembly      │ 0.2589 │     5.68 % │ 0.2336 │
 Linear solve  │ 0.4318 │     7.58 % │ 0.3118 │
 Linear setup  │ 2.2839 │    40.11 % │ 1.6490 │
 Precond apply │ 0.1878 │    33.51 % │ 1.3780 │
 Update        │ 0.0788 │     1.38 % │ 0.0569 │
 Convergence   │ 0.0697 │     1.53 % │ 0.0629 │
 Input/Output  │ 0.0290 │     0.13 % │ 0.0052 │
 Other         │ 0.0548 │     0.96 % │ 0.0396 │
├───────────────┼────────┼────────────┼────────┤
 Total         │ 5.6949 │   100.00 % │ 4.1117 │
╰───────────────┴────────┴────────────┴────────╯
Jutul: Simulating 9 years, 44.69 weeks as 123 report steps
╭────────────────┬───────────┬───────────────┬───────────╮
 Iteration type   Avg/step   Avg/ministep      Total 
 123 steps  170 ministeps   (wasted) 
├────────────────┼───────────┼───────────────┼───────────┤
 Newton         │   5.15447 │       3.72941 │  634 (30) │
 Linearization  │   6.53659 │       4.72941 │  804 (32) │
 Linear solver  │   12.6585 │       9.15882 │ 1557 (30) │
 Precond apply  │   25.3171 │       18.3176 │ 3114 (60) │
╰────────────────┴───────────┴───────────────┴───────────╯
╭───────────────┬────────┬────────────┬────────╮
 Timing type      Each    Relative   Total 
     ms  Percentage       s 
├───────────────┼────────┼────────────┼────────┤
 Properties    │ 0.2470 │     5.56 % │ 0.1566 │
 Equations     │ 0.2136 │     6.10 % │ 0.1717 │
 Assembly      │ 0.2570 │     7.34 % │ 0.2066 │
 Linear solve  │ 0.2702 │     6.08 % │ 0.1713 │
 Linear setup  │ 2.0960 │    47.17 % │ 1.3289 │
 Precond apply │ 0.1915 │    21.17 % │ 0.5962 │
 Update        │ 0.0769 │     1.73 % │ 0.0488 │
 Convergence   │ 0.0682 │     1.95 % │ 0.0548 │
 Input/Output  │ 0.2718 │     1.64 % │ 0.0462 │
 Other         │ 0.0566 │     1.27 % │ 0.0359 │
├───────────────┼────────┼────────────┼────────┤
 Total         │ 4.4433 │   100.00 % │ 2.8171 │
╰───────────────┴────────┴────────────┴────────╯
Jutul: Simulating 9 years, 44.69 weeks as 123 report steps
╭────────────────┬───────────┬───────────────┬────────────╮
 Iteration type   Avg/step   Avg/ministep       Total 
 123 steps  174 ministeps    (wasted) 
├────────────────┼───────────┼───────────────┼────────────┤
 Newton         │   5.30081 │       3.74713 │   652 (15) │
 Linearization  │   6.71545 │       4.74713 │   826 (16) │
 Linear solver  │   25.9024 │       18.3103 │  3186 (87) │
 Precond apply  │   51.8049 │       36.6207 │ 6372 (174) │
╰────────────────┴───────────┴───────────────┴────────────╯
╭───────────────┬────────┬────────────┬────────╮
 Timing type      Each    Relative   Total 
     ms  Percentage       s 
├───────────────┼────────┼────────────┼────────┤
 Properties    │ 0.2487 │     4.48 % │ 0.1621 │
 Equations     │ 0.2160 │     4.93 % │ 0.1784 │
 Assembly      │ 0.2575 │     5.88 % │ 0.2127 │
 Linear solve  │ 0.4186 │     7.54 % │ 0.2730 │
 Linear setup  │ 2.2315 │    40.21 % │ 1.4550 │
 Precond apply │ 0.1860 │    32.76 % │ 1.1852 │
 Update        │ 0.0782 │     1.41 % │ 0.0510 │
 Convergence   │ 0.0702 │     1.60 % │ 0.0580 │
 Input/Output  │ 0.0301 │     0.14 % │ 0.0052 │
 Other         │ 0.0577 │     1.04 % │ 0.0376 │
├───────────────┼────────┼────────────┼────────┤
 Total         │ 5.5494 │   100.00 % │ 3.6182 │
╰───────────────┴────────┴────────────┴────────╯
Jutul: Simulating 9 years, 44.69 weeks as 123 report steps
╭────────────────┬───────────┬───────────────┬─────────────╮
 Iteration type   Avg/step   Avg/ministep        Total 
 123 steps  372 ministeps     (wasted) 
├────────────────┼───────────┼───────────────┼─────────────┤
 Newton         │      21.0 │       6.94355 │ 2583 (1335) │
 Linearization  │   24.0244 │       7.94355 │ 2955 (1424) │
 Linear solver  │    31.252 │       10.3333 │ 3844 (1391) │
 Precond apply  │   62.5041 │       20.6667 │ 7688 (2782) │
╰────────────────┴───────────┴───────────────┴─────────────╯
╭───────────────┬────────┬────────────┬────────╮
 Timing type      Each    Relative   Total 
     ms  Percentage       s 
├───────────────┼────────┼────────────┼────────┤
 Properties    │ 0.2443 │     6.81 % │ 0.6309 │
 Equations     │ 0.2308 │     7.36 % │ 0.6820 │
 Assembly      │ 0.2581 │     8.23 % │ 0.7626 │
 Linear solve  │ 0.2344 │     6.53 % │ 0.6054 │
 Linear setup  │ 1.7789 │    49.58 % │ 4.5948 │
 Precond apply │ 0.1909 │    15.84 % │ 1.4678 │
 Update        │ 0.0762 │     2.12 % │ 0.1969 │
 Convergence   │ 0.0709 │     2.26 % │ 0.2095 │
 Input/Output  │ 0.0239 │     0.10 % │ 0.0089 │
 Other         │ 0.0421 │     1.17 % │ 0.1087 │
├───────────────┼────────┼────────────┼────────┤
 Total         │ 3.5878 │   100.00 % │ 9.2673 │
╰───────────────┴────────┴────────────┴────────╯
Jutul: Simulating 9 years, 44.69 weeks as 123 report steps
╭────────────────┬───────────┬───────────────┬──────────╮
 Iteration type   Avg/step   Avg/ministep     Total 
 123 steps  159 ministeps  (wasted) 
├────────────────┼───────────┼───────────────┼──────────┤
 Newton         │   4.42276 │       3.42138 │  544 (0) │
 Linearization  │   5.71545 │       4.42138 │  703 (0) │
 Linear solver  │    17.065 │       13.2013 │ 2099 (0) │
 Precond apply  │   34.1301 │       26.4025 │ 4198 (0) │
╰────────────────┴───────────┴───────────────┴──────────╯
╭───────────────┬────────┬────────────┬────────╮
 Timing type      Each    Relative   Total 
     ms  Percentage       s 
├───────────────┼────────┼────────────┼────────┤
 Properties    │ 0.2512 │     4.71 % │ 0.1367 │
 Equations     │ 0.2221 │     5.38 % │ 0.1562 │
 Assembly      │ 0.2611 │     6.32 % │ 0.1836 │
 Linear solve  │ 0.3637 │     6.81 % │ 0.1979 │
 Linear setup  │ 2.3583 │    44.18 % │ 1.2829 │
 Precond apply │ 0.1917 │    27.71 % │ 0.8048 │
 Update        │ 0.0824 │     1.54 % │ 0.0449 │
 Convergence   │ 0.0744 │     1.80 % │ 0.0523 │
 Input/Output  │ 0.0375 │     0.21 % │ 0.0060 │
 Other         │ 0.0713 │     1.34 % │ 0.0388 │
├───────────────┼────────┼────────────┼────────┤
 Total         │ 5.3379 │   100.00 % │ 2.9038 │
╰───────────────┴────────┴────────────┴────────╯
Jutul: Simulating 9 years, 44.69 weeks as 123 report steps
╭────────────────┬───────────┬───────────────┬──────────╮
 Iteration type   Avg/step   Avg/ministep     Total 
 123 steps  153 ministeps  (wasted) 
├────────────────┼───────────┼───────────────┼──────────┤
 Newton         │   4.03252 │       3.24183 │  496 (0) │
 Linearization  │   5.27642 │       4.24183 │  649 (0) │
 Linear solver  │   17.1789 │       13.8105 │ 2113 (0) │
 Precond apply  │   34.3577 │       27.6209 │ 4226 (0) │
╰────────────────┴───────────┴───────────────┴──────────╯
╭───────────────┬────────┬────────────┬────────╮
 Timing type      Each    Relative   Total 
     ms  Percentage       s 
├───────────────┼────────┼────────────┼────────┤
 Properties    │ 0.2476 │     4.61 % │ 0.1228 │
 Equations     │ 0.2148 │     5.23 % │ 0.1394 │
 Assembly      │ 0.2573 │     6.27 % │ 0.1670 │
 Linear solve  │ 0.4623 │     8.61 % │ 0.2293 │
 Linear setup  │ 2.2151 │    41.25 % │ 1.0987 │
 Precond apply │ 0.1860 │    29.51 % │ 0.7860 │
 Update        │ 0.0781 │     1.45 % │ 0.0387 │
 Convergence   │ 0.0698 │     1.70 % │ 0.0453 │
 Input/Output  │ 0.0302 │     0.17 % │ 0.0046 │
 Other         │ 0.0631 │     1.18 % │ 0.0313 │
├───────────────┼────────┼────────────┼────────┤
 Total         │ 5.3694 │   100.00 % │ 2.6632 │
╰───────────────┴────────┴────────────┴────────╯

Plot the oil rate at the producer over the ensemble

julia
using Statistics
fig = Figure()
ax = Axis(fig[1, 1])
for i = 1:N
    ws = wells[i]
    q = -ws[:Producer][:orat]
    lines!(ax, report_time, q)
end
xlims!(ax, [mean(report_time), report_time[end]])
ylims!(ax, 0, 0.0075)
fig

Plot the average saturation over the ensemble

julia
avg = mean(saturations)
fig = Figure()
h = nothing
ax = Axis(fig[1, 1])
h = contourf!(ax, avg)
fig

Plot the isocontour lines over the ensemble

julia
fig = Figure()
h = nothing
ax = Axis(fig[1, 1])
for s in saturations
    contour!(ax, s, levels = 0:0.1:1)
end
fig

Example on GitHub

If you would like to run this example yourself, it can be downloaded from the JutulDarcy.jl GitHub repository as a script, or as a Jupyter Notebook

This example took 50.622666306 seconds to complete.

This page was generated using Literate.jl.