MRST - Transforming research on reservoir simulation

Release Notes for MRST 2020a

Highlights of MRST 2020a

The new release comes with a number of bug fixes, improvements and increased documentation coverage. We especially highlight:

  • Several new modules for domain decomposition and standardized setup of examples.
  • A much faster direct-assembly AD backend with a corresponding linear solver for fast simulations.
  • The full historical Norne model can now be simulated with MRST after update to many features for initialization and wells.
  • Utilities that make it easier to set up non-standard schemes such as explicit, WENO, AIM and MPFA for the AD-based solvers.
  • A new MPFA-O implementation that handles general polyhedra.
  • New ensemble functionality and GUIs for flow diagnostics analysis of ensembles.
  • Functionality to add wells by specifying a well trajectory.

Changes and updates to core functionality

The ‘processGRDECL’ function now knows how to preserve the input data’s original corner points. This feature is enabled by the ‘PreserveCpNodes’ option. If available, such corner points will enable the ‘computeGeometry’ function to calculate cell- and face centroids that coincide with those of ECLIPSE which in turn enable calculating ECLIPSE-compatible transmissibility values. We also capture transmissibility multiplier arrays such as MULTX and MULTZ- in the ‘grdecl2Rock’ function and incorporate these into MRST’s transmissibility calculation.

To that end the ‘getFaceTransmissibility’ function has been moved from the ‘ad-core’ module and into MRST’s core feature set. This function is aware of corner-point geometry and knows how to incorporate ECLIPSE-style transmissibility multipliers. Consequently, this function enables calculating ECLIPSE-compatible transmissibility values. For backwards compatibility with earlier versions of MRST this function still recognises a ‘deck’-style third input argument and will extract multiplier information if available and not already stored in the ‘rock’. This syntax will issue a diagnostic message however and will be removed in a future version of MRST. We recommend that ECLIPSE-style multipliers be passed in the ‘rock’ parameter and that any ‘deck’ arguments be removed at the call site.

The ADI class now supports sine and cosine as well as inverse sine (asin) and inverse cosine (acos) operations.

We have added a new ‘tensor assembly’ package that contains tools to handle multiple indices in a generic framework based on a tensorial approach. These tools are used in the new implementation of the MPFA-O method and in the VAG module.

The new function ‘computeVerticalGridIntersection’ computes all polygons resulting from cutting a grid with a piecewise linear vertical surface. This is a foundational tool for expanded capabilities in the wellpaths module.


New modules


Module example-suite

The example-suite implements a way to structure and maintain examples. The suite consists of a set of functions defining a number of different examples. These can either be called as stand-alone functions that return an initial state, model and simulation schedule, or in conjunction with the new MRSTExample class, which implements functionality for setting up a simulation problem with a reasonable solver configuration, and functionality for generating nice plots. The class also provides access to the entire example setup, facilitating easy setup of different configurations of the same example. In the future, the class can be used for e.g. setting up ensembles, and for regression testing.

Module domain-decomposition

This module implements nonlinear variable domain decomposition. A full model is partitioned into subdomains, each defining a subdomain model with the class SubdomainModel. The module has utility functions for subtracting the corresponding substate and subschedule, so that the entire simulation can be carried out in the subdomain only. Nonlinear domain decomposition is implemented in the DomainDecompositionModel, which at each timestep solves each subdomain independently, treating all other subdomains as fixed (additive) or serially, with updated values in subdomains that are already solved (multiplicative). Convergence of the full problem is ensured by checking the full residual after the subdomain solves. Instead of defining boundary conditions, SubdomainModel evaluates the model equations in its subdomain cells, in addition to one layer of external cells along boundaries against adjacent subdomains, and the linearized problem is altered so that values in external cells remain unchanged. This way, we avoid imposing potentially complicated and error-prone boundary conditions between subdomains.

The module also supports dynamic subdomain partitions based on the current reservoir state through the Partition class. For instance, optimal reordering of transport equations based on the intercell flux graph can be achieved with multiplicative nonlinear domain decomposition using the TopologicalFluxPartition.


