examples.buckley_leverett

  1"""The Buckley--Leverett solution -- the one *exact* answer we can be checked against.
  2
  3Every other example here is an *illustration*: we plot what the code does, and
  4argue that it looks right. This one is a **verification**: in 1D, with constant
  5total velocity, the water equation has an analytic solution
  6(Buckley & Leverett, 1942), so the numerical profile can be compared with the
  7truth, and the error made to shrink under refinement.
  8
  9**The construction.** With $ ∇ ⋅ v = 0 $ the transport equation collapses to the
 10scalar conservation law
 11$$ φ \\, ∂s/∂t + v \\, ∂f(s)/∂x = 0 \\,, $$
 12whose characteristics carry a given $s$ at speed $ v f'(s)/φ $. Since $f$ is
 13S-shaped, $f'$ is *not* monotone, so those characteristics cross: the profile
 14would become multivalued, and a **shock** forms instead. Its saturation,
 15$S_f$, is fixed by requiring that the shock speed (from mass balance across it)
 16equal the characteristic speed just behind it -- which is the **Welge (1952)
 17tangent** construction: the chord from the initial state to $(S_f, f(S_f))$ must
 18be tangent to $f$,
 19$$ \\frac{f(S_f)}{S_f - s_\\mathrm{wc}} = f'(S_f) \\,. $$
 20Everything else follows: the shock travels at that chord's slope, breakthrough
 21occurs when it reaches the outlet, and thereafter the outlet saturation is read
 22off $ f'(s) = 1/t_D $ (which is Welge's production forecast).
 23
 24**Closed form for this model.** The relative permeabilities
 25(`minires.fluids.Fluid.RelPerm`) being quadratic, the tangent condition can
 26be solved by hand. In terms of the normalized saturation (ref
 27`minires.fluids.Fluid.rescale_sat`) and the endpoint mobility ratio $ M = v_o/v_w $,
 28$$ S_f^* = \\frac{1}{\\sqrt{1 + M}} \\,, \\qquad
 29   t_D^\\mathrm{bt} = (1 - s_\\mathrm{wc} - s_\\mathrm{or})
 30       \\, \\frac{2 \\, (1 + M - \\sqrt{1+M})}{M \\, \\sqrt{1+M}} \\,, $$
 31so for the default unit-viscosity fluids $ S_f = 1/\\sqrt{2} ≈ 0.7071 $ and
 32breakthrough comes at $ 2(\\sqrt{2}-1) ≈ 0.8284 $ pore volumes injected. Below,
 33the tangent is located *numerically* (from the model's own
 34`minires.fluids.Fluid.fractional_flow`) and asserted to agree with these -- so the check cuts both ways:
 35it validates the analytic solution we then compare the simulation against.
 36
 37Notes on the setup:
 38
 39- Time is measured in **pore volumes injected** ($t_D$), and distance in
 40  fractions of the length ($x_D$). The rate, pore volume and length are all
 41  $1$ here, so $ t_D = t $ and $ x_D = y $ -- no scaling clutters the plots.
 42- The 1D domain is a single **column** of cells (`Nx=1`), i.e. the flow is
 43  along $y$. A **row** (`Ny=1`) works just the same (`tests/test_transport.py`
 44  checks that the two agree); the column is merely what the plots below assume.
 45- The ends are **wells** (a source and a sink), not boundary conditions, which
 46  costs two $O(h)$ discrepancies with the textbook problem: the inlet cell only
 47  approaches $ s = 1 - s_\\mathrm{or} $ (it is filled at a finite rate, not
 48  held at a value), and the outlet is sampled half a cell short of $ x_D = 1 $.
 49  Both are visible below, and both vanish under refinement.
 50
 51In the figures:
 52
 53- "fractional flow": the construction itself. The tangent from the initial
 54  state finds $ S_f = 0.71 $ for the unit-viscosity fluids, but only $0.44$ for
 55  the case with $ M = 5 $ -- a *more adverse* mobility ratio gives a *weaker*
 56  shock, arriving sooner (at $ t_D = 0.35 $ rather than $0.83$), which is the
 57  whole reason mobility ratio matters to a waterflood. Where the tangent
 58  reaches $ f = 1 $ (open squares) is Welge's second reading of it: the
 59  *average* saturation behind the front -- and hence, nothing having been
 60  produced yet, the recovery at breakthrough. Asserted against the simulation
 61  below.
 62- "saturation profile": the numerical and analytic profiles, at a time before
 63  breakthrough. The explicit scheme follows the rarefaction to within 0.017 in
 64  saturation (0.005 on average), and smears the shock over some 4 cells. The
 65  implicit scheme is four times as far off in the rarefaction, and spreads the
 66  shock over five times as many cells -- as in `examples.quarter_five_spot`, it is
 67  the more diffusive of the two. Note that neither *overshoots*.
 68- "verification" (left): the water cut at the producer, against Welge's
 69  forecast. It breaks through a shade early (the producer being half a cell
 70  inside), then follows the analytic curve to within 0.004.
 71- "verification" (right): the $L_1$ error of the profile, which shrinks as
 72  $ h^{0.87} $ -- first-order convergence, all but the last 13% of it, the
 73  shortfall being the usual one for a monotone first-order scheme resolving a
 74  discontinuity. (The zig-zag about that fit is also expected: the error
 75  depends on where the shock happens to fall within a cell.) The rate is what
 76  matters: it certifies that the scheme is consistent, i.e. that the error is
 77  *discretization*, not *bug*.
 78"""
 79
 80from mpl_tools.place import freshfig
 81import numpy as np
 82from scipy.optimize import minimize_scalar
 83
 84from minires import ResSim
 85from minires.plotting import show
 86
 87## Setup
 88
 89
 90def make_model(N: int, fluid: dict) -> ResSim:
 91    """A 1D column of `N` cells: injector at the bottom, producer at the top.
 92
 93    Unit length, unit pore volume, unit rate -- so that time *is* $t_D$
 94    (pore volumes injected), and position *is* $x_D$.
 95    """
 96    return ResSim(Lx=1, Ly=1, Nx=1, Ny=N, fluid=fluid,
 97                  wells=[dict(xy=[0, 0], rate=+1),
 98                         dict(xy=[0, 1], rate=-1)])
 99
