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.2473 │     5.59 % │ 0.1012 │
 Equations     │ 0.1981 │     6.02 % │ 0.1089 │
 Assembly      │ 0.2566 │     7.81 % │ 0.1411 │
 Linear solve  │ 0.2922 │     6.61 % │ 0.1195 │
 Linear setup  │ 1.8158 │    41.07 % │ 0.7427 │
 Precond apply │ 0.1958 │    27.30 % │ 0.4937 │
 Update        │ 0.0768 │     1.74 % │ 0.0314 │
 Convergence   │ 0.0687 │     2.09 % │ 0.0378 │
 Input/Output  │ 0.0307 │     0.24 % │ 0.0043 │
 Other         │ 0.0676 │     1.53 % │ 0.0277 │
├───────────────┼────────┼────────────┼────────┤
 Total         │ 4.4213 │   100.00 % │ 1.8083 │
╰───────────────┴────────┴────────────┴────────╯

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  149 ministeps  (wasted) 
├────────────────┼───────────┼───────────────┼──────────┤
 Newton         │   3.86179 │       3.18792 │  475 (0) │
 Linearization  │   5.07317 │       4.18792 │  624 (0) │
 Linear solver  │   14.5285 │       11.9933 │ 1787 (0) │
 Precond apply  │   29.0569 │       23.9866 │ 3574 (0) │
╰────────────────┴───────────┴───────────────┴──────────╯
╭───────────────┬────────┬────────────┬────────╮
 Timing type      Each    Relative   Total 
     ms  Percentage       s 
├───────────────┼────────┼────────────┼────────┤
 Properties    │ 0.2476 │     4.90 % │ 0.1176 │
 Equations     │ 0.2145 │     5.57 % │ 0.1339 │
 Assembly      │ 0.2595 │     6.74 % │ 0.1620 │
 Linear solve  │ 0.3357 │     6.64 % │ 0.1594 │
 Linear setup  │ 2.1915 │    43.35 % │ 1.0410 │
 Precond apply │ 0.1878 │    27.95 % │ 0.6712 │
 Update        │ 0.0800 │     1.58 % │ 0.0380 │
 Convergence   │ 0.0713 │     1.85 % │ 0.0445 │
 Input/Output  │ 0.0287 │     0.18 % │ 0.0043 │
 Other         │ 0.0625 │     1.24 % │ 0.0297 │
├───────────────┼────────┼────────────┼────────┤
 Total         │ 5.0558 │   100.00 % │ 2.4015 │
╰───────────────┴────────┴────────────┴────────╯
Jutul: Simulating 9 years, 44.69 weeks as 123 report steps
╭────────────────┬───────────┬───────────────┬───────────╮
 Iteration type   Avg/step   Avg/ministep      Total 
 123 steps  164 ministeps   (wasted) 
├────────────────┼───────────┼───────────────┼───────────┤
 Newton         │   4.98374 │        3.7378 │  613 (15) │
 Linearization  │   6.31707 │        4.7378 │  777 (16) │
 Linear solver  │   12.8049 │       9.60366 │ 1575 (15) │
 Precond apply  │   25.6098 │       19.2073 │ 3150 (30) │
╰────────────────┴───────────┴───────────────┴───────────╯
╭───────────────┬────────┬────────────┬────────╮
 Timing type      Each    Relative   Total 
     ms  Percentage       s 
├───────────────┼────────┼────────────┼────────┤
 Properties    │ 0.2536 │     4.78 % │ 0.1554 │
 Equations     │ 0.2248 │     5.37 % │ 0.1747 │
 Assembly      │ 0.2655 │     6.34 % │ 0.2063 │
 Linear solve  │ 0.2691 │     5.07 % │ 0.1650 │
 Linear setup  │ 2.1208 │    39.94 % │ 1.3000 │
 Precond apply │ 0.1909 │    18.48 % │ 0.6014 │
 Update        │ 0.0836 │     1.58 % │ 0.0513 │
 Convergence   │ 0.6507 │    15.53 % │ 0.5056 │
 Input/Output  │ 0.0295 │     0.15 % │ 0.0048 │
 Other         │ 0.1475 │     2.78 % │ 0.0904 │
├───────────────┼────────┼────────────┼────────┤
 Total         │ 5.3098 │   100.00 % │ 3.2549 │
╰───────────────┴────────┴────────────┴────────╯
Jutul: Simulating 9 years, 44.69 weeks as 123 report steps
╭────────────────┬───────────┬───────────────┬──────────╮
 Iteration type   Avg/step   Avg/ministep     Total 
 123 steps  146 ministeps  (wasted) 
