examples.egg

  1"""The Egg model, flattened to 2D, against its published 3D solution. **In metric units.**
  2
  3The [Egg model](https://doi.org/10.1002/gdj3.21) (Jansen et al., 2014) is a synthetic
  4benchmark reservoir, much used in history-matching and optimization studies: an
  5egg-shaped region of 60 × 60 × 7 cells (8 × 8 × 4 m, 18553 of them active) of
  6**channelized** permeability (50-7000 mD; an ensemble of 101 realizations, of which
  7this is the first) and porosity 0.2, flooded for 3600 days by 8 water injectors at
  879.5 m³/day into 4 producers held at 395 bar, from an initial 400 bar. Oil 5 cP,
  9water 1 cP, both at 1e-5/bar; Corey relative permeabilities with exponents 3 and 4.
 10The data (`egg.npz`, ref the "Data" section below) are JutulDarcy's copy of the
 11ECLIPSE deck, which comes with the deck's solution by ECLIPSE 100 -- in 3D, with
 12gravity -- as its reference.
 13
 14## The 2D reduction
 15
 16This simulator is areal, so the 7 layers are **vertically averaged**: for flow along
 17the layers, the exact upscaling of a stack is the arithmetic mean of their
 18permeabilities (each carries flux in proportion to its own), the inactive ones
 19counting as 0 over the 28 m column, and the porosity likewise; the egg's outline
 20becomes the `ResSim.active` mask. What is lost is
 21gravity (a mere 0.27 bar of hydrostatic contrast over the column, against a 5 bar
 22drawdown) and the layers' *individual* channel patterns, through the fastest of
 23which the water arrives first in 3D. The wells, being vertical through all 7
 24layers, are just one completion each here; the rate, per unit thickness (ref
 25`ResSim.cdarcy`), is the 79.5 m³/day over the 28 m.
 26
 27The deck's own quirks are kept: the water starts at $S_w = 0.1$, *below* the 0.2 at
 28which it becomes mobile (so the first tenth of a pore volume injected only fills up
 29the immobile saturation), and the injectors' 420 bar limit, which never binds. The
 30relative permeabilities are the deck's Corey curves (`minires.fluids.Fluid`): exponents 3
 31and 4 with end-points 0.6 and 0.8 at the residuals $ S_{wc} = 0.2 $, $ S_{or} = 0.15 $
 32(the table runs on to $ S_w = 0.9 $, but nothing but water moves there). Their
 33fractional flow is steeper than the quadratic default's (its maximal slope 5.65,
 34against 3.5), which the CFL estimate follows (`ResSim.estimate_1CFL`).
 35
 36## Validation
 37
 38The 2D model reproduces the 3D reference to within an RMS of **0.01** in the water
 39cut of every producer, and to within 4-6% on average in their oil rates; the
 40recovery after 3600 days is 0.57 of the oil in place, against 0.59. Breakthrough
 41comes 15-30 days later than in 3D at three of the four producers -- the vertical
 42averaging smoothing the fastest layer's channel away, as said -- and PROD2, the
 43strongest well, produces 7% less oil.
 44
 45That a 2D areal model gets this close is because the reservoir is thin (28 m
 46against 480 m across), the wells penetrate all of it, and the channels are largely
 47stacked, so that the flow is nearly layer-parallel. It is *not* generic: had the
 48layers had separate channel systems, the average would have connected wells that
 49none of the layers connects.
 50
 51The same flattened model run in JutulDarcy (a fully implicit code, reading the
 52reduction as an ECLIPSE deck) agrees with this one to within 0.007 RMS in the water
 53cut, and lies as far from the 3D original as this one does (0.008 against 0.009):
 54the residual is the flattening, not the simulator.
 55
 56In the figures:
 57
 58- "fields" (left): the vertically averaged permeability, its channels
 59  (the yellow bands) running roughly N-S; injectors are black triangles, producers
 60  white. Outside the egg the cells are inactive (`ResSim.active`), hence blank.
 61- "fields" (right three): the oil being displaced -- fastest along the channels,
 62  and first towards PROD2 (the producer sitting on the widest one, near the centre),
 63  which breaks through after 6 months.
 64- "production" (left): the water cut, this model (solid) against the 3D reference
 65  (dashed), per producer. The curves lie on top of each other at the scale of the
 66  plot; the 2D breakthroughs come slightly later.
 67- "production" (right): the oil rate, on a log axis, likewise; the 3D reference's
 68  early points are its finer initial report steps. PROD2's initial rate is set by
 69  its well index -- the same Peaceman formula in both codes -- and its later
 70  decline by how fast the water arrives.
 71
 72## Data
 73
 74`egg.npz` (73 KB) holds, `x`-first as the simulator orders its fields: `perm`, the
 75deck's `PERMX` in mD, `(60, 60, 7)` = `(Nx, Ny, Nz)`, `float32` (`PERMY = PERMX` and
 76`PERMZ = 0.1 PERMX` in the deck, irrelevant in 2D); `active`, its `ACTNUM`, `bool`;
 77and the reference solution's 135 report times (`ref_time`, in days) with its oil and
 78water production rates (`ref_orat`, `ref_wrat`, `(135, 4)`, m³/day, `PROD1..4`).
 79Everything else the model needs -- cell size, porosity, fluids, wells, controls -- is
 80a handful of numbers, kept in this script.
 81
 82The source is the deck as JutulDarcy ships it (`GeoEnergyIO`'s `EGG` test input:
 83`EGG.DATA` with `MDARCY.INC`, `ACTIVE.INC`, and `REFERENCE.CSV`, the ECLIPSE 100
 84solution in SI, production signed negative), i.e. realization 1 of the 101 published
 85by TU Delft (doi:10.4121/uuid:916c86cd-3558-4672-829a-105c62985ab2, 4TU General
 86Terms of Use). With `d` that directory, the file was made by
 87
 88    import re, numpy as np
 89    def keyword(path, kw):                  # a GRID keyword's numbers, as (Nx, Ny, Nz)
 90        body = re.sub(r"--.*", "", open(path).read().split(kw, 1)[1].split("/", 1)[0])
 91        vals = [float(v) for tok in body.split() for n, _, v in [tok.rpartition("*")]
 92                for _ in range(int(n) if n else 1)]              # ECLIPSE's `n*value`
 93        return np.array(vals).reshape(7, 60, 60).transpose(2, 1, 0)   # I fastest
 94    perm = keyword(f"{d}/MDARCY.INC", "PERMX")
 95    active = keyword(f"{d}/ACTIVE.INC", "ACTNUM").astype(bool)
 96    ref = np.genfromtxt(f"{d}/REFERENCE.CSV", delimiter=",", names=True)
 97    prods = [f"PROD{i}" for i in range(1, 5)]
 98    orat = -np.column_stack([ref[w + "orat"] for w in prods]) * 86400
 99    wrat = -np.column_stack([ref[w + "wrat"] for w in prods]) * 86400
100    wrat[np.abs(wrat) < 1e-9] = 0                               # the reference's round-off
101    np.savez_compressed("egg.npz", perm=perm.astype(np.float32), active=active,
102                        ref_time=ref["time"] / 86400, ref_orat=orat, ref_wrat=wrat)
103"""
104
105from pathlib import Path
106
107from mpl_tools.place import freshfig
108import numpy as np
109
110from minires import ResSim
111from minires.plotting import show
112
113## The data: realization 1 of the Egg ensemble (ref the module docstring)
114data = np.load(Path(__file__).with_name("egg.npz"))
115perm, active = data["perm"], data["active"]       # (Nx, Ny, Nz) = (60, 60, 7)
116Nx, Ny, Nz = perm.shape
117h, dz = 8.0, 4.0                                  # cell size [m]
118H = Nz*dz                                         # thickness [m]
119
120## The 2D reduction: vertical averaging
121# For flow along the layers, the exact upscaling of a stack of layers is the
122# arithmetic mean of their permeabilities (each carries flux in proportion to
123# its own), with the inactive layers (`K = 0`) counted in the thickness. The pore
124# volume is likewise the column's. The egg-shaped footprint -- the cells active in
125# any layer -- becomes the model's `active` mask; outside it, `K` and `por` are
126# immaterial (but `K` must be finite and positive, so the mean is filled in).
127K   = (perm * active).sum(-1) / Nz
128por = 0.2 * active.mean(-1)
129footprint = active.any(-1)
130K   = np.where(footprint, K, K[footprint].mean())
131
132## Wells: 8 injectors on rate, 4 producers on BHP, from the deck's `WELSPECS`
133# (1-based cell indices) and `SCHEDULE`. The rate is per unit thickness (ref
134# `ResSim.cdarcy`), the whole column being one cell here.
135ij = dict(INJECT1=(5, 57), INJECT2=(30, 53), INJECT3=(2, 35), INJECT4=(27, 29),
136          INJECT5=(50, 35), INJECT6=(8, 9), INJECT7=(32, 2), INJECT8=(57, 6),
137          PROD1=(16, 43), PROD2=(35, 40), PROD3=(23, 16), PROD4=(43, 18))
138q_inj, p_prod, p_max, rw = 79.5/H, 395, 420, 0.2   # m²/day, bar, bar, m
139wells = {name: dict(xy=[(i - .5)*h, (j - .5)*h], rw=rw,
140                    **(dict(rate=q_inj) if name.startswith("INJ") else dict(bhp=p_prod)))
141         for name, (i, j) in ij.items()}
142producers = [i for i, name in enumerate(ij) if name.startswith("PROD")]
143
144## Fluids: 1 / 5 cP, both 1e-5/bar (rock incompressible), and the deck's `SWOF`
145C = 86400 * 9.869233e-16 * 1e5 / 1e-3               # m, day, bar, mD, cP: 0.008527
146Sw0 = 0.1                                           # initial water saturation
147p0 = 400                                            # initial pressure [bar]
148
149
150# The deck's `SWOF` table is Corey (ref `minires.fluids.Fluid`): exponents 3 (water) and 4
151# (oil), end-points 0.6 and 0.8, water immobile below Sw = 0.2 and oil below So = 0.15.
152fluid = dict(vw=1, vo=5, swc=0.2, sor=0.15, nw=3, no=4, krw0=0.6, kro0=0.8)
153
154model = ResSim(Lx=Nx*h, Ly=Ny*h, Nx=Nx, Ny=Ny, cdarcy=C, K=K, por=por, active=footprint,
155               ct=1e-5, fluid=fluid, wells=wells)
156
157## Simulate 3600 days in 30-day steps
158dt, nSteps = 30, 120
159tt = dt * np.arange(1, nSteps + 1)
160SS, PP = model.sim(dt, nSteps, np.full(model.Nxy, Sw0), P0=np.full(model.Nxy, p0), pbar=False)
161
162# The injectors' 420 bar limit (which the deck also has) never binds, so rate
163# control is the whole story, as in the reference.
164bhp_inj = model.wells.actual_bhp[:8]
165assert bhp_inj.max() < p_max, "the injectors' BHP limit binds"
166
167## Production, per well: water cut and oil rate
168cells = model.xy2ind(*model.wells.xy[producers].T)
169cut = model.fluid.fractional_flow(SS[1:][:, cells])   # (nSteps, 4)
170rate = -model.wells.actual_rates[producers].T * H        # m³/day, total, positive
171oil = rate * (1 - cut)
172
173## The reference: the 3D model (ECLIPSE 100, via JutulDarcy's copy of the deck)
174ref_t, ref_o, ref_w = data["ref_time"], data["ref_orat"], data["ref_wrat"]
175ref_cut = ref_w / (ref_w + ref_o)
176sel = ref_t >= dt                                        # our first report time
177cut_at_ref = np.column_stack([np.interp(ref_t[sel], tt, cut[:, j]) for j in range(4)])
178oil_at_ref = np.column_stack([np.interp(ref_t[sel], tt, oil[:, j]) for j in range(4)])
179rms_cut = np.sqrt(np.mean((cut_at_ref - ref_cut[sel])**2, axis=0))
180rel_oil = np.abs(oil_at_ref - ref_o[sel]).mean(axis=0) / ref_o[sel].mean(axis=0)
181ooip = 0.2 * h*h*dz * active.sum() * (1 - Sw0)                      # m³
182recovery = oil.sum()*dt / ooip
183recovery_ref = np.trapezoid(ref_o, ref_t, axis=0).sum() / ooip
184assert (rms_cut < 0.02).all() and (rel_oil < 0.1).all(), "the 2D model has drifted from the 3D reference"
185
186print(f"RMS water-cut misfit per producer: {rms_cut.round(3)}")
187print(f"Mean relative oil-rate misfit per producer: {rel_oil.round(3)}")
188print(f"Recovery (fraction of OOIP): {recovery:.3f}, 3D reference {recovery_ref:.3f}")
189print(f"Injector BHP: {bhp_inj.min():.1f}-{bhp_inj.max():.1f} bar")
190
191## Plot: the permeability, and the water's advance
192fig, axs = freshfig("Egg -- fields", ncols=4, figsize=(15, 4), sharex=True, sharey=True)
193kws: dict = dict(colorbar=False, finalize=False, wells=dict(size=.4, text=False))
194model.plt_field(axs[0], np.log10(np.where(footprint, K, np.nan)).ravel(), cmap="viridis",
195                levels=np.linspace(2, 3.7, 18), title="$\\log_{10} K$ [mD]", **kws)
196for ax, k in zip(axs[1:], [20, 60, nSteps]):
197    model.plt_field(ax, SS[k], "oil", title=f"Oil saturation, t = {k*dt} days", labels=False, **kws)
198fig.tight_layout()
199
200## Plot: production, against the 3D reference
201fig, (ax1, ax2) = freshfig("Egg -- production", ncols=2, figsize=(12, 4))
202names = [list(ij)[i] for i in producers]
203for j, name in enumerate(names):
204    ax1.plot(tt, cut[:, j], c=f"C{j}", label=name)
205    ax1.plot(ref_t, ref_cut[:, j], "--", c=f"C{j}", lw=1)
206    ax2.plot(tt, oil[:, j], c=f"C{j}", label=name)
207    ax2.plot(ref_t, ref_o[:, j], "--", c=f"C{j}", lw=1)
208ax1.plot([], [], "k-", label="2D (this)")
209ax1.plot([], [], "k--", lw=1, label="3D (ECLIPSE 100)")
210ax1.set(title="Water cut", xlabel="Time [day]", ylabel="$f_w$")
211ax1.legend(fontsize="small", ncols=2)
212ax2.set(title="Oil rate", xlabel="Time [day]", ylabel="[m³/day]", yscale="log")
213fig.tight_layout()
214
215# Regression values, checked by `tests/test_examples.py`.
216__digest__ = dict(water_cut = cut,
217                  oil_rate  = oil,
218                  S_final   = SS[-1],
219                  bhp_inj   = bhp_inj[:, -1],
220                  rms_cut   = rms_cut,
221                  rel_oil   = rel_oil,
222                  recovery  = [recovery, recovery_ref])
223
224if __name__ == "__main__":
225    show()