The initialization component of the program is executed once per simulation. The routines in the initialization component are executed sequentially, as shown in the flow chart, from the program start to the start of the first time step.

During initialization, the input file is read twice.The first read of the input file is used to allocate and initialize memory for all the global variables (i.e., those variables passed among routines via modules).The second read of the input file is used to define the problem.The input file readers are not foolproof but do report error messages when input errors are noted.Once the input file has been successfully read twice, the initial states are checked for thermodynamic consistency, phase conditions are established, and the primary variables are selected for every grid cell.The next stage of initialization involves computing all the secondary variables at grid cells and boundary surfaces.If geomechanics calculations are included in the simulation, then the initial stress state of the domain is established.If coupled wells are included in the simulation, then the well trajectories are determined and the wells are equilibrated with the formation.Once the well trajectories are defined, the coupled conservation equations and coupled well equations are sequenced such that the band width of the Jacobian matrix is minimized.If geochemistry calculations are included in the simulation, then the initial chemical state of the system is established.With the initial primary and secondary variables computed, the next stage of initialization is to compute the initial surface fluxes at both internal surfaces and boundary surfaces (e.g., volumetric phase flux, component advective/diffusive flux, heat flux).The final initialization stage is to record user requested output to the various output files (i.e., output, plot.xxx, and surface).If a zero-time-step simulation is requested, the simulator stops at the end of the initialization, after recording output.This option is very useful for computing property values or examining initial states.

The iteration component of the program contains a pair of nested loops. The outer loop performs time stepping and the inner loop performs Newton-Raphson iteration.

• Termination of the Newton-Raphson loop occurs upon successful convergence or with an iteration limit violation.
• Termination of the time-stepping loop occurs due to the simulation limit or a time step reduction limit violation.
• It should be noted that "Solute Transport" is shown as a single routine on the STOMP flow diagram; however, it comprises several transport routines within a solute loop.

The time stepping component is where the simulator loops over time, with each time-step loop solving for the conditions at a new point in time.Geomechanics, reactive transport (geochemistry), and solute transport are solved sequentially with the coupled multifluid flow and heat transport.Each new time step begins with the assignment of a time step, which is determined from the previous time step parameters (e.g., time step and convergence) and user input (e.g., maximum time step, minimum time step, time-step acceleration factor, output times, boundary condition times, source times, coupled well times).The next step is to store the old time step values of the primary and secondary variables.Old-time-step values of fluxes are not stored.Next, boundary conditions and sources at the new time step are assigned.Fluxes across the interior and boundary surfaces are then computed.The Jacobian matrix, solution vector, and problem vector are all set to zero.The Jacobian matrix and problem vector are then loaded for all active nodes and interior surfaces (i.e., non-boundary surfaces).The Jacobian matrix and problem vector is then corrected for boundary surface fluxes.At this point the Jacobian matrix and problem vector are complete and the linear system is either solved directly with a banded solver or iteratively.The resulting solution vector from the linear system solve contains corrections (updates) to the primary variables.The primary variables are then updated and residuals (errors) in the conservation equations are determined.The updated primary variables are then used to establish new phase conditions and a new primary variable set for each grid cell.Increments to the primary variables for the numerical derivatives are determined, and the secondary variables are computed using native and incremented primary variables.Convergence is then established using the conservation equation residuals; where the highest relative residual across the computational domain determines convergence.

The closure routines are executed at the completion of a simulation (regardless of the termination cause).

If convergence has occurred for all conservation equations and for all grid cells, then the simulation proceeds to the geomechanics, solute transport and reactive transport routines. Prior to calling these routines, integrals for sources are computed and surface fluxes are updated using the converged values of the primary and secondary variables. Before starting a new time step, output is recorded and the simulation limits are checked. If the simulation time or number of time steps has reached a user defined limit, then the simulation moves to closure at which point output is recorded, all simulation files are closed, and the program stops. If convergence has not occurred, then the iteration limit is checked. If the number of Newton-Raphson iterations exceeds the user-defined limit, then the time step is reduced and the time step is repeated at a reduced time step amount. If the number of Newton-Raphson iterations is less than the user-defined limit, then another iteration loop occurs, starting with a calculation of the boundary conditions. If the iteration limit is exceeded, then the new reduced time step is checked against the minimum time step or the number of allowed consecutive time-step reductions is checked. If either of these checks fail, then the simulation transfers to closure, resulting in recording of the output and cessation of the simulation.