<jupyter_start><jupyter_text>Well placement optimization in SPE 10 reservoir modelThis is my experimental notebook to conduct optimization of well placement for an inverted quarter-five (Q5) waterflooding in SPE 10 reservoir model. The MATLAB script of MRST is ported with Python in Google Colab, working in this way; an optimizer in Python optimizes a function that carries the whole MATLAB code and returns objective function.The objective function is to minimize OIP (=OOIP-cumulative production). Example of optimizer here is Bayesian optimization.   Install Octave in Colab<jupyter_code>!apt install octave<jupyter_output>Reading package lists... Done
Building dependency tree       
Reading state information... Done
The following additional packages will be installed:
  aglfn epstool fonts-droid-fallback fonts-noto-mono ghostscript gnuplot-data
  gnuplot-qt gsfonts imagemagick-6-common info install-info libamd2
  libauthen-sasl-perl libcamd2 libccolamd2 libcholmod3 libcolamd2
  libcupsfilters1 libcupsimage2 libcxsparse3 libdata-dump-perl libemf1
  libencode-locale-perl libfftw3-single3 libfile-listing-perl libfltk-gl1.3
  libfltk1.3 libfont-afm-perl libgail-common libgail18 libglpk40
  libgraphicsmagick++-q16-12 libgraphicsmagick-q16-3 libgs9 libgs9-common
  libgtk2.0-0 libgtk2.0-bin libgtk2.0-common libhtml-form-perl
  libhtml-format-perl libhtml-parser-perl libhtml-tagset-perl
  libhtml-tree-perl libhttp-cookies-perl libhttp-daemon-perl libhttp-date-perl
  libhttp-message-perl libhttp-negotiate-perl libijs-0.35 libio-html-perl
  libio-socket-ssl-perl libjbig2dec0 liblqr-1-0 liblua5.3-0
  liblwp-mediat[...]<jupyter_text>Cloning repositories<jupyter_code>!git clone https://github.com/yohanesnuwara/MRST-simulations
