Doc and README pass

README.md

@@ -3,43 +3,53 @@
 [![Build Status](https://github.com/sintefmath/JutulDarcy.jl/actions/workflows/CI.yml/badge.svg?branch=main)](https://github.com/sintefmath/JutulDarcy.jl/actions/workflows/CI.yml?query=branch%3Amain)
 
 # JutulDarcy.jl
-Darcy-scale and subsurface flow using [Jutul.jl](https://github.com/sintefmath/Jutul.jl) developed by the [Computational Geosciences group](https://www.sintef.no/en/digital/departments-new/applied-mathematics/computational-geoscience/) at [SINTEF Digital](https://www.sintef.no/en/digital/).
+Darcy-scale and subsurface flow (CO2 sequestration, gas/H2 storage, oil/gas fields) using [Jutul.jl](https://github.com/sintefmath/Jutul.jl) developed by the [Computational Geosciences group](https://www.sintef.no/en/digital/departments-new/applied-mathematics/computational-geoscience/) at [SINTEF Digital](https://www.sintef.no/en/digital/).
 
 ## Key features
-- Written in pure Julia, with automatic differentiation
+- Written in pure Julia, with automatic differentiation and dynamic sparsity detection
+- Support for sensitivities with respect to any model parameters using the adjoint method
 - High performance assembly and linear solvers, with support for two-stage CPR BILU(0)-CPR Krylov solvers
 - Equation-of-state compositional, immiscible and black oil flow is supported and validated against existing simulators
 - Unstructured grids and complex cases input from [the Matlab Reservoir Simulation Toolbox (MRST)](https://www.mrst.no) using the `jutul` module.
 - Support for general multisegment wells with rigorous mass balance, complex well limits and time-dependent controls
-- 3D visualization of grids and wells
-- Interactive plotting of well curves (somewhat experimental)
+- 3D visualization of grids and wells in [JutulViz.jl](https://github.com/sintefmath/JutulViz.jl)
+- Interactive plotting of well curves
 
 The compositional simulator has been matched against commercial offerings, AD-GPRS and MRST. The blackoil simulator has been validated on the standard SPE benchmarks (SPE1, SPE9, ...).
 
+## Example run times on benchmarks
+| Name      | Cells | Report steps | Preconditioner   | Time [s] |
+|-----------|-------|--------------|------------------|----------|
+| SPE1CASE2 | 300   | 120          | block-ILU(0)     | 0.85     |
+| SPE9      | 9000  | 35           | block-ILU(0)     | 9.30     |
+| Egg       | 18553 | 123          | CPR-block-ILU(0) | 22.5     |
+Simulated with `julia -O2`, no threads.
+
 ## A few of the packages used by Jutul and JutulDarcy
 Jutul builds upon many of the excellent packages in the Julia ecosystem. Here are a few of them, and what they are used for:
 - [ForwardDiff.jl](https://github.com/JuliaDiff/ForwardDiff.jl) implements the Dual number class used throughout the code
+- [SparsityTracing.jl](https://github.com/PALEOtoolkit/SparsityTracing.jl/) provides sparsity detection inside Jutul
 - [Krylov.jl](https://github.com/JuliaSmoothOptimizers/Krylov.jl) provides the iterative linear solvers
 - [ILUZero.jl](https://github.com/mcovalt/ILUZero.jl/blob/master/src/ILUZero.jl) for ILU(0) preconditioners
 - [AlgebraicMultigrid.jl](https://github.com/JuliaLinearAlgebra/AlgebraicMultigrid.jl) for AMG preconditioners
 - [Tullio.jl](https://github.com/mcabbott/Tullio.jl) for automatically optimized loops and [Polyester.jl]() for lightweight threads
 - [TimerOutputs.jl](https://github.com/KristofferC/TimerOutputs.jl) and [ProgressMeter.jl](https://github.com/timholy/ProgressMeter.jl) gives nice output to terminal 
-- [Makie.jl](https://makie.juliaplots.org/) is used for the visualization features
+- [Makie.jl](https://makie.juliaplots.org/) is used for the visualization features found in [JutulViz.jl](https://github.com/sintefmath/JutulViz.jl)
 - [MultiComponentFlash.jl](https://github.com/moyner/MultiComponentFlash.jl) provides many of the compositional features
 
 ...and many more, both directly in the Project.toml file and indirectly!
 
 ## Getting started
-Install [Julia](https://julialang.org/) and add the package to your environment of choice.
+Install [Julia](https://julialang.org/) and add the package to your environment of choice (see details below).
 
-Here is a self-contained example that demonstrates a conceptual multiphase flow simulation of CO2 injection with multisegment wells. To run the example, you need the following packages: `Jutul` for the grid structure, `JutulDarcy` for the reservoir simulator and `Plots` for the simple plotting used here. To add these, you can run the following in Julia
+Here is a self-contained example that demonstrates a conceptual multiphase flow simulation of CO2 injection with multisegment wells. To run the example, you need the following packages: `Jutul` for the grid structure, `JutulDarcy` for the reservoir simulator and `Plots` for the simple plotting used here. To add these, you can run the following in Julia:
 ```julia
 using Pkg
 Pkg.add("Jutul")
 Pkg.add("JutulDarcy")
 Pkg.add("Plots")
 ```
-Once the packages are added, you can run this example, for example by saving it in a file `example.jl` and running `include("example.jl")`. Note that the first time run will take some time due to compilation, but subsequent runs in the same session should be quick.
+Once the packages are added and precompiled, you can run this example, for example by saving it in a file `example.jl` and running `include("example.jl")`. Note that the first time run will take some time due to compilation, but subsequent runs in the same session should be quick.
 ```julia
 using Jutul, JutulDarcy, Plots
 nx = ny = 10

src/linsolve.jl

@@ -1,3 +1,23 @@
+"""
+    reservoir_linsolve(model; <keyword arguments>)
+
+Set up iterative linear solver for a reservoir model from [`setup_reservoir_model`](@ref).
+
+# Arguments
+- `model`: Reservoir model that will linearize the equations for the linear solver
+- `precond=:cpr`: Preconditioner type to use: Either :cpr (Constrained-Pressure-Residual) or :ilu0 (block-incomplete-LU) (no effect if `solver = :direct`).
+- `v=0`: verbosity (can lead to a large amount of output)
+- `solver=:bicgstab`: the symbol of a Krylov.jl solver (typically :gmres or :bicgstab)
+- `update_interval=:once`: how often the CPR AMG hierarchy is reconstructed (:once, :iteration, :ministep, :step)
+- `update_interval_partial=:iteration`: how often the pressure system is updated in CPR
+- `max_coarse`: max size of coarse level if using AMG
+- `cpr_type=nothing`: type of CPR (`:true_impes`, `:quasi_impes` or `nothing` for automatic)
+- `partial_update=true`: perform partial update of CPR preconditioner outside of AMG update (see above)
+- `rtol=1e-3`: relative tolerance for the linear solver
+- `max_iterations=100`: limit for linear solver iterations
+
+Additional keywords are passed onto the linear solver constructor.
+"""
 function reservoir_linsolve(model,  precond = :cpr;
                                     rtol = nothing,
                                     v = 0,

src/mrst_input.jl

@@ -1005,6 +1005,10 @@ Simulate a MRST case from `file_name` as exported by `writeJutulInput` in MRST.
 - `output_path = nothing`: Directory for output files. Files will be written under this directory. Defaults to the folder of `file_name`.
 - `write_mrst = true`: Write MRST compatible output after completed simulation that can be read by `readJutulOutput` in MRST.
 - `backend=:csc`: choice of backend for linear systems. :csc for default Julia sparse, :csr for experimental parallel CSR.
+- `verbose=true`: print some extra information specific to this routine upon calling
+- `nthreads=Threads.nthreads()`: number of threads to use
+- `linear_solver=:bicgstab`: name of Krylov.jl solver to use, or :direct (for small cases only)
+- `info_level=0`: standard Jutul info_level. 0 for minimal printing, -1 for no printing, 1-5 for various levels of verbosity
 
 Additional input arguments are passed onto [`setup_reservoir_simulator`](@ref) and [`simulator_config`](@ref) if applicable.
 """

src/utils.jl

@@ -6,6 +6,14 @@ reservoir_storage(model, storage) = storage
 reservoir_storage(model::MultiModel, storage) = storage.Reservoir
 
 export setup_reservoir_model
+"""
+    setup_reservoir_model(reservoir, system; wells = [], <keyword arguments>)
+    setup_reservoir_model(reservoir, system; wells = [], context = DefaultContext(), reservoir_context = nothing, backend = :csc, <keyword arguments>)
+
+Set up a reservoir `MultiModel` for a given reservoir `SimulationModel` and an optional vector of wells.
+
+The routine automatically sets up a facility and couples the wells with the reservoir and that facility.
+"""
 function setup_reservoir_model(reservoir, system; wells = [], context = DefaultContext(), reservoir_context = nothing, backend = :csc, kwarg...)
     # List of models (order matters)
     models = OrderedDict{Symbol, Jutul.AbstractSimulationModel}()
@@ -177,6 +185,15 @@ function reservoir_multimodel(models::AbstractDict; specialize = false)
 end
 
 export setup_reservoir_state
+"""
+    setup_reservoir_state(model, <keyword arguments>)
+    # Ex: For immiscible two-phase
+    setup_reservoir_state(model, Pressure = 1e5, Saturations = [0.2, 0.8])
+
+Convenience constructor that initializes a state for a `MultiModel` set up using [`setup_reservoir_model`](@ref).
+The main convenience over [`setup_state`](@ref) is only the reservoir initialization values need be provided: wells
+are automatically initialized from the connected reservoir cells.
+"""
 function setup_reservoir_state(model; kwarg...)
     rmodel = reservoir_model(model)
     pvars = [k for k in keys(Jutul.get_primary_variables(rmodel))]
@@ -236,6 +253,11 @@ function setup_reservoir_state(model; kwarg...)
 end
 
 export setup_reservoir_forces
+"""
+    setup_reservoir_forces(model; control = nothing, limits = nothing, set_default_limits = true, <keyword arguments>)
+
+Set up driving forces for a reservoir model with wells
+"""
 function setup_reservoir_forces(model::MultiModel; control = nothing, limits = nothing, set_default_limits = true, kwarg...)
     @assert (isnothing(control) && isnothing(limits)) || haskey(model.models, :Facility) "Model must have facility."
     facility = model.models.Facility
@@ -246,6 +268,11 @@ end
 
 export full_well_outputs, well_output, well_symbols, wellgroup_symbols, available_well_targets
 
+"""
+    full_well_outputs(model, states, forces; targets = available_well_targets(model.models.Reservoir), shortname = false)
+
+Get the full set of well outputs after a simulation has occured, for plotting or other post-processing.
+"""
 function full_well_outputs(model, states, forces; targets = available_well_targets(model.models.Reservoir), shortname = false)
     out = Dict()
     if shortname
@@ -263,6 +290,11 @@ function full_well_outputs(model, states, forces; targets = available_well_targe
     return out
 end
 
+"""
+    well_output(model, states, well_symbol, forces, target = BottomHolePressureTarget)
+
+Get a specific well output from a valid operational target once a simulation is completed an `states` are available.
+"""
 function well_output(model::MultiModel, states, well_symbol, forces, target = BottomHolePressureTarget)
     n = length(states)
     d = zeros(n)