Understand Compiling#

Author: Zakariya Abugrin | Date: May 2025

Introduction#

In this tutorial, we see how a model can be compiled to calculate different types of solutions such as: analytical, numerical, and neurical which is one of the main research topics in reservoirflow.

Hint

Compiling solutions is the most interesting idea introduced in reservoirflow which allows to solve the same model using different solutions so we can compare them with each other and/or combine them together.

../../_images/tutorial_understand_compiling_solutions_diagram_dark.svg

../../_images/tutorial_understand_compiling_solutions_diagram_light.svg

\

The following points are important to understand the compiling concept in reservoirflow:

Every model from models module has a compile() method which can be used to select a solution.
Every model come with an empty dictionary of solutions and compilers accessed by model.solutions and model.compilers.
Every model must be compiled before the solution can be computed.
After compiling a model, a new attribute called model.solution is added.
Every model can be solved using different solutions where model.set_solution() can be used to switch to previous solutions.
Solutions can be computed using model.solution.solve() for a single time step and model.solution.run() for multiple time steps.
Solution functions are mapped directly to the corresponding model and can be accessed directly using model.solve() and model.run().
Specifically, neurical solutions (i.e. based on neural networks) have model.solution.fit() to train based on physics-loss and model.solution.predict().
Specifically, neurical solution functions are mapped directly to the corresponding model and can be accessed directly using model.train() and model.predict().

Now, let’s see how we can apply this concept.

Import `reservoirflow`#

We start with importing reservoirflow as rf. The abbreviation rf refers to reservoirflow where all modules under this library can be accessed. rf is also used throughout the API documentation. We recommend our users to stick with this convention.

import reservoirflow as rf

print(rf.__version__)

0.1.0b3

# Check what is available in solutions:
dir(rf.solutions)

['Compiler',
 'Solution',
 '__all__',
 '__builtins__',
 '__cached__',
 '__doc__',
 '__file__',
 '__loader__',
 '__name__',
 '__package__',
 '__path__',
 '__spec__',
 'compiler',
 'solution']

Warning

Modules under solutions are not loaded by default in reservoirflow. As a result, rf.solutions.numerical.FDM will fail. As can be seen below, only classes Compiler and Solution are loaded by default.

# rf.solutions.numerical.fdm.FDM # this will not work
# rf.solutions.numerical.FDM # this will not work

# But, we can load these modules or classes explicitly:
# from reservoirflow.solutions import numerical
# from reservoirflow.solutions.numerical import FDM

Note

By default, modules numerical, analytical, and neurical are not available under solutions.

Tip

Specific solutions modules such as rf.solutions.numerical will only be available once a solution for that module (e.g. rf.solutions.numerical.FDM) was used to compile a model. Of course, this does not prevent loading these modules explicitly (e.g. from reservoirflow.solutions import numerical, or from reservoirflow.solutions.numerical import FDM).

Build a model#

A reservoir simulation model requires two objects: Grid and Fluid. The function below create_model() starts by creating these objects which are used to initiate a Model object using BlackOil class.

def create_model():
    # Grid:
    grid = rf.grids.RegularCartesian(
        nx=4,
        ny=1,
        nz=1,
        dx=300,
        dy=350,
        dz=40,
        phi=0.27,
        kx=270,
        dtype="double",
    )

    # Fluid:
    fluid = rf.fluids.SinglePhase(mu=0.5, B=1, dtype="double")

    # Model:
    model = rf.models.BlackOil(
        grid,
        fluid,
        pi=4000,
        dt=1,
        start_date="10.10.2018",
        verbose=False,
        dtype="double",
    )

    # Production well:
    model.set_well(cell_id=4, q=-600, s=1.5, r=3.5)

    # Boundaries:
    model.set_boundaries({0: ("pressure", 4000), 5: ("rate", 0)})

    return model


model = create_model()

Simulation Run#

To perform the simulation run, method model.run() can be used. The code below performs a simulation run for nsteps=10 (i.e. number of steps) and using isolver=cgs:

model.solutions

{}

model.solve()  # same as: model.run()

The model is not compiled. Use model.compile() to add solve() and run() methods.

help(model.solve)  # same as: help(model.run)

Help on method solve in module reservoirflow.models.model:

solve(**kwargs) method of reservoirflow.models.black_oil.BlackOil instance
    Solve a single timestep.
    
    .. attention::
        This method is not available until the model is compiled
        using ``model.compile()``.
    
    Once the model is compiled, the documentation of the assigned
    solution can be accessed using one of the following methods:
    
    >>> help(model.solve) # or help(model.solution.solve)
    >>> print(model.solve.__doc__) # or print(model.solution.solve.__doc__)

