Buckley-Leverett two-phase problem
Introduction Validation ImmiscibleThe Buckley-Leverett test problem is a classical reservoir simulation benchmark that demonstrates the nonlinear displacement process of a viscous fluid being displaced by a less viscous fluid, typically taken to be water displacing oil.
Problem definition
This is a simple model without wells, where the flow is driven by a simple source term and a simple constant pressure boundary condition at the outlet. We define a function that sets up a two-phase system, a simple 1D domain and replaces the default relative permeability functions with quadratic functions:
In addition, the phase viscosities are treated as constant parameters of 1 and 5 centipoise for the displacing and resident fluids, respectively.
The function is parametrized on the number of cells and the number of time-steps used to solve the model. This function, since it uses a relatively simple setup without wells, uses the Jutul functions directly.
using JutulDarcy, Jutul
function solve_bl(;nc = 100, time = 1.0, nstep = nc)
T = time
tstep = repeat([T/nstep], nstep)
domain = get_1d_reservoir(nc)
nc = number_of_cells(domain)
timesteps = tstep*3600*24
bar = 1e5
p0 = 100*bar
sys = ImmiscibleSystem((LiquidPhase(), VaporPhase()))
model = SimulationModel(domain, sys)
kr = BrooksCoreyRelativePermeabilities(sys, [2.0, 2.0], [0.2, 0.2])
replace_variables!(model, RelativePermeabilities = kr)
tot_time = sum(timesteps)
pv = pore_volume(domain)
irate = 500*sum(pv)/tot_time
src = SourceTerm(1, irate, fractional_flow = [1.0, 0.0])
bc = FlowBoundaryCondition(nc, p0/2)
forces = setup_forces(model, sources = src, bc = bc)
parameters = setup_parameters(model, PhaseViscosities = [1e-3, 5e-3]) # 1 and 5 cP
state0 = setup_state(model, Pressure = p0, Saturations = [0.0, 1.0])
states, report = simulate(state0, model, timesteps,
forces = forces, parameters = parameters)
return states, model, report
endsolve_bl (generic function with 1 method)Run the base case
We solve a small model with 100 cells and 100 steps to serve as the baseline.
n, n_f = 100, 1000
states, model, report = solve_bl(nc = n)
print_stats(report)Jutul: Simulating 1 day as 100 report steps
╭────────────────┬───────────┬───────────────┬──────────╮
│ Iteration type │ Avg/step │ Avg/ministep │ Total │
│ │ 100 steps │ 100 ministeps │ (wasted) │
├────────────────┼───────────┼───────────────┼──────────┤
│ Newton │ 4.15 │ 4.15 │ 415 (0) │
│ Linearization │ 5.15 │ 5.15 │ 515 (0) │
│ Linear solver │ 4.15 │ 4.15 │ 415 (0) │
│ Precond apply │ 0.0 │ 0.0 │ 0 (0) │
╰────────────────┴───────────┴───────────────┴──────────╯
╭───────────────┬────────┬────────────┬────────╮
│ Timing type │ Each │ Relative │ Total │
│ │ ms │ Percentage │ s │
├───────────────┼────────┼────────────┼────────┤
│ Properties │ 0.0142 │ 0.23 % │ 0.0059 │
│ Equations │ 0.4465 │ 8.92 % │ 0.2299 │
│ Assembly │ 0.2141 │ 4.28 % │ 0.1102 │
│ Linear solve │ 2.2996 │ 37.04 % │ 0.9543 │
│ Linear setup │ 0.0000 │ 0.00 % │ 0.0000 │
│ Precond apply │ 0.0000 │ 0.00 % │ 0.0000 │
│ Update │ 0.0127 │ 0.20 % │ 0.0053 │
│ Convergence │ 0.1723 │ 3.44 % │ 0.0888 │
│ Input/Output │ 0.2987 │ 1.16 % │ 0.0299 │
│ Other │ 2.7765 │ 44.72 % │ 1.1523 │
├───────────────┼────────┼────────────┼────────┤
│ Total │ 6.2086 │ 100.00 % │ 2.5766 │
╰───────────────┴────────┴────────────┴────────╯
╭────────────────┬───────────┬───────────────┬──────────╮
│ Iteration type │ Avg/step │ Avg/ministep │ Total │
│ │ 100 steps │ 100 ministeps │ (wasted) │
├────────────────┼───────────┼───────────────┼──────────┤
│ Newton │ 4.15 │ 4.15 │ 415 (0) │
│ Linearization │ 5.15 │ 5.15 │ 515 (0) │
│ Linear solver │ 4.15 │ 4.15 │ 415 (0) │
│ Precond apply │ 0.0 │ 0.0 │ 0 (0) │
╰────────────────┴───────────┴───────────────┴──────────╯
╭───────────────┬────────┬────────────┬────────╮
│ Timing type │ Each │ Relative │ Total │
│ │ ms │ Percentage │ s │
├───────────────┼────────┼────────────┼────────┤
│ Properties │ 0.0142 │ 0.23 % │ 0.0059 │
│ Equations │ 0.4465 │ 8.92 % │ 0.2299 │
│ Assembly │ 0.2141 │ 4.28 % │ 0.1102 │
│ Linear solve │ 2.2996 │ 37.04 % │ 0.9543 │
│ Linear setup │ 0.0000 │ 0.00 % │ 0.0000 │
│ Precond apply │ 0.0000 │ 0.00 % │ 0.0000 │
│ Update │ 0.0127 │ 0.20 % │ 0.0053 │
│ Convergence │ 0.1723 │ 3.44 % │ 0.0888 │
│ Input/Output │ 0.2987 │ 1.16 % │ 0.0299 │
│ Other │ 2.7765 │ 44.72 % │ 1.1523 │
├───────────────┼────────┼────────────┼────────┤
│ Total │ 6.2086 │ 100.00 % │ 2.5766 │
╰───────────────┴────────┴────────────┴────────╯Run refined version (1000 cells, 1000 steps)
Using a grid with 100 cells will not yield a fully converged solution. We can increase the number of cells at the cost of increasing the runtime a bit. Note that most of the time is spent in the linear solver, which uses a direct sparse LU factorization by default. For larger problems it is recommended to use an iterative solver. The high-level interface used in later examples automatically sets up an iterative solver with the appropriate preconditioner.
states_refined, _, report_refined = solve_bl(nc = n_f);
print_stats(report_refined)Jutul: Simulating 23 hours, 60 minutes as 1000 report steps
╭────────────────┬────────────┬────────────────┬──────────╮
│ Iteration type │ Avg/step │ Avg/ministep │ Total │
│ │ 1000 steps │ 1000 ministeps │ (wasted) │
├────────────────┼────────────┼────────────────┼──────────┤
│ Newton │ 4.12 │ 4.12 │ 4120 (0) │
│ Linearization │ 5.12 │ 5.12 │ 5120 (0) │
│ Linear solver │ 4.12 │ 4.12 │ 4120 (0) │
│ Precond apply │ 0.0 │ 0.0 │ 0 (0) │
╰────────────────┴────────────┴────────────────┴──────────╯
╭───────────────┬────────┬────────────┬─────────╮
│ Timing type │ Each │ Relative │ Total │
│ │ ms │ Percentage │ s │
├───────────────┼────────┼────────────┼─────────┤
│ Properties │ 0.0807 │ 2.70 % │ 0.3324 │
│ Equations │ 0.0809 │ 3.36 % │ 0.4144 │
│ Assembly │ 0.0345 │ 1.43 % │ 0.1767 │
│ Linear solve │ 2.6155 │ 87.43 % │ 10.7758 │
│ Linear setup │ 0.0000 │ 0.00 % │ 0.0000 │
│ Precond apply │ 0.0000 │ 0.00 % │ 0.0000 │
│ Update │ 0.0464 │ 1.55 % │ 0.1913 │
│ Convergence │ 0.0320 │ 1.33 % │ 0.1640 │
│ Input/Output │ 0.0592 │ 0.48 % │ 0.0592 │
│ Other │ 0.0514 │ 1.72 % │ 0.2116 │
├───────────────┼────────┼────────────┼─────────┤
│ Total │ 2.9916 │ 100.00 % │ 12.3253 │
╰───────────────┴────────┴────────────┴─────────╯
╭────────────────┬────────────┬────────────────┬──────────╮
│ Iteration type │ Avg/step │ Avg/ministep │ Total │
│ │ 1000 steps │ 1000 ministeps │ (wasted) │
├────────────────┼────────────┼────────────────┼──────────┤
│ Newton │ 4.12 │ 4.12 │ 4120 (0) │
│ Linearization │ 5.12 │ 5.12 │ 5120 (0) │
│ Linear solver │ 4.12 │ 4.12 │ 4120 (0) │
│ Precond apply │ 0.0 │ 0.0 │ 0 (0) │
╰────────────────┴────────────┴────────────────┴──────────╯
╭───────────────┬────────┬────────────┬─────────╮
│ Timing type │ Each │ Relative │ Total │
│ │ ms │ Percentage │ s │
├───────────────┼────────┼────────────┼─────────┤
│ Properties │ 0.0807 │ 2.70 % │ 0.3324 │
│ Equations │ 0.0809 │ 3.36 % │ 0.4144 │
│ Assembly │ 0.0345 │ 1.43 % │ 0.1767 │
│ Linear solve │ 2.6155 │ 87.43 % │ 10.7758 │
│ Linear setup │ 0.0000 │ 0.00 % │ 0.0000 │
│ Precond apply │ 0.0000 │ 0.00 % │ 0.0000 │
│ Update │ 0.0464 │ 1.55 % │ 0.1913 │
│ Convergence │ 0.0320 │ 1.33 % │ 0.1640 │
│ Input/Output │ 0.0592 │ 0.48 % │ 0.0592 │
│ Other │ 0.0514 │ 1.72 % │ 0.2116 │
├───────────────┼────────┼────────────┼─────────┤
│ Total │ 2.9916 │ 100.00 % │ 12.3253 │
╰───────────────┴────────┴────────────┴─────────╯Plot results
We plot the saturation front for the base case at different times together with the final solution for the refined model. In this case, refining the grid by a factor 10 gave us significantly less smearing of the trailing front.
using GLMakie
x = range(0, stop = 1, length = n)
x_f = range(0, stop = 1, length = n_f)
f = Figure()
ax = Axis(f[1, 1], ylabel = "Saturation", title = "Buckley-Leverett displacement")
for i in 1:6:length(states)
lines!(ax, x, states[i][:Saturations][1, :], color = :darkgray)
end
lines!(ax, x_f, states_refined[end][:Saturations][1, :], color = :red)
f
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 20.166483645 seconds to complete.This page was generated using Literate.jl.