├────────────────┼───────────┼───────────────┼──────────┤
 Newton         │   3.60976 │        3.0411 │  444 (0) │
 Linearization  │   4.79675 │        4.0411 │  590 (0) │
 Linear solver  │   16.6423 │       14.0205 │ 2047 (0) │
 Precond apply  │   33.2846 │       28.0411 │ 4094 (0) │
╰────────────────┴───────────┴───────────────┴──────────╯
╭───────────────┬────────┬────────────┬────────╮
 Timing type      Each    Relative   Total 
     ms  Percentage       s 
├───────────────┼────────┼────────────┼────────┤
 Properties    │ 0.2491 │     4.36 % │ 0.1106 │
 Equations     │ 0.2156 │     5.02 % │ 0.1272 │
 Assembly      │ 0.2588 │     6.02 % │ 0.1527 │
 Linear solve  │ 0.3920 │     6.86 % │ 0.1740 │
 Linear setup  │ 2.2129 │    38.74 % │ 0.9825 │
 Precond apply │ 0.1859 │    30.00 % │ 0.7609 │
 Update        │ 0.3349 │     5.86 % │ 0.1487 │
 Convergence   │ 0.0724 │     1.69 % │ 0.0427 │
 Input/Output  │ 0.0362 │     0.21 % │ 0.0053 │
 Other         │ 0.0706 │     1.24 % │ 0.0314 │
├───────────────┼────────┼────────────┼────────┤
 Total         │ 5.7118 │   100.00 % │ 2.5360 │
╰───────────────┴────────┴────────────┴────────╯
Jutul: Simulating 9 years, 44.69 weeks as 123 report steps
╭────────────────┬───────────┬───────────────┬──────────╮
 Iteration type   Avg/step   Avg/ministep     Total 
 123 steps  151 ministeps  (wasted) 
├────────────────┼───────────┼───────────────┼──────────┤
 Newton         │   3.97561 │       3.23841 │  489 (0) │
 Linearization  │   5.20325 │       4.23841 │  640 (0) │
 Linear solver  │   15.6504 │       12.7483 │ 1925 (0) │
 Precond apply  │   31.3008 │       25.4967 │ 3850 (0) │
╰────────────────┴───────────┴───────────────┴──────────╯
╭───────────────┬────────┬────────────┬────────╮
 Timing type      Each    Relative   Total 
     ms  Percentage       s 
├───────────────┼────────┼────────────┼────────┤
 Properties    │ 0.2529 │     4.80 % │ 0.1237 │
 Equations     │ 0.2218 │     5.51 % │ 0.1419 │
 Assembly      │ 0.2592 │     6.44 % │ 0.1659 │
 Linear solve  │ 0.3604 │     6.84 % │ 0.1763 │
 Linear setup  │ 2.2673 │    43.02 % │ 1.1087 │
 Precond apply │ 0.1899 │    28.36 % │ 0.7309 │
 Update        │ 0.0813 │     1.54 % │ 0.0398 │
 Convergence   │ 0.0747 │     1.86 % │ 0.0478 │
 Input/Output  │ 0.0383 │     0.22 % │ 0.0058 │
 Other         │ 0.0749 │     1.42 % │ 0.0366 │
├───────────────┼────────┼────────────┼────────┤
 Total         │ 5.2707 │   100.00 % │ 2.5774 │
╰───────────────┴────────┴────────────┴────────╯
Jutul: Simulating 9 years, 44.69 weeks as 123 report steps
╭────────────────┬───────────┬───────────────┬──────────╮
 Iteration type   Avg/step   Avg/ministep     Total 
 123 steps  149 ministeps  (wasted) 
├────────────────┼───────────┼───────────────┼──────────┤
 Newton         │   3.78862 │       3.12752 │  466 (0) │
 Linearization  │       5.0 │       4.12752 │  615 (0) │
 Linear solver  │   18.7317 │       15.4631 │ 2304 (0) │
 Precond apply  │   37.4634 │       30.9262 │ 4608 (0) │
╰────────────────┴───────────┴───────────────┴──────────╯
╭───────────────┬────────┬────────────┬────────╮
 Timing type      Each    Relative   Total 
     ms  Percentage       s 
├───────────────┼────────┼────────────┼────────┤
 Properties    │ 0.2578 │     4.49 % │ 0.1201 │
 Equations     │ 0.2343 │     5.38 % │ 0.1441 │
 Assembly      │ 0.2682 │     6.16 % │ 0.1650 │
 Linear solve  │ 0.4305 │     7.50 % │ 0.2006 │
 Linear setup  │ 2.2661 │    39.46 % │ 1.0560 │
 Precond apply │ 0.1858 │    31.98 % │ 0.8560 │
 Update        │ 0.0890 │     1.55 % │ 0.0415 │
 Convergence   │ 0.0786 │     1.81 % │ 0.0484 │
 Input/Output  │ 0.0427 │     0.24 % │ 0.0064 │
 Other         │ 0.0821 │     1.43 % │ 0.0382 │