model.compile(stype="numerical", method="FDM", sparse=False)
# model.compile(stype="analytical", method="1D1P", sparse=False)

[info] FDM was assigned as model.solution.

model.compiler

rf.solutions.Compiler(model='BlackOil', stype='numerical', method='FDM', sparse=False)

print(model.run.__doc__)

Perform a simulation run for nsteps.

        Parameters
        ----------
        nsteps : int, optional
            _description_
        threading : bool, optional
            _description_
        check_MB : bool, optional
            _description_
        isolver : str, optional
            iterative solver for sparse matrices. Available solvers are
            ["bicg", "bicgstab", "cg", "cgs", "gmres", "lgmres",
            "minres", "qmr", "gcrotmk", "tfqmr"].
            If None, direct solver is used. Only relevant when argument
            sparse=True. Direct solver is recommended for more accurate
            calculations. To improve performance, "cgs" is recommended
            to increase performance while option "minres" is not recommended due to
            high MB error. For more information check [1][2].

        References
        ----------
        - SciPy: `Solving Linear Problems <https://docs.scipy.org/doc/scipy/reference/sparse.linalg.html#solving-linear-problems>`_.
        - SciPy: `Iterative Solvers <https://scipy-lectures.org/advanced/scipy_sparse/solvers.html#iterative-solvers>`_.
        

help(model.run)

Help on method run in module reservoirflow.solutions.numerical.fdm:

run(nsteps=10, threading=True, vectorize=True, check_MB=True, print_arrays=False, isolver=None) method of reservoirflow.solutions.numerical.fdm.FDM instance
    Perform a simulation run for nsteps.
    
    Parameters
    ----------
    nsteps : int, optional
        _description_
    threading : bool, optional
        _description_
    check_MB : bool, optional
        _description_
    isolver : str, optional
        iterative solver for sparse matrices. Available solvers are
        ["bicg", "bicgstab", "cg", "cgs", "gmres", "lgmres",
        "minres", "qmr", "gcrotmk", "tfqmr"].
        If None, direct solver is used. Only relevant when argument
        sparse=True. Direct solver is recommended for more accurate
        calculations. To improve performance, "cgs" is recommended
        to increase performance while option "minres" is not recommended due to
        high MB error. For more information check [1][2].
    
    References
    ----------
    - SciPy: `Solving Linear Problems <https://docs.scipy.org/doc/scipy/reference/sparse.linalg.html#solving-linear-problems>`_.
    - SciPy: `Iterative Solvers <https://scipy-lectures.org/advanced/scipy_sparse/solvers.html#iterative-solvers>`_.

# print(rf.solutions.numerical.fdm.FDM.run.__doc__)
# print(rf.solutions.numerical.FDM.run.__doc__)
# print(model.solution.run.__doc__)
print(model.run.__doc__)

Perform a simulation run for nsteps.

        Parameters
        ----------
        nsteps : int, optional
            _description_
        threading : bool, optional
            _description_
        check_MB : bool, optional
            _description_
        isolver : str, optional
            iterative solver for sparse matrices. Available solvers are
            ["bicg", "bicgstab", "cg", "cgs", "gmres", "lgmres",
            "minres", "qmr", "gcrotmk", "tfqmr"].
            If None, direct solver is used. Only relevant when argument
            sparse=True. Direct solver is recommended for more accurate
            calculations. To improve performance, "cgs" is recommended
            to increase performance while option "minres" is not recommended due to
            high MB error. For more information check [1][2].

        References
        ----------
        - SciPy: `Solving Linear Problems <https://docs.scipy.org/doc/scipy/reference/sparse.linalg.html#solving-linear-problems>`_.
        - SciPy: `Iterative Solvers <https://scipy-lectures.org/advanced/scipy_sparse/solvers.html#iterative-solvers>`_.
        

dir(rf.solutions)

['Compiler',
 'Solution',
 '__all__',
 '__builtins__',
 '__cached__',
 '__doc__',
 '__file__',
 '__loader__',
 '__name__',
 '__package__',
 '__path__',
 '__spec__',
 'compiler',
 'numerical',
 'solution']

dir(rf.solutions.numerical)

['FDM',
 'FEM',
 'FVM',
 '__all__',
 '__builtins__',
 '__cached__',
 '__doc__',
 '__file__',
 '__loader__',
 '__name__',
 '__package__',
 '__path__',
 '__spec__',
 'fdm',
 'fem',
 'fvm',
 'solvers']