100
101## The analytic solution
102# NB: built from the model's *own* `fluid.fractional_flow`, so that it cannot
103# drift from the simulator's notion of the fluids -- only from its numerics.
104
105
106def welge_tangent(model) -> tuple:
107    """Locate the shock: `(S_f, shock speed)`, ref the docstring.
108
109    Rather than solving $ f(S)/(S - s_\\mathrm{wc}) = f'(S) $, we *maximize*
110    the chord slope -- the same point, but needing no derivative, and with no
111    root-bracketing to get wrong.
112    """
113    lo, hi = model.fluid.swc, 1 - model.fluid.sor
114    chord = lambda S: model.fluid.fractional_flow(S) / (S - lo)  # noqa: E731
115    opt = minimize_scalar(lambda S: -chord(S), method="bounded",
116                          bounds=(lo + 1e-9, hi), options=dict(xatol=1e-12))
117    return opt.x, chord(opt.x)
118
119
120def analytic(model, tD, xD, S_f):
121    """The saturation at `xD` at time `tD`: the profile, and its shock.
122
123    Behind the shock, $ x_D = t_D \\, f'(s) $ -- monotone in $s$ over
124    $ [S_f, 1 - s_\\mathrm{or}] $ (the tangent point lying beyond $f$'s
125    inflection), so it is inverted by interpolation. The shock itself needs no
126    special treatment: the tangent condition puts its position at exactly the
127    end of that range, so anything ahead of it is simply `right=swc`.
128    """
129    ss = np.linspace(S_f, 1 - model.fluid.sor, 10001)
130    xx = tD * model.fluid.dfractional_flow(ss)  # speed of each s
131    return np.interp(xD, xx[::-1], ss[::-1], right=model.fluid.swc)
132
133
134## Verify the analytic solution against the closed form
135cases: dict = dict(
136    A=dict(),                              # defaults: vw = vo = 1, swc = sor = 0
137    B=dict(vo=5., swc=.2, sor=.2),         # a contrast in both viscosity and endpoints
138)
139
140models = {case: make_model(200, fluid) for case, fluid in cases.items()}
141
142for case, model in models.items():
143    S_f, speed = welge_tangent(model)
144    # The closed form (ref the docstring), in terms of `M` and the endpoints
145    M = model.fluid.vo / model.fluid.vw
146    span = 1 - model.fluid.swc - model.fluid.sor
147    S_f_exact = model.fluid.swc + span / np.sqrt(1 + M)
148    tD_bt_exact = span * 2*(1 + M - np.sqrt(1 + M)) / (M * np.sqrt(1 + M))
149    assert np.isclose(S_f, S_f_exact, rtol=1e-8), "Welge tangent misplaced."
150    assert np.isclose(1/speed, tD_bt_exact, rtol=1e-8), "Closed form disagrees."
151    print(f"Case {case}: M = {M:g}, S_f = {S_f:.4f}, breakthrough at {1/speed:.4f} PVI")
152
153## Simulate: the profile, both schemes, both cases
154tD_snap = dict(A=.5, B=.3)  # a time before breakthrough, for each case
155profiles: dict = {}
156
157for case, model in models.items():
158    S_f, speed = welge_tangent(model)
159    xD = model.mesh[1].ravel()
160    S0 = np.full(model.Nxy, model.fluid.swc)
161
162    nSteps = 50
163    dt = tD_snap[case] / nSteps
164    S_exp, _ = model.sim(dt, nSteps, S0, pbar=False)
165    S_imp, _ = model.sim(dt, nSteps, S0, pbar=False, implicit=True)
166
167    profiles[case] = dict(xD=xD, S_f=S_f, speed=speed,
168                          exact=analytic(model, tD_snap[case], xD, S_f),
169                          explicit=S_exp[-1], implicit=S_imp[-1])
170
171    # Neither scheme may overshoot the physical range: the analytic solution is
172    # bounded by its data, and a monotone scheme must be too.
173    for scheme, S in [("explicit", S_exp), ("implicit", S_imp)]:
174        assert model.fluid.swc - 1e-12 <= S.min() and S.max() <= 1 - model.fluid.sor + 1e-12, (
175            f"Case {case}, {scheme} scheme: saturation out of bounds.")
176
177## Simulate: the production history (case A, past breakthrough)
178model = models["A"]
179S_f, speed = welge_tangent(model)
180tD_bt = 1 / speed
181
182nSteps = 150
183dt = 1.5 / nSteps
184tt = dt * np.arange(nSteps + 1)
185SS, _ = model.sim(dt, nSteps, np.full(model.Nxy, model.fluid.swc), pbar=False)
186
187# The water cut is the fractional flow of the producer's cell -- ref
188# `ResSim.assemble_wells`, which is what draws the produced fluid at that ratio.
189i_prd = model.xy2ind(*model.wells.xy[1])
190water_cut = model.fluid.fractional_flow(SS[:, i_prd])
191# Welge's forecast: the outlet saturation is the one whose characteristic has
192# just arrived. `analytic` returns `swc` (whence a zero water cut) before that.
193water_cut_exact = np.array([model.fluid.fractional_flow(analytic(model, t, 1., S_f))
194                            for t in tt])
195water_cut_exact[0] = 0  # `tD = 0` puts the whole profile at the inlet
196
197# Mass balance: what was injected is either still in place, or was produced.
198# (The tolerance is set by the trapezoidal integration of the jump at
199# breakthrough, not by the scheme, which conserves mass exactly.)
200in_place = (SS[-1] - model.fluid.swc).mean()
201produced = np.trapezoid(water_cut, tt)
202assert np.isclose(in_place + produced, tt[-1], rtol=2e-3), "Water unaccounted for."
203
204# Before breakthrough that balance is *exact*, nothing having been produced.
205# NB: "nothing" is not quite `0`: the explicit scheme's stencil advances one
206# cell per sub-step, so a (multiplicatively vanishing) tail of the front runs
207# ahead of it -- reaching the producer at some $ 10^{-150} $, long before the
208# water does. Hence the threshold, which also defines breakthrough below.
209DRY = 1e-4  # water cut counting as "no water", i.e. $ s ⪅ 0.01 $
210k_pre = round(.7 / dt)
211assert water_cut[k_pre] < DRY, "Breakthrough far too early."
212assert np.isclose((SS[k_pre] - model.fluid.swc).mean(), tt[k_pre]), "Injected water lost."
213
214# Breakthrough should be *slightly* early, the producer sitting half a cell
215# short of the outlet -- i.e. by `(hy/2) / speed`, which is under one `dt` here.
216tD_bt_sim = tt[(water_cut > DRY).argmax()]
217assert 0 <= tD_bt - tD_bt_sim < 2*dt, "Breakthrough mistimed."
218
219# Evaluated at breakthrough, it is the *other* reading of the tangent (ref the
220# "fractional flow" figure): the mean saturation is then the average behind the
221# front, $ s_wc + t_D^bt $ -- to within the smearing.
222k_bt = round(tD_bt / dt)
223assert np.isclose(SS[k_bt].mean(), model.fluid.swc + tD_bt, rtol=1e-2), (
224    "Welge average is off.")
225
226## Simulate: convergence under refinement (case A, explicit scheme)
227NN = np.array([50, 100, 200, 400, 800])
228L1 = np.zeros(len(NN))
229
230for i, N in enumerate(NN):
231    m = make_model(N, cases["A"])
232    S, _ = m.sim(tD_snap["A"]/50, 50, np.full(m.Nxy, m.fluid.swc), pbar=False)
233    # NB: the analytic solution is grid-independent -- only sampled anew
234    exact = analytic(m, tD_snap["A"], m.mesh[1].ravel(), profiles["A"]["S_f"])
235    L1[i] = abs(S[-1] - exact).mean()
236
237fit = np.polyfit(np.log(NN), np.log(L1), 1)
238rate = -fit[0]
239print(f"Convergence: L1 error ~ h^{rate:.2f}")
240assert .7 < rate < 1.1, "Lost (near-)first-order convergence."
241
242## Plot: the fractional-flow curve and the Welge tangent
243fig, ax = freshfig("Buckley-Leverett -- fractional flow", figsize=(6, 5))
244
245for case, p in profiles.items():
246    m = models[case]
247    S_f, tD_bt_ = p["S_f"], 1 / p["speed"]
248    ss = np.linspace(m.fluid.swc, 1 - m.fluid.sor, 201)
249    (h,) = ax.plot(ss, m.fluid.fractional_flow(ss), lw=2,
250                   label=f"$M$ = {m.fluid.vo/m.fluid.vw:g}, "
251                         f"$s_\\mathrm{{wc}}$ = {m.fluid.swc:g}, "
252                         f"$s_\\mathrm{{or}}$ = {m.fluid.sor:g}")
253    # The tangent, from the initial state up to `f = 1`, which it reaches at
254    # the *average* saturation behind the front (Welge's other reading of it).
255    ax.plot([m.fluid.swc, m.fluid.swc + tD_bt_], [0, 1], ":", c=h.get_color(), lw=1)
256    ax.plot(S_f, m.fluid.fractional_flow(S_f), "o", c=h.get_color(), ms=8,
257            label=f"$S_f$ = {S_f:.3f},  $t_D^\\mathrm{{bt}}$ = {tD_bt_:.3f}")
258    ax.plot(m.fluid.swc + tD_bt_, 1, "s", c=h.get_color(), ms=6, mfc="none",
259            label=f"$\\bar{{s}}$ = {m.fluid.swc + tD_bt_:.3f} (at breakthrough)")
260
261ax.set(title="The Welge tangent construction", xlabel="Water saturation, $s$",
262       ylabel="Fractional flow, $f(s)$", xlim=(0, 1), ylim=(0, 1.08))
263ax.legend(fontsize="small", loc="lower right")
264fig.tight_layout()
265
266## Plot: the saturation profiles, numerical vs. analytic
267fig, axs = freshfig("Buckley-Leverett -- saturation profile", ncols=2,
268                    sharey=True, figsize=(10, 4.5))
269
270for ax, (case, p) in zip(axs, profiles.items()):
271    ax.plot(p["xD"], p["exact"], "k-", lw=2, label="Analytic (Buckley-Leverett)")
272    ax.plot(p["xD"], p["explicit"], "C0.", ms=4, label="Explicit (upwind)")
273    ax.plot(p["xD"], p["implicit"], "C1.", ms=4, label="Implicit (Newton)")
274    ax.axvline(tD_snap[case] * p["speed"], c="k", ls=":", lw=1,
275               label="Shock position, $t_D \\, f'(S_f)$")
276    m = models[case]
277    ax.set(title=f"Case {case}:  $M$ = {m.fluid.vo/m.fluid.vw:g},"
278                 f"  $t_D$ = {tD_snap[case]}", xlabel="$x_D$")
279axs[0].set_ylabel("Water saturation, $s$")
280axs[0].legend(fontsize="small")
281fig.tight_layout()
282
283## Plot: water cut, and the convergence of the profile
284fig, (ax1, ax2) = freshfig("Buckley-Leverett -- verification", ncols=2,
285                           figsize=(10, 4.5))
286
287ax1.plot(tt, water_cut_exact, "k-", lw=2, label="Welge forecast")
288ax1.plot(tt, water_cut, "C0.", ms=4, label="Simulated (explicit)")
289ax1.axvline(tD_bt, c="k", ls=":", lw=1,
290            label=f"Breakthrough, $2(\\sqrt{{2}}-1)$ = {tD_bt:.4f}")
291ax1.set(title="Water cut at the producer", xlabel="$t_D$ (pore volumes injected)",
292        ylabel="$f_w$", ylim=(-.03, 1))
293ax1.legend(fontsize="small", loc="lower right")
294
295ax2.loglog(NN, L1, "C0-o", label="$L_1$ error")
296ax2.loglog(NN, np.exp(np.polyval(fit, np.log(NN))), "C0--", lw=1,
297           label=f"Fit: $\\propto h^{{{rate:.2f}}}$")
298ax2.loglog(NN, L1[0] * NN[0]/NN, "k:", lw=1, label="$O(h)$, for reference")
299ax2.set(title=f"Convergence of the profile at $t_D$ = {tD_snap['A']}",
300        xlabel="$N_y$", ylabel="Mean $|s - s_\\mathrm{exact}|$",
301        xticks=NN, xticklabels=[str(N) for N in NN])
302ax2.minorticks_off()
303ax2.legend(fontsize="small")
304fig.tight_layout()
305
306# Regression values, checked by `tests/test_examples.py`.
307__digest__ = dict(explicit  = profiles["A"]["explicit"],
308                  implicit  = profiles["A"]["implicit"],
309                  case_B    = profiles["B"]["explicit"],
310                  welge     = [profiles[c][k] for c in "AB"
311                               for k in ["S_f", "speed"]],
312                  water_cut = water_cut[water_cut > DRY],
313                  bt        = [tD_bt_sim, tD_bt],
314                  L1        = L1)
315
316if __name__ == "__main__":
317    show()