├───────────────┼────────┼────────────┼────────┤
 Total         │ 5.7428 │   100.00 % │ 2.6762 │
╰───────────────┴────────┴────────────┴────────╯
Jutul: Simulating 9 years, 44.69 weeks as 123 report steps
╭────────────────┬───────────┬───────────────┬──────────╮
 Iteration type   Avg/step   Avg/ministep     Total 
 123 steps  147 ministeps  (wasted) 
├────────────────┼───────────┼───────────────┼──────────┤
 Newton         │   3.65041 │       3.05442 │  449 (0) │
 Linearization  │   4.84553 │       4.05442 │  596 (0) │
 Linear solver  │   10.7398 │       8.98639 │ 1321 (0) │
 Precond apply  │   21.4797 │       17.9728 │ 2642 (0) │
╰────────────────┴───────────┴───────────────┴──────────╯
╭───────────────┬────────┬────────────┬────────╮
 Timing type      Each    Relative   Total 
     ms  Percentage       s 
├───────────────┼────────┼────────────┼────────┤
 Properties    │ 0.4385 │     8.85 % │ 0.1969 │
 Equations     │ 0.2190 │     5.87 % │ 0.1305 │
 Assembly      │ 0.2600 │     6.96 % │ 0.1550 │
 Linear solve  │ 0.2928 │     5.91 % │ 0.1315 │
 Linear setup  │ 2.2066 │    44.53 % │ 0.9908 │
 Precond apply │ 0.1896 │    22.51 % │ 0.5008 │
 Update        │ 0.0811 │     1.64 % │ 0.0364 │
 Convergence   │ 0.0725 │     1.94 % │ 0.0432 │
 Input/Output  │ 0.0453 │     0.30 % │ 0.0067 │
 Other         │ 0.0744 │     1.50 % │ 0.0334 │
├───────────────┼────────┼────────────┼────────┤
 Total         │ 4.9557 │   100.00 % │ 2.2251 │
╰───────────────┴────────┴────────────┴────────╯
Jutul: Simulating 9 years, 44.69 weeks as 123 report steps
╭────────────────┬───────────┬───────────────┬────────────╮
 Iteration type   Avg/step   Avg/ministep       Total 
 123 steps  186 ministeps    (wasted) 
├────────────────┼───────────┼───────────────┼────────────┤
 Newton         │   6.13821 │       4.05914 │   755 (30) │
 Linearization  │   7.65041 │       5.05914 │   941 (32) │
 Linear solver  │   28.1626 │       18.6237 │ 3464 (137) │
 Precond apply  │   56.3252 │       37.2473 │ 6928 (274) │
╰────────────────┴───────────┴───────────────┴────────────╯
╭───────────────┬────────┬────────────┬────────╮
 Timing type      Each    Relative   Total 
     ms  Percentage       s 
├───────────────┼────────┼────────────┼────────┤
 Properties    │ 0.2546 │     4.66 % │ 0.1922 │
 Equations     │ 0.2214 │     5.05 % │ 0.2083 │
 Assembly      │ 0.2602 │     5.94 % │ 0.2448 │
 Linear solve  │ 0.3985 │     7.30 % │ 0.3009 │
 Linear setup  │ 2.2502 │    41.20 % │ 1.6989 │
 Precond apply │ 0.1862 │    31.29 % │ 1.2903 │
 Update        │ 0.0823 │     1.51 % │ 0.0621 │
 Convergence   │ 0.0753 │     1.72 % │ 0.0708 │
 Input/Output  │ 0.0355 │     0.16 % │ 0.0066 │
 Other         │ 0.0646 │     1.18 % │ 0.0488 │
├───────────────┼────────┼────────────┼────────┤
 Total         │ 5.4621 │   100.00 % │ 4.1239 │
╰───────────────┴────────┴────────────┴────────╯
Jutul: Simulating 9 years, 44.69 weeks as 123 report steps
╭────────────────┬───────────┬───────────────┬──────────╮
 Iteration type   Avg/step   Avg/ministep     Total 
 123 steps  152 ministeps  (wasted) 
├────────────────┼───────────┼───────────────┼──────────┤
 Newton         │    3.7561 │       3.03947 │  462 (0) │
 Linearization  │   4.99187 │       4.03947 │  614 (0) │
 Linear solver  │   16.7073 │       13.5197 │ 2055 (0) │
 Precond apply  │   33.4146 │       27.0395 │ 4110 (0) │
╰────────────────┴───────────┴───────────────┴──────────╯
╭───────────────┬────────┬────────────┬────────╮
 Timing type      Each    Relative   Total 
     ms  Percentage       s 