Changes to existing modules



  • The LinearizedProblem’s ‘numel’ method has been renamed to ‘numeq’ to avoid conflicting with Octave which treats ‘numel’ as the way to retrieve the number of elements in an array of ‘LinearizedProblem’ objects. This was not widely used in MRST, but please take care if you rely on the functionality.
  • Several new options for diagonal automatic differentiation backends for fast assembly. This includes a large number of performance improvements, an option to use rowMajor storage for better memory locality and an option to defer assembly to make it easier to assemble into new formats. As an example of this functionality, we include the possibility of outputting block CSR matrices instead of the usual scalar CSC. In addition, two AMGCL solvers that take in the CSC format directly are included. With the caveat that wells should be perforated in a single cell, this can be significantly faster for larger cases.
  • blockAssemblyBigModelExample: New example demonstrating the new block solvers and different backends
  • MultipliedPoreVolume state function has been split into a base class PoreVolume for static volumes, and BlackOilPoreVolume for the version that uses pvMultR.
  • StateFunctionGrouping can be instantiated on its own.
  • FlowPropertyFunctions has been split off, forming a new class PVTPropertyFunctions that hold PVT-related functions.
  • ReservoirModel: Now instantiates FacilityModel by default.
  • State functions: Many improvements, documentation. Some functions have been given more appropriate names. Please see full git log at
  • NonLinearSolver: oscillationThreshold property can be adjusted from the default value of 1 to make dampening/linesearch trigger more easily.
  • AIM/Explicit solvers have been given a tune-up.
  • WrapperModel has been updated with many missing parent functions.
  • Pressure reduction factors are now scaled with 1/pv to improve scaling of resulting pressure equation
  • ensemblePackedProblemsExample: New example demonstrating how to simulate an ensemble of cases in parallel
  • Improved initialization logic for corner cases (zero pc curves, immobile phases, support for SWATINIT pc scaling)
  • Rewritten support for threshold pressures
  • printStateFunctionGroupingTikz: New function for making a LaTeX/pgfplots version of a state function diagram
  • Bug fixes to WrapperModel
  • Rel. perm. Scaling is now implemented by a virtual class that also takes care of capillary scaling
  • Much better logic for ordering of wells when schedule is complex (should match better with Eclipse, also includes several variations useful for reordering long wells)
  • MPFAvsTPFAwithADExample: Example demonstrating how to run MPFA with generic solvers.
  • convertDeckScheduleToMRST: Option to silence warnings from well processing. A few of the more noisy default warnings have been disabled.
  • ComponentTotalMass now has option to set minimum derivatives on the major diagonal bands, which can help for cases where an iterative solver is being used
  • getWellOutput: Support for getting multi-valued fields from wellsols
  • NumericalPressureReductionFactors added that automatically computes correct pressure weights for any multicomponent model.
  • initEclipsePackedProblemAD: Version of initEclipseProblemAD that outputs a packed problem with result handlers and default nonlinear solver setup.
  • copyPackedProblem: Function to copy a packed problem while updating the storage locations
  • Utilities buildMexOperators, buildLinearSolvers make it easy to quickly build parts of MRST’s MEX acceleration
  • AMGCL wrappers now have nicer output to terminal, with explanation of the options that are set.
  • Added FactorTimeStepSelector: Simple step selector that increases and decreases timesteps by multiplying with a constant factor
  • Ad-core: Bug-fixes in the computeSensitivitiesAdjoint
  • Many smaller optimizations in framework to avoid computing redundant quantities
  • Rename DiagonalSubset -> FixedWidthJacobian
  • Wells can be solved separately from reservoir equations before each iteration by passing a FacilityNonLinearSolver to model constructor
  • plotPackedProblem: New function to plot output from a packed problem
  • getReportOutput: New function for getting aggregated data from a simulation report
  • simulatePackedProblemBackground: Store path setup when simulating in the background
  • Utilities setTimeDiscretization, setWENODiscretization, setMPFADiscretization for quickly changing to alternative schemes


  • Bug fixes to component density for multiple PVT regions
  • fieldModelNorneExample that demonstrates how to run the Norne model in MRST.


  • GenericSurfactantPolymerModel: A generic model to handle the surfactant and polymer EOR simulation with black-oil model based on the new AD-OO framework with state functions. With this model, simulations with different combinations of surfactant, polymer, oil, water and gas can be achieved.


  • Improved support for MEX interpolation in input tables.
  • Experimental options to speed up table lookup even further if the tables are given on uniform grids.
  • More keywords now give out function handles on the new format for multiple regions


  • compositionalHybridUpwind: New example that demonstrates hybrid upwind for a multicomponent problem with gravity.


  • Improved mass-conservation checks when water is present.



  • Now contains functionality from the opm_gridprocessing module.


  • New diagnosticsViewer GUI which calculates and displays flow diagnostics calculations for an ensemble of models sharing the same grid and well structure. This allows for interactive viewing in 3D of flow diagnostics results prior to any multiphase simulations being run.
  • New ensembleGUI to show flow diagnostics properties for an ensemble of models, calculated using residence time distributions.
  • New ensemble class to easily create, store and retrieve an ensemble of models with accompanying flow diagnostics results and full simulation results.


  • Updated to comply with latest changes in autodiff.
  • Use simulatePackedProblem in geothermalExampleHTATES.


  • The fluid object construction functions ‘init-Fluid-’ now form structures that are compatible with MRST’s AD-based solver routines. In particular the structures now support queries of the form.
    • krO = fluid.krO(so)
    • krW = fluid.krW(sw)
    • muO = fluid.muO(po)
  • For backwards compatibility the structures still support the original ‘properties’ and ‘relperm’ functions, but we expect to deprecate or remove these functions in a future release of MRST.


  • buildLinearSolvers: Dedicated build routine for all linear solvers
  • Support for building AMGCL in GNU Octave.
  • New AMGCL interface that uses block CSR as input. See ad-core changelog.
  • Support for initial guess in AMGCL solvers


  • A new alternative implementation of MPFA that uses the new tensor assembly tools from core. The legacy implementation (called by default) is limited to a certain type of grid cells. The new version can handle general polyhedrals. There is still a limitation for hanging nodes.
  • Large systems in the new solver can use assembly by blocks to deal with memory limitations in Matlab.
  • Examples linearPressureTestMPFA.m and mpfatest.m compare the different implementations.