print(
    " stype:",
    model.compiler.stype,
    "\n",
    "method:",
    model.compiler.method,
    # "\n",
    # "mode:",
    # model.compiler.mode,
    # "\n",
    # "solver:",
    # model.compiler.solver,
    "\n",
    "solution:",
    # model.compiler.solution,
    model.compiler.model.solution,
    model.compiler.model.name,
    "\n",
    model.solutions,
)

 stype: numerical 
 method: FDM 
 solution: <reservoirflow.solutions.numerical.fdm.FDM object at 0x7f5b37118150> BlackOil 
 {'BlackOil-numerical-FDM-False': <reservoirflow.solutions.numerical.fdm.FDM object at 0x7f5b37118150>}

model.run(
    nsteps=10,
    vectorize=True,
    threading=True,
    isolver="cgs",
)

[info] Simulation run started: 10 timesteps.
[info] Simulation run of 10 steps finished in 0.04 seconds.
[info] Material Balance Error: 1.693933882052079e-11.

model.get_df()

	Time	Q0	Q4	P0	P1	P2	P3	P4	Qw4	Pwf4
Step
0	0	0.0	0.0	4000.0	4000.000000	4000.000000	4000.000000	4000.000000	0.0	4000.000000
1	1	600.0	-600.0	4000.0	3989.436768	3968.310305	3947.183842	3926.057379	-600.0	3898.992647
2	2	600.0	-600.0	4000.0	3989.436768	3968.310305	3947.183842	3926.057379	-600.0	3898.992647
3	3	600.0	-600.0	4000.0	3989.436768	3968.310305	3947.183842	3926.057379	-600.0	3898.992647
4	4	600.0	-600.0	4000.0	3989.436768	3968.310305	3947.183842	3926.057379	-600.0	3898.992647
5	5	600.0	-600.0	4000.0	3989.436768	3968.310305	3947.183842	3926.057379	-600.0	3898.992647
6	6	600.0	-600.0	4000.0	3989.436768	3968.310305	3947.183842	3926.057379	-600.0	3898.992647
7	7	600.0	-600.0	4000.0	3989.436768	3968.310305	3947.183842	3926.057379	-600.0	3898.992647
8	8	600.0	-600.0	4000.0	3989.436768	3968.310305	3947.183842	3926.057379	-600.0	3898.992647
9	9	600.0	-600.0	4000.0	3989.436768	3968.310305	3947.183842	3926.057379	-600.0	3898.992647
10	10	600.0	-600.0	4000.0	3989.436768	3968.310305	3947.183842	3926.057379	-600.0	3898.992647

model.solution.get_cells_eq()

{np.int64(1): (defaultdict(int, {p1: 85.2012000000000, p2: -28.4004000000000}),
  227203.200000000),
 np.int64(2): (defaultdict(int,
              {p1: 28.4004000000000,
               p3: 28.4004000000000,
               p2: -56.8008000000000}),
  0.0),
 np.int64(3): (defaultdict(int,
              {p2: 28.4004000000000,
               p4: 28.4004000000000,
               p3: -56.8008000000000}),
  0.0),
 np.int64(4): (defaultdict(int, {p3: 28.4004000000000, p4: -28.4004000000000}),
  600.000000000000)}

model.solution.get_matrices_symb(True)

(array([[ 85.2012, -28.4004,   0.    ,   0.    ],
        [ 28.4004, -56.8008,  28.4004,   0.    ],
        [  0.    ,  28.4004, -56.8008,  28.4004],
        [  0.    ,   0.    ,  28.4004, -28.4004]]),
 array([[227203.2],
        [     0. ],
        [     0. ],
        [   600. ]]))

model.solution.get_matrices_vect(True)

(array([[-85.2012,  28.4004,   0.    ,   0.    ],
        [ 28.4004, -56.8008,  28.4004,   0.    ],
        [  0.    ,  28.4004, -56.8008,  28.4004],
        [  0.    ,   0.    ,  28.4004, -28.4004]]),
 array([[-227203.2],
        [      0. ],
        [      0. ],
        [    600. ]]))

model.solve(print_arrays=True)

step: 10
[[ 85.2012 -28.4004   0.       0.    ]
 [ 28.4004 -56.8008  28.4004   0.    ]
 [  0.      28.4004 -56.8008  28.4004]
 [  0.       0.      28.4004 -28.4004]
 [-85.2012  28.4004   0.       0.    ]
 [ 28.4004 -56.8008  28.4004   0.    ]
 [  0.      28.4004 -56.8008  28.4004]
 [  0.       0.      28.4004 -28.4004]
 [  0.       0.       0.       0.    ]
 [  0.       0.       0.       0.    ]
 [  0.       0.       0.       0.    ]
 [  0.       0.       0.       0.    ]]
[[ 227203.2 -227203.2       0. ]
 [      0.        0.        0. ]
 [      0.        0.        0. ]
 [    600.      600.        0. ]]

Comments 💬#

Feel free to make a comment, ask a question, or share your opinion about this specific content. Please keep in mind the Commenting Guidelines ⚖.