├───────────────┼────────┼────────────┼────────┤
 Properties    │ 0.2510 │     4.54 % │ 0.1160 │
 Equations     │ 0.2199 │     5.29 % │ 0.1350 │
 Assembly      │ 0.2608 │     6.27 % │ 0.1602 │
 Linear solve  │ 0.3793 │     6.86 % │ 0.1752 │
 Linear setup  │ 2.2352 │    40.44 % │ 1.0327 │
 Precond apply │ 0.1873 │    30.14 % │ 0.7697 │
 Update        │ 0.0818 │     1.48 % │ 0.0378 │
 Convergence   │ 0.0723 │     1.74 % │ 0.0444 │
 Input/Output  │ 0.0340 │     0.20 % │ 0.0052 │
 Other         │ 0.1675 │     3.03 % │ 0.0774 │
├───────────────┼────────┼────────────┼────────┤
 Total         │ 5.5269 │   100.00 % │ 2.5534 │
╰───────────────┴────────┴────────────┴────────╯
Jutul: Simulating 9 years, 44.69 weeks as 123 report steps
╭────────────────┬───────────┬───────────────┬─────────────╮
 Iteration type   Avg/step   Avg/ministep        Total 
 123 steps  258 ministeps     (wasted) 
├────────────────┼───────────┼───────────────┼─────────────┤
 Newton         │   11.4959 │       5.48062 │   1414 (30) │
 Linearization  │   13.5935 │       6.48062 │   1672 (32) │
 Linear solver  │   52.0976 │       24.8372 │  6408 (152) │
 Precond apply  │   104.195 │       49.6744 │ 12816 (304) │
╰────────────────┴───────────┴───────────────┴─────────────╯
╭───────────────┬────────┬────────────┬────────╮
 Timing type      Each    Relative   Total 
     ms  Percentage       s 
├───────────────┼────────┼────────────┼────────┤
 Properties    │ 0.2504 │     4.63 % │ 0.3541 │
 Equations     │ 0.2192 │     4.79 % │ 0.3665 │
 Assembly      │ 0.2600 │     5.69 % │ 0.4347 │
 Linear solve  │ 0.3893 │     7.20 % │ 0.5504 │
 Linear setup  │ 2.2900 │    42.35 % │ 3.2380 │
 Precond apply │ 0.1864 │    31.23 % │ 2.3883 │
 Update        │ 0.0812 │     1.50 % │ 0.1148 │
 Convergence   │ 0.0748 │     1.64 % │ 0.1251 │
 Input/Output  │ 0.0275 │     0.09 % │ 0.0071 │
 Other         │ 0.0479 │     0.89 % │ 0.0677 │
├───────────────┼────────┼────────────┼────────┤
 Total         │ 5.4079 │   100.00 % │ 7.6467 │
╰───────────────┴────────┴────────────┴────────╯
Jutul: Simulating 9 years, 44.69 weeks as 123 report steps
╭────────────────┬───────────┬───────────────┬────────────╮
 Iteration type   Avg/step   Avg/ministep       Total 
 123 steps  192 ministeps    (wasted) 
├────────────────┼───────────┼───────────────┼────────────┤
 Newton         │   6.72358 │       4.30729 │   827 (90) │
 Linearization  │   8.28455 │       5.30729 │  1019 (96) │
 Linear solver  │   26.0081 │       16.6615 │ 3199 (248) │
 Precond apply  │   52.0163 │       33.3229 │ 6398 (496) │
╰────────────────┴───────────┴───────────────┴────────────╯
╭───────────────┬────────┬────────────┬────────╮
 Timing type      Each    Relative   Total 
     ms  Percentage       s 
├───────────────┼────────┼────────────┼────────┤
 Properties    │ 0.2469 │     4.86 % │ 0.2041 │
 Equations     │ 0.2147 │     5.21 % │ 0.2188 │
 Assembly      │ 0.2581 │     6.26 % │ 0.2630 │
 Linear solve  │ 0.3428 │     6.74 % │ 0.2835 │
 Linear setup  │ 2.1874 │    43.04 % │ 1.8090 │
 Precond apply │ 0.1862 │    28.34 % │ 1.1914 │
 Update        │ 0.0784 │     1.54 % │ 0.0648 │
 Convergence   │ 0.1185 │     2.87 % │ 0.1207 │
 Input/Output  │ 0.0287 │     0.13 % │ 0.0055 │
 Other         │ 0.0514 │     1.01 % │ 0.0425 │
├───────────────┼────────┼────────────┼────────┤
 Total         │ 5.0826 │   100.00 % │ 4.2033 │
╰───────────────┴────────┴────────────┴────────╯

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 44.671602026 seconds to complete.

This page was generated using Literate.jl.