!git clone https://github.com/yohanesnuwara/reservoir_datasets
!git clone https://bitbucket.org/mrst/mrst-core.git
!git clone https://yohanesnuwara@bitbucket.org/mrst/mrst-autodiff.git<jupyter_output>Cloning into 'mrst-autodiff'...
remote: Counting objects: 40863, done.[K
remote: Compressing objects: 100% (12732/12732), done.[K
remote: Total 40863 (delta 29924), reused 37400 (delta 27804)[K
Receiving objects: 100% (40863/40863), 15.26 MiB | 8.31 MiB/s, done.
Resolving deltas: 100% (29924/29924), done.<jupyter_text>Main MATLAB program<jupyter_code>%%writefile addpaths.m
addpath '/content/mrst-core/utils/units'
addpath '/content/mrst-core/gridprocessing'
addpath '/content/mrst-core/utils'
addpath '/content/mrst-core/modules'
addpath '/content/mrst-core/solvers'
addpath '/content/mrst-core/utils/gridtools'
addpath '/content/mrst-core/params/rock'
addpath '/content/mrst-core/params/wells_and_bc'
addpath '/content/mrst-core/params'
addpath '/content/MRST-simulations/modules/incomp/fluid/incompressible'
addpath '/content/MRST-simulations/modules/incomp'
addpath '/content/MRST-simulations/modules/incomp/fluid/utils'
addpath '/content/MRST-simulations/modules/incomp/transport'
%%writefile SPE10_waterflooding.m
addpaths

%% Set up model and parameters
T           = 5*year;
nlayer      = 1; % number of layers of SPE 10 model to be simulated
startlayer  = 4; % the SPE 10 model starts from layer i-th defined
dims        = [60 220 nlayer];
domain      = dims.*[20 10 2]*ft;
G           = computeGeometry(cartGrid(dims,domain));

rock        = getSPE10rock((1:dims(1)), (1:dims(2)), ... 
              startlayer:startlayer+(nlayer-1));
            
rock.poro = max(rock.poro,.0005);

hT = computeTrans(G, rock); % compute transmissibility 
%x0 = initState(G,W,100*barsa, [0 1]); % initialize

%% Set fluid model
fluid = initSimpleFluid('mu' , [   1,    1] .* centi*poise     , ...
                        'rho', [1000,  850] .* kilogram/meter^3, ...
                        'n'  , [   2,    2]);

%% Set parameters describing the well
wtype    = {'bhp', 'bhp', 'bhp', 'bhp', 'bhp'};
wtarget  = [200,   200,   200,   200,   500] .* barsa();
wrad     = [0.125, 0.125, 0.125, 0.125, 0.125] .* meter;
wloc     = [  1,   60,     1,   60,  30;
              1,    1,   220,  220, 111];
wname    = {'P1', 'P2', 'P3', 'P4', 'I'};

W = [];
for w = 1 : numel(wtype)
  if w < 5
    % production well is oil
    phase = [0 1];
  else
    % production well is water
    phase = [1 0];
  endif
  W = verticalWell(W, G, rock, wloc(1,w), wloc(2,w), 1, ...
      'Type', wtype{w}, 'Val', wtarget(w), ...
      'Radius', wrad(w), 'Name', wname{w}, ...
      'InnerProduct', 'ip_tpf', 'Comp_i', phase);
end

x0 = initState(G,W,100*barsa, [0 1]); % initialize

%% Simulation
M = 50;
x = x0;
[dt,dT] = deal(zeros(M,1), T/M);
wellSol = cell(M,3);
oip     = zeros(M,3);
n = 1;

for i=1:M
    x  = incompTPFA(x, G, hT, fluid, 'wells', W);
    x  = implicitTransport(x, G, dT, rock, fluid, 'wells', W);    
        
    dt(i) = dT;
    oip(i,n) = sum(x.s(:,2).*poreVolume(G,rock));
    wellSol{i,n} = getWellSol(W,x, fluid);
end

%% Print the Oil in Place
display(oip(50));<jupyter_output>Overwriting SPE10_waterflooding.m<jupyter_text>Functions<jupyter_code>%%writefile /content/mrst-core/modules/getSPE10rock.m
function rock = getSPE10rock(varargin)
%Define rock properties for Model 2 of tenth SPE CSP
%
% SYNOPSIS:
%   rock = getSPE10rock
%   rock = getSPE10rock(layers)
%   rock = getSPE10rock(I, J, K)
%
% PARAMETERS:
%   layers - Which of the 85 model layers to include in a specific test.
%            OPTIONAL.  If unspecified, all 85 layers (for a total of
%            60-by-220-by-85 == 1,122,000 cells) are included.
%
%            Some possible values are
%               layers = ( 1 : 35).';  %  Tarbert formation
%               layers = (36 : 85).';  %  Upper Ness formation
%
%   I,J,K  - Global Cartesian indices identifying subset of rock data.
%            Useful, e.g., in extracting lateral sections.  Arrays must
%            satisfy
%
%               ALL((0 < I) & (I <=  60)) && ...
%               ALL((0 < J) & (J <= 220)) && ...
%               ALL((0 < K) & (K <=  85))
%
%            OPTIONAL.  If unspecified, treated as I=1:60,J=1:220,K=1:85
%            (i.e., the entire model; 1,122,000 cells).  In other words,
%            equivalently to an unspecified 'layers' input.
%
% RETURNS:
%   rock - Rock structure having fields 'perm' and 'poro' pertaining to
%          the specified layers or cell subset.
%
% NOTE:
%   The permeability data is returned in strict SI units (metres squared).
%
% EXAMPLE:
%   rock = getSPE10rock(85)
%
% SEE ALSO:
%   `getSPE10setup`, `make_spe10_data`.

%{
Copyright 2009-2019 SINTEF Digital, Mathematics & Cybernetics.

This file is part of The MATLAB Reservoir Simulation Toolbox (MRST).

MRST is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.

MRST is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
GNU General Public License for more details.

You should have received a copy of the GNU General Public License
along with MRST.  If not, see <http://www.gnu.org/licenses/>.
%}

   assert (all(cellfun(@isnumeric, varargin)), ...
           'Input arrays must be numeric.');

   % Load data in pre-processed form.
   rock = load_mat_file();

   % Extract caller's requested subset from 'rock' data
   %
   % Return only PERM and PORO data and exclude any other information that
   % might be stored in on-disk representation of rock data.
   %
   ix   = define_subset(varargin{:});

   rock = struct('perm', rock.perm(ix, :), ...
                 'poro', rock.poro(ix   ));
end

%--------------------------------------------------------------------------

function ix = define_subset(varargin)
   [Nx, Ny, Nz] = deal(60, 220, 85);
                                          % ()
   [I, J, K]    = deal(1:Nx, 1:Ny, 1:Nz); % Default (entire dataset)

   if nargin == 1
                                          % (layers)
      K = varargin{1};                    % Caller specified ind. layers

   elseif nargin == 3
                                          % (I, J, K)
      [I, J, K] = deal(varargin{:});      % Caller specified box

   elseif nargin ~= 0
      file = mfilename();
      error(['Syntax is\n\t'              , ...
             'rock = %s         %% or\n\t', ...
             'rock = %s(layers) %% or\n\t', ...
             'rock = %s(I, J, K)'], file, file, file);
   end

   validate_range(I, Nx);
   validate_range(J, Ny);
   validate_range(K, Nz);

   [I, J, K] = ndgrid(I, J, K);

   ix = sub2ind([Nx, Ny, Nz], I(:), J(:), K(:));
end

%--------------------------------------------------------------------------

function rock = load_mat_file()
   %rdir = getDatasetPath('spe10', 'skipAvailableCheck', true);
   rdir = '/content/mrst-core/data/SPE10';
   rock_file = fullfile(rdir, 'spe10_rock');

   if ~exist([rock_file, '.mat'], 'file')
      ok = makeSPE10DataAvailable();

      if ~ok
         error('SPE10Download:Fail', 'Failed to download SPE 10 dataset');
      end
   end

   data = load(rock_file);
   rock = data.rock; % Fa-Fa-Fa
end

%--------------------------------------------------------------------------

function validate_range(i, n)
   assert (all((0 < i) & (i <= n)), ...
           '%s outside valid range 1:%d', inputname(1), n);
end
%%writefile /content/mrst-core/modules/getWellSol.m
function wellSol = getWellSol(W, x, fluid)

mu = fluid.properties();
wellSol(numel(W))=struct;
for i=1:numel(W)
    out = min(x.wellSol(i).flux,0); iout = out<0;
    in  = max(x.wellSol(i).flux,0); iin  = in>0;
    
    swc  = x.s(W(i).cells,:);
    lamc = bsxfun(@rdivide,fluid.relperm(swc), mu); 
    fc   = lamc(:,1)./sum(lamc,2);
    
    sww  = W(i).compi;
    lamw = bsxfun(@rdivide,fluid.relperm(sww), mu);
    fw   = lamw(:,1)./sum(lamw,2);
    
    wellSol(i).name = W(i).name;
    wellSol(i).bhp  = x.wellSol(i).pressure;
    wellSol(i).wcut = iout.*fc + iin.*fw;
    wellSol(i).Sw   = iout.*swc(:,1) + iin.*sww(:,1);
    wellSol(i).qW   = out.*fc + in.*fw;
    wellSol(i).qWs  = sum(wellSol(i).qW);
    wellSol(i).qO   = out.*(1-fc) + in.*(1-fw);
    wellSol(i).qOs  = sum(wellSol(i).qO);
end<jupyter_output>Writing /content/mrst-core/modules/getWellSol.m<jupyter_text>Access SPE 10 model data<jupyter_code>!mkdir /content/mrst-core/data
!mv /content/reservoir_datasets/MRST/SPE10 /content/mrst-core/data<jupyter_output><empty_output><jupyter_text>Execute main program The result is OIP after 5 years of production<jupyter_code>!octave -W SPE10_waterflooding.m
out = !octave -W SPE10_waterflooding.m
np.float64(out)<jupyter_output><empty_output><jupyter_text>Well placement optimization with Python<jupyter_code>import numpy as np
!pip install bayesian-optimization
from bayes_opt import BayesianOptimization
%%writefile calculate_oip_spe10.m
%% Script that reads the well location from input text file, and then calculates
%% OIP after several years of production from "SPE10_waterflooding.m" program

%% Read well location from input text file
fileID = fopen('input.txt','r');
formatSpec = '%f';
P = fscanf(fileID,formatSpec);

%% Output is 1D array. Reshape to x and y coords
wloc = reshape(P,2,[]);

%% The whole SPE10_waterflooding.m program
addpaths

%% Set up model and parameters
T           = 5*year;
nlayer      = 1; % number of layers of SPE 10 model to be simulated
startlayer  = 4; % the SPE 10 model starts from layer i-th defined
dims        = [60 220 nlayer];
domain      = dims.*[20 10 2]*ft;
G           = computeGeometry(cartGrid(dims,domain));

rock        = getSPE10rock((1:dims(1)), (1:dims(2)), ... 
              startlayer:startlayer+(nlayer-1));
            
rock.poro = max(rock.poro,.0005);

hT = computeTrans(G, rock); % compute transmissibility 
%x0 = initState(G,W,100*barsa, [0 1]); % initialize

%% Set fluid model
fluid = initSimpleFluid('mu' , [   1,    1] .* centi*poise     , ...
                        'rho', [1000,  850] .* kilogram/meter^3, ...
                        'n'  , [   2,    2]);

%% Set parameters describing the well. wloc has been updated above.
wtype    = {'bhp', 'bhp', 'bhp', 'bhp', 'bhp'};
wtarget  = [200,   200,   200,   200,   500] .* barsa();
wrad     = [0.125, 0.125, 0.125, 0.125, 0.125] .* meter;
wname    = {'P1', 'P2', 'P3', 'P4', 'I'};

W = [];
for w = 1 : numel(wtype)
  if w < 5
    % production well is oil
    phase = [0 1];
  else
    % production well is water
    phase = [1 0];
  endif
  W = verticalWell(W, G, rock, wloc(1,w), wloc(2,w), 1, ...
      'Type', wtype{w}, 'Val', wtarget(w), ...
      'Radius', wrad(w), 'Name', wname{w}, ...
      'InnerProduct', 'ip_tpf', 'Comp_i', phase);
end

x0 = initState(G,W,100*barsa, [0 1]); % initialize

%% Simulation
M = 50;
x = x0;
[dt,dT] = deal(zeros(M,1), T/M);
wellSol = cell(M,3);
oip     = zeros(M,3);
n = 1;

for i=1:M
    x  = incompTPFA(x, G, hT, fluid, 'wells', W);
    x  = implicitTransport(x, G, dT, rock, fluid, 'wells', W);    
        
    dt(i) = dT;
    oip(i,n) = sum(x.s(:,2).*poreVolume(G,rock));
    wellSol{i,n} = getWellSol(W,x, fluid);
end

%% Print the Oil in Place
display(oip(50));
def objective(x, y):
  """
  The objective function returns the OIP (calculated from the MATLAB program
  calculate_oip_spe10.m) that will be minimized

  An optimizer will optimize the injector locations ONE BY ONE. This function
  should be put in a loop.

  Input:

  x, y coordinates of one of the injector wells
  """
  # change location of one of the injectors
  wloc[i] = np.array([x, y])

  # save the locations into input text file
  np.savetxt('input.txt', wloc)

  # execute calculate_oip_spe10.m MATLAB program
  oip = !octave -W calculate_oip_spe10.m

  return np.float64(oip)
# Initial injector well locations
wloc = np.array([[1,1],[60,1],[1,220],[60,220],[30,111]])
x = wloc[:,0]
y = wloc[:,1]
objective(11.67,100.76)
for i in range(2):
  # Define search bound for each injector
  if i == 0:
    # First well (SW)
    pbounds = {"x": (0, 30), "y": (0, 111)}
  if i == 1:
    # Second well (SE)
    pbounds = {"x": (30, 60), "y": (0, 111)}
  if i == 2:
    # Third well (NW)
    pbounds = {"x": (0, 30), "y": (111, 220)}
  if i == 3:
    # Fourth well (NE)
    pbounds = {"x": (30, 60), "y": (111, 220)}
  
  # Optimization
  optimizer = BayesianOptimization(
      f=objective,
      pbounds=pbounds,
      verbose=2,  # verbose = 1 prints only when a maximum is observed, verbose = 0 is silent
      random_state=42,
  )

  print('Optimization at injector well:', i+1)
  optimizer.maximize(init_points=5, n_iter=20)

  x_reloc = optimizer.max['params']['x']
  y_reloc = optimizer.max['params']['y']<jupyter_output>Optimization at injector well: 1
|   iter    |  target   |     x     |     y     |
-------------------------------------------------
{'fraction_non_alphanumeric': 0.1610512129380054, 'fraction_numerical': 0.04016172506738545, 'mean_word_length': 2.8034341363403383, 'pattern_counts': {'":': 8, '<': 36, '<?xml version=': 0, '>': 26, 'https://': 4, 'lorem ipsum': 0, 'www.': 1, 'xml': 0}, 'pii_count': 1, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 8, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '7', 'dir': 'jupyter-structured-clean-dedup', 'id': '564067', 'license': 'MIT', 'n_tokens_mistral': 5754, 'n_tokens_neox': 5244, 'n_words': 1433, 'path': 'MRST_SPE10_optimization_python.ipynb', 'repo_name': 'yohanesnuwara/MRST-simulations', 'size': '14840'}