This module, which was first released with MRST 2019b, was developed by W. Zhang and M. Al Kobaisi at Khalifa University of Science and Technology. It adds support for nonlinear finite volume discretizations (NFVM). It follows the development of paper W. Zhang, M. Al Kobaisi, Cell-centered nonlinear finite volume methods with improved robustness, SPE Journal (, 2019. The current release includes:

  • Support the new object oriented discretization framework. The previous interface is discarded, and with that the old examples except runHAP.
  • runHAP.m is rewritten to support the new discretization framework.
  • An example comparing TPFA, MPFA and NFVM is added.
  • Miscellaneous performance improvements, including a mex routine for determining if a point is inside the convex hull of a set of points (3D only).


  • Now released as a supplementary module which must be downloaded from the MRST website (
  • Updated with several new patches to improve MRST compatibility with Octave.
  • Improved support for mex-compilation.
  • Improved support for class-based architecture of the MRST automatic-differentiation based solvers.


  • This module has been subsumed into the “deckformat” module and no longer exists separately. If you directly use either of the ‘mprocessGRDECL’ or ‘processgrid’ functions, then please ensure that you also activate the ‘deckformat’ module.


  • New class OptimizationProblem that can handle optimization over ensembles.
  • Several utility-functions that enable a version of the generalized reduced gradient (GRG) method. This means that for objective functions including simulations, both well targets and well limits can be optimized simultaneously (hence providing a control policy rather than a single trajectory).
  • Two examples are included to illustrate usage: one small (optimizeSimpleEnsemble1D) and one moderate (optimizeEggEnsemble).
  • The OptimizationProblem-class supports running evaluations and/or the entire optimization in multiple background processes for speed-up.


  • Remove scaling of solvent equation, since the model uses CNV convergence.
  • Improved spe10WAGexample with iterative linear solver.


  • Now uses the latest version of the tensor assembly functionalities.


  • The new function ‘computeTraversedCells’ traces a piecewise linear trajectory through a grid and records the intersected cells and the coordinates of its face intersections. This enables calculating more accurate well indices (connection transmissibility factors) using the ‘lineSegments’ option of the core function ‘addWell’.
  • New function called ‘addWellFromTrajectory’ that mimics the behaviour of ‘addWell’, but which is defined in terms of a well trajectory rather than a sequence of cell indices.

Published December 16, 2019