Documentation from Jutul.jl

AMG on CPU (Julia native)

Same as EntityMajorLayout, but the system is a sparse matrix where each entry is a small dense matrix.

For a test system with primary variables P, S and equations E1, E2 and two cells this will give a diagonal of length 2: [(∂E1/∂p)₁ (∂E1/∂S)₁ ; (∂E2/∂p)₁ (∂E2/∂S)₁] [(∂E1/∂p)₂ (∂E1/∂S)₂ ; (∂E2/∂p)₂ (∂E2/∂S)₂]

Entity for faces on the boundary (faces that are only connected to a single Cells)

CartesianMesh(dims, [Δ, [origin]])

Create a Cartesian mesh with dimensions specified by the Tuple dims.

Arguments

• dims::Tuple: Number of grid cells in each direction. For example, (nx, ny) will give a 2D grids with nx cells in the x-direction.
• Δ::Tuple=ones(length(dims)): Equal length to dims. First option: A Tuple of scalars where each entry is the length of each cell in that direction. For

example, specifying (Δx, Δy) for a uniform grid with each grid cell having area ofΔx*Δy. Second option: ATuple` of vectors where each entry contains the cell sizes in the direction.

• origin=zeros(length(dims)): The origin of the first corner in the grid.

Examples

Generate a uniform 3D mesh that discretizes a domain of 2 by 3 by 5 units with 3 by 5 by 2 cells:

julia> CartesianMesh((3, 5, 2), (2.0, 3.0, 5.0))
CartesianMesh (3D) with 3x5x2=30 cells

Generate a non-uniform 2D mesh:

julia> CartesianMesh((2, 3), ([1.0, 2.0], [0.1, 3.0, 2.5]))
CartesianMesh (3D) with 3x5x2=30 cells

Entity for Cells (closed volumes with averaged properties for a finite-volume solver)

CoarseMesh(G::JutulMesh, p)

Construct a coarse mesh from a given JutulMesh that can be converted to an UnstructuredMesh instance. The second argument p should be a partition Vector with one entry per cell in the original grid that assigns that cell to a coarse block. Should be one-indexed and the numbering should be sequential and contain at least one fine cell for each coarse index. This is tested by the function.

Cache that holds an AD vector/matrix together with their positions.

DataDomain(domain::JutulDomain; property1 = p1, property2 = p2, ...)

A wrapper around a domain that allows for storing of entity-associated data.

Example:

# Grid with 6 cells and 7 interior faces
g = CartesianMesh((2, 3))
d = DataDomain(g)
d[:cell_vec] = rand(6) #ok, same as:
d[:cell_vec, Cells()] = rand(6) #ok
d[:cell_vec, Faces()] = rand(6) #not ok!
d[:face_vec, Faces()] = rand(7) #ok!
# Can also add general arrays if last dimension == entity dimension
d[:cell_vec, Cells()] = rand(10, 3, 6) #ok
# Can add general data too, but needs to be specified
d[:not_on_face_or_cell, nothing] = rand(3) # also ok

Default context

DiscretizedDomain(domain, disc = nothing)

A type for a discretized domain of some other domain or mesh. May contain one or more discretizations as-needed to write equations.

Equations are grouped by entity, listing all equations and derivatives for entity 1 before proceeding to entity 2 etc.

For a test system with primary variables P, S and equations E1, E2 and two cells this will give the following ordering on the diagonal: (∂E1/∂p)₁, (∂E2/∂S)₁, (∂E1/∂p)₂, (∂E2/∂S)₂

Equations are stored sequentially in rows, derivatives of same type in columns:

For a test system with primary variables P, S and equations E1, E2 and two cells this will give the following ordering on the diagonal: (∂E1/∂p)₁, (∂E1/∂p)₂, (∂E2/∂S)₁, (∂E2/∂S)₂

Set of a variable where equations are defined

Struct that contains mappings for a set of faces that are made up of nodes and are part of cells.

Entity for Faces (intersection between pairs of Cells)

Discretization of kgradp + upwind

Abstract type for fraction variables (vector variables that sum up to unity over each entity).

By default, these are limited to the [0, 1] range through maximum_value and minimum_value default implementations.

The global set of variables

Multi-model preconditioners

Entity for half-faces (face associated with a single Cells)

HelperSimulator(model::M, T = Float64; state0 = setup_state(model), executor::E = Jutul.default_executor()) where {M, E}

Construct a helper simulator that can be used to compute the residuals and/or accumulation terms for a given type T. Useful for coupling Jutul to other solvers and types of automatic differentiation.

ILU(0) preconditioner on CPU

IndirectionMap(vals::Vector{V}, pos::Vector{Int}) where V

Create a indirection map that encodes a variable length dense vector.

pos is assumed to be a Vector{Int} of length n+1 where n is the number of dense vectors that is encoded. The vals array holds the entries for vector i in the range pos[i]:(pos[i+1]-1) for fast lookup. Indexing into the indirection map with index k will give a view into the values for vector k.

Damped Jacobi preconditioner on CPU

An AutoDiffCache is a type that holds both a set of AD values and a map into some global Jacobian.

JutulCase(model, dt = [1.0], forces = setup_forces(model); state0 = nothing, parameters = nothing, kwarg...)

Set up a structure that holds the complete specification of a simulation case.

JutulConfig(name = nothing)

A configuration object that acts like a Dict{Symbol,Any} but contains additional data to limit the valid keys and values to those added by add_option!

Abstract type for the context Jutul should execute in (matrix formats, memory allocation, etc.)

d = disc(i, Cells())

Ask discretization for entry i when discretizing some equation on the chosen entity (e.g. Cells)

Abstract type for a Jutul discretization

Abstract type for domains where equations can be defined

Super-type for all entities where JutulVariables can be defined.

Abstract type for all residual equations

Abstract type for driving forces

Abstract type for matrix layouts. A layout determines how primary variables and equations are ordered in a sparse matrix representation. Note that this is different from the matrix format itself as it concerns the ordering itself: For example, if all equations for a single cell come in sequence, or if a single equation is given for all entities before the next equation is written.

Different layouts does not change the solution of the system, but different linear solvers support different layouts.

A mesh is a type of domain that has been discretized. Abstract subtype.

Abstract type for the physical system to be solved.

Abstract type for all variables in Jutul.

A variable is associated with a JutulEntity through the associated_entity function. A variable is local to that entity, and cannot depend on other entities. Variables are used by models to define:

• primary variables: Sometimes referred to as degrees of freedom, primary unknowns or solution variables
• parameters: Static quantities that impact the solution
• secondary variables: Can be computed from a combination of other primary and secondary variables and parameters.

Full LU factorization as preconditioner (intended for smaller subsystems)

LUSolver(; reuse_memory = true, check = true, max_size = 50000)

Direct solver that calls lu directly. Direct solvers are highly accurate, but are costly in terms of memory usage and execution speed for larger systems.

MRSTWrapMesh(G, N = nothing)

Mesh that adapts an exported MRST mesh to the Jutul interface. G is assumed to be read directly from file using MAT.matread. The raw exported grid can be found under the data field.

MultiModel(models)

A model variant that is made up of many named submodels, each a fully realized SimulationModel.

models should be a NamedTuple or Dict{Symbol, JutulModel}.

An entity for something that isn't associated with an entity

Entity for Nodes (intersection between multiple Faces)

A context that uses a CSR sparse matrix format together with threads. Experimental.

Sparse Approximate Inverse preconditioner of lowest order – SPAI(0)

Single-point upwinding.

Abstract type for scalar variables (one entry per entity, e.g. pressure or temperature in each cell of a model)

SimulationModel(domain, system; <kwarg>)

Instantiate a model for a given system discretized on the domain.

SimulationModel(g::JutulMesh, system; discretization = nothing, kwarg...)

Type that defines a simulation model - everything needed to solve the discrete equations.

The minimal setup requires a JutulMesh that defines topology together with a JutulSystem that imposes physical laws.

Simulator(model; <kwarg>)

Set up a simulator object for a model that can be used by simulate!. To avoid manually instantiating the simulator, the non-mutating simulate interface can be used instead.

SparsityTracingWrapper(x::AbstractArray{T, N}) where {T, N}

Create a sparsity tracing wrapper for a numeric array. This wrapped array produces outputs that have the same value as the wrapped type, but contains a SparsityTracing seeded value with seed equal to the column index (if matrix) or linear index (if vector).

Two-point flux approximation.

TwoPointFiniteVolumeGeometry(neighbors, areas, volumes, normals, cell_centers, face_centers)

Store two-point geometry information for a given list of neighbors specified as a 2 by n matrix where n is the number of faces such that face i connectes cells N[1, i] and N[2, i].

The two-point finite-volume geometry contains the minimal set of geometry information required to compute standard finite-volume discretizations.

UnstructuredMesh(g::CartesianMesh)

Convert CartesianMesh instance to unstructured grid (3D only)

Set of a variable where variables are defined

Abstract type for vector variables (more than one entry per entity, for example saturations or displacements)

Absolute allowable change for variable during a nonlinear update.

add_option!(opts::JutulConfig, :my_cool_option, 3, "My option has this brief description")

Add an option to existing JutulConfig structure. Additional currently undocumented keyword arguments can be used to restrict valid types and values.

Update an equation so that it knows where to store its derivatives in the Jacobian representation.

Update positions of law's derivatives in global Jacobian

allocate_array_ad(v::AbstractMatrix, ...)

Convert matrix to AD matrix.

allocate_array_ad(v::AbstractVector, ...)

Convert vector to AD vector.

allocate_array_ad(n[, m]; <keyword arguments>)

Allocate vector or matrix as AD with optionally provided context and a specified non-zero on the diagonal.

Arguments

• n::Integer: number of entries in vector, or number of rows if m is given.
• m::Integer: number of rows (optional)

Keyword arguments

• npartials = 1: Number of partials derivatives to allocate for each element
• diag_pos = nothing: Indices of where to put entities on the diagonal (if any)

Other keyword arguments are passed onto get_ad_entity_scalar.

Examples:

Allocate a vector with a single partial:

julia> allocate_array_ad(2)
2-element Vector{ForwardDiff.Dual{nothing, Float64, 1}}:
 Dual{nothing}(0.0,0.0)
 Dual{nothing}(0.0,0.0)

Allocate a vector with two partials, and set the first to one:

julia> allocate_array_ad(2, diag_pos = 1, npartials = 2)
2-element Vector{ForwardDiff.Dual{nothing, Float64, 2}}:
 Dual{nothing}(0.0,1.0,0.0)
 Dual{nothing}(0.0,1.0,0.0)

Set up a matrix with two partials, where the first column has partials [1, 0] and the second [0, 1]:

julia> allocate_array_ad(2, 2, diag_pos = [1, 2], npartials = 2)
2×2 Matrix{ForwardDiff.Dual{nothing, Float64, 2}}:
 Dual{nothing}(0.0,1.0,0.0)  Dual{nothing}(0.0,1.0,0.0)
 Dual{nothing}(0.0,0.0,1.0)  Dual{nothing}(0.0,0.0,1.0)

Apply a set of forces to all equations. Equations that don't support a given force will just ignore them, thanks to the power of multiple dispatch.

apply_forces_to_equation!(diag_part, storage, model, eq, eq_s, force, time)

Update an equation with the effect of a force. The default behavior for any force we do not know about is to assume that the force does not impact this particular equation.

Create a mapped array that produces only the values when indexed.

Only useful for AD arrays, otherwise it does nothing.

Return the domain entity the equation is associated with

The entity a variable is associated with, and can hold partial derivatives with respect to.

cell_dims(g, pos)::Tuple

Get physical cell dimensions of cell with index pos for grid g.

cell_index(g, pos)

Get linear (scalar) index of mesh cell from provided IJK tuple pos.

compress_timesteps(timesteps, forces = nothing; max_step = Inf)

Compress a set of timesteps and forces to the largest possible steps that still covers the same interval and changes forces at exactly the same points in time, while being limited to a maximum size of max_step.

compress_timesteps(case::JutulCase; max_step = Inf)

Compress time steps for a Jutul case. See compress_timesteps for the general case.

compute_boundary_trans(d::DataDomain, perm)

Compute the boundary half face transmissibilities for perm. The input perm can either be the symbol of some data defined on Cells(), a vector of numbers for each cell or a matrix with number of columns equal to the number of cells.

computefacetrans(g::DataDomain, perm)

Compute face trans for the interior faces. The input perm can either be the symbol of some data defined on Cells(), a vector of numbers for each cell or a matrix with number of columns equal to the number of cells.

compute_half_face_trans(g::DataDomain, perm)

Compute half-face trans for the interior faces. The input perm can either be the symbol of some data defined on Cells(), a vector of numbers for each cell or a matrix with number of columns equal to the number of cells.

convergence_criterion(model, storage, eq, eq_s, r; dt = 1)

Get the convergence criterion values for a given equation. Can be checked against the corresponding tolerances.

Arguments

• model: model that generated the current equation.
• storage: global simulator storage.
• eq::JutulEquation: equation implementation currently being checked
• eq_s: storage for eq where values are contained.
• r: the local residual part corresponding to this model, as a matrix with column index equaling entity index

convert_from_si(value, unit_name::Union{Symbol, String})

Convert value from SI representation to the unit in unit_symbol.

Examples

julia> convert_from_si(3600.0, :hour) # Get 3600 s represented as hours
1.0

Convert a state containing variables as arrays of doubles to a state where those arrays contain the same value as Dual types. The dual type is currently taken from ForwardDiff.

convert_to_si(value, unit_name::String)

Convert value to SI representation from value in the unit given by unit_symbol.

Examples

julia> convert_to_si(1.0, :hour) # Get 1 hour represented as seconds
3600.0

Lower corner for one dimension, without any transforms applied

Give out source, target arrays of equal length for a given equation attached to the given model.

Give out I, J arrays of equal length for a given equation attached to the given model.

Number of independent primary variables / degrees of freedom per computational entity.

descalarize_primary_variable!(dest_array, model, V, var::Jutul.ScalarVariable, index)

Descalarize a primary variable, overwriting destarray at entity index. The AD status of entries in `destarray` will be retained.

descalarize_primary_variables!(state, model, V, pvars::NamedTuple = (; pairs(model.primary_variables)...), ind = eachindex(V))

Replace valeus in state by the scalarized values found in V.

dim(g)::Integer

Get the dimension of a mesh.

extra_debug_output!(report, storage, model, config, iteration, dt)

Add extra debug output to report during a nonlinear iteration.

extract_submesh(g::UnstructuredMesh, cells)

Extract a subgrid for a given mesh and a iterable of cells to keep.

Inner cell to local cell (full set)

get_1d_interpolator(xs, ys; <keyword arguments>)

Get a 1D interpolator F(x) ≈ y for a table xs, ys that by default does constant extrapolation

Arguments

• xs: sorted list of parameter points.
• ys: list of function values with equal length to xs
• method=LinearInterpolant: constructor for the interpolation. Defaults to LinearInterpolant which does simple linear interpolation.
• cap_endpoints = true: Add values so that the endpoints are capped (constant extrapolation). Otherwise, the extrapolation will match the method.
• cap_start = cap_endpoints: Fine-grained version of cap_endpoints for the start of the interval only (extrapolation for x < xs[1])
• cap_end = cap_endpoints:Fine-grained version of cap_endpoints for the end of the interval only (extrapolation for x > xs[end])

Additional keyword arguments are passed onto the interpolator constructor.

get_2d_interpolator(xs, ys, fs; method = BilinearInterpolant, cap_endpoints = true)

For xs of length nx and ys of length ny generate a 2D interpolation for values given as a nx by ny matrix. By default cap_endpoints=true, and constant extrapolation is used. Fine-grined control over extrapolation can be achieved by setting the keywords arguments cap_x = (cap_low_x, cap_high_x) and analogously for cap_y.

get_ad_entity_scalar(v::Real, npartials, diag_pos = nothing; <keyword_arguments>)

Get scalar with partial derivatives as AD instance.

Arguments

• v::Real: Value of AD variable.
• npartials: Number of partial derivatives each AD instance holds.
• diag_pos = nothing: Position(s) of where to set 1 as the partial derivative instead of zero.

Keyword arguments

• tag = nothing: Tag for AD instance. Two AD values of the different tag cannot interoperate to avoid perturbation confusion (see ForwardDiff documentation).

Get dependencies of variable when viewed as a secondary variable. Normally autogenerated with @jutul_secondary

get_diagonal_entries(eq::JutulEquation, eq_s)

Get the diagonal entries of a cache, i.e. the entries where entity type and index equals that of the governing equation.

Note: Be very careful about modifications to this array, as this is a view into the internal AD buffers and it is very easy to create inconsistent Jacobians.

get_entity_tag(basetag, entity)

Combine a base tag (which can be nothing) with a entity to get a tag that captures base tag + entity tag for use with AD initialization.

Get entries of autodiff cache. Entries are AD vectors that hold values and derivatives.

Get the entries of the main autodiff cache for an equation.

Note: This only gets the .equation field's entries.

get_mesh_entity_tag(met::JutulMesh, entity::JutulEntity, tag_group::Symbol, tag_value = missing; throw = true)

Get the indices tagged for entity in group tag_group, optionally for the specific tag_value. If ismissing(tag_value), the Dict containing the tag group will be returned.

get_parameters(model::SimulationModel)

Get the parameter definitions (as OrderedDict) for a given model.

Parameters are defined as static values in a forward simulation that combine with the primary variables to compute secondary variables and model equations.

get_primary_variable_ordered_entities(model::SimulationModel)

Get only the entities where primary variables are present, sorted by their order in the primary variables.

get_primary_variables(model::SimulationModel)

Get the primary variable definitions (as OrderedDict) for a given model.

Primary variables are sometimes referred to as solution variables or primary unknowns. The set of primary variables completely determines the state of the system together with the parameters.

get_secondary_variables(model::SimulationModel)

Get the secondary variable definitions (as OrderedDict) for a given model.

Secondary variables are variables that can be computed from the primary variables together with the parameters.

get_sparse_arguments(storage, model)

Get the [SparsePattern]@ref for the Jacobian matrix of a given simulator storage and corresponding model.

get_tstr(dT, lim = 3)

Get formatted time string of dT given in seconds, limited to lim number of units.

get_variable(model::SimulationModel, name::Symbol)

Get implementation of variable or parameter with name name for the model.

get_variables(model::SimulationModel)

Get all variable definitions (as OrderedDict) for a given model.

This is the union of get_secondary_variables and get_primary_variables.

Local cell -> global cell (full set)

Local face -> global face (full set)

Initialize context when setting up a model

initialize_extra_state_fields!(state, model::JutulModel)

Add model-dependent changing variables that need to be in state, but are never AD variables themselves (for example status flags).

initialize_storage!(storage, model::JutulModel; initialize_state0 = true)

Initialize the already allocated storage at the beginning of a simulation. Use this to e.g. set up extra stuff in state0 needed for initializing the simulation loop.

Initializer for the value of non-scalar primary variables

Local cell in full set -> inner cell (or zero)

interpolation_constant_lookup(X, constant_dx = missing)

Generate a lookup table for linear interpolation when dx is evenly spaced.

Note: Setting constant_dx=true can lead to incorrect interpolations if the data is not evenly spaced.

pth = jutul_output_path(name = missing; subfolder = "jutul", basedir = missing, create = true)

Get path for output. The final path will be found in /basedir/<subfolder/name. If subfolder=missing, the path will be set to /basedir/name instead. name will be autogenerated if not provided.

Pass the optional input create = false to avoid making the directory. To globally set the default output dir, set ENV["JUTUL_OUTPUT_PATH"]to your desiredbasedir``.

linear_timestep_selection(x, x0, x1, dt0, dt1)

Produce linear estimate of timestep dt for some value x from observed observations. If the observations have the same x or dt values, a simple scaling based on the x1 value is used.

load_balanced_endpoint(block_index, nvals, nblocks)

Endpoint for interval block_index that subdivides nvals into nblocks in a load balanced manner. This is done by adding one element to the first set of blocks whenever possible.

load_balanced_interval(b, n, m)

Create UnitRange for block b ∈ [1, m] for interval of total length n

local_ad(state::T, index::I, ad_tag::∂T) where {T, I<:Integer, ∂T}

Create localad for state for index I of AD tag of type adtag

Global cell -> local cell (full set)

Global face -> local face (full set)

local_residual_view(r_buf, model, eq, equation_offset)

Get a matrix view of the residual so that, independent of ordering, the column index corresponds to the entity index for the given equation eq starting at equation_offset in the global residual buffer r_buf.

Upper (inclusive) limit for variable.

merge_step_report_errors(data; fn = max)

Merge step reports errors of the same type using a pair wise reduction (default: max)

Lower (inclusive) limit for variable.

model_accumulation!(acc, sim::HelperSimulator, x, dt = 1.0;
    forces = setup_forces(sim.model),
    update_secondary = true,
    kwarg...
)

Compute the accumulation term into Vector acc.

model_residual!(r, sim, x, x0 = missing, dt = 1.0;
    forces = setup_forces(sim.model),
    include_accumulation = true,
    kwarg...
)

Fill in the model residual into Vector r.

model_residual(sim::HelperSimulator, x, y = missing; kwarg...)

Out of place version of model_residual!

number_of_cells(g)::Integer

Get the number of cells in a mesh.

Total number of degrees of freedom for a model, over all primary variables and all entities.

Get the number of entities (e.g. the number of cells) that the equation is defined on.

Number of entities (e.g. Cells, Faces) a variable is defined on. By default, each primary variable exists on all cells of a discretized domain

Get number of entities a cache is defined on.

Number of entities for vector stored in state (just the number of elements)

Number of entities for matrix stored in state (convention is number of columns)

Get the total number of equations on the domain of model.

Get the number of equations per entity. For example, mass balance of two components will have two equations per grid cell (= entity)

number_of_faces(g)::Integer

Get the number of faces in a mesh.

number_of_partials_per_entity(model::SimulationModel, entity::JutulEntity)

Get the number of local partial derivatives per entity in a model for a given JutulEntity. This is the sum of degrees_of_freedom_per_entity for all primary variables defined on entity.

Total number of values for a model, for a given type of variables over all entities

numerical_eltype(x::AbstractArray{T}) where T

Get the numerical eltype (i.e. the inner type of the element type that could potentially be AD)

numerical_type(::T) where T

Get the numerical eltype (i.e. the inner type of the element type that could potentially be AD). This function should be overloaded if you have a custom type that wraps a numeric/potentially AD type.

partition(N::AbstractMatrix, num_coarse, weights = ones(size(N, 2)); partitioner = MetisPartitioner(), groups = nothing, n = maximum(N), group_by_weights = false, buffer_group = true)

Partition based on neighborship (with optional groups kept contigious after partitioning)

partition_hypergraph(g, n::Int, partitioner = MetisPartitioner(); expand = true)

Partition a hypergraph from setup_partitioner_hypergraph using a given partitioner. If the optional expand parameter is set to true the result will be expanded to the full graph (i.e. where groups are not condensed).

physical_representation(x)

Get the physical representation of an object. The physical representation is usually some kind of mesh or domain that represents a physical domain.

physical_representation(x::DataDomain)

Get the underlying physical representation (domain or mesh) that is wrapped.

physical_representation(x::DiscretizedDomain)

Get the underlying physical representation (domain or mesh) that was discretized.

physical_representation(m::SimulationModel)

Get the underlying physical representation for the model (domain or mesh)

pick_next_timestep(sel::IterationTimestepSelector, sim, config, dt_prev, dT, forces, reports, current_reports, step_index, new_step)

Pick the next time-step for IterationTimestepSelector. This function uses the number of iterations from previous timesteps to estimate the relationship between the last and the new time step.

states, reports = readresults(pth; readstates = true, read_reports = true)

Read results from a given output_path provded to simulate or simulator_config.

Relative allowable change for variable during a nonlinear update. A variable with value |x| and relative limit 0.2 cannot change more than |x|*0.2.

replace_variables!(model, throw = true, varname = vardef, varname2 = vardef2)

Replace one or more variables that already exists (primary, secondary or parameters).

Arguments

• model: instance where variables is to be replaced
• varname=vardef::JutulVariables: replace variable with varname by vardef
• throw=true: throw an error if the named variable definition is not found in primary or secondary, otherwise silently return

scalarize_primary_variable(model, source_vec, var::Jutul.ScalarVariable, index)

Scalarize a primary variable. For scalars, this means getting the value itself.

scalarize_primary_variables(model, state, pvars = model.primary_variables)

Create a vector where each entry corresponds to a tuple of values that minimally defines the given variables. All variables must belong to the same type of entity. This is checked by this function.

scalarize_primary_variables!(V::Vector{T}, model, state, pvars::NamedTuple) where T

Scalarize into array. See scalarizeprimaryvariables for more details.

scalarized_primary_variable_type(model, var::Jutul.ScalarVariable)

Get the type of a scalarized numerical variable (=Float64 for variables that are already represented as scalars)

set_default_tolerances(model)

Set default tolerances for the nonlinear convergence check of the governing equations.

set_parameters!(model, parname = pardef)

Set a parameter with name varname to the definition vardef (adding if it does not exist)

set_primary_variables!(model, varname = vardef)
set_primary_variables!(model, varname1 = vardef1, varname2 = vardef2)

Set a primary variable with name varname to the definition vardef (adding if it does not exist)

set_secondary_variables!(model, varname = vardef)
set_secondary_variables!(model, varname1 = vardef1, varname2 = vardef2)

Set a secondary variable with name varname to the definition vardef (adding if it does not exist)

setup_adjoint_storage(model; state0 = setup_state(model), parameters = setup_parameters(model))

Set up storage for use with solve_adjoint_sensitivities!.

setup_forces(model::JutulModel; force_name = force_value)

Set up forces for a given model. Keyword arguments varies depending on what the model supports.

setup_parameter_optimization(model, state0, param, dt, forces, G, opt_cfg = optimization_config(model, param);
                                                        grad_type = :adjoint,
                                                        config = nothing,
                                                        print = 1,
                                                        copy_case = true,
                                                        param_obj = false,
                                                        kwarg...)

Set up function handles for optimizing the case defined by the inputs to simulate together with a per-timestep objective function G.

Generally calling either of the functions will mutate the data Dict. The options are: Fo(x) -> evaluate objective dFo(dFdx, x) -> evaluate gradient of objective, mutating dFdx (may trigger evaluation of Fo) Fand_dF(F, dFdx, x) -> evaluate F and/or dF. Value of nothing will mean that the corresponding entry is skipped.

setup_parameters(model::JutulModel; name = value)

Set up a parameter storage for a given model with values for the parameter defined in the model.

Arguments

• name=value: The name of the parameter together with the value(s) of the parameter.

A scalar (or short vector of the right size for VectorVariables) will be repeated over the entire domain, while a vector (or matrix for VectorVariables) with length (number of columns for VectorVariables) equal to the entity count (for example, number of cells for a cell variable) will be used directly.

setup_partitioner_hypergraph(N::Matrix{Int};
    num_nodes::Int = maximum(N),
    num_edges::Int = size(N, 2),
    node_weights::Vector{Int} = ones(Int, num_nodes),
    edge_weights::Vector{Int} = ones(Int, num_edges),
    groups = [Int[]]
)

Set up a hypergraph structure for a given neighborship matrix. N should be a matrix with two rows, with one pair of cells in each column. Optionally node and edge weights can be provided. If a list of groups are provided, these nodes will be accumulated together in the hypergraph.

setup_state!(state, model::JutulModel, init_values::AbstractDict = Dict())

Initialize primary variables and other state fields, given initial values as a Dict

setup_state(model::JutulModel, name1 = value1, name2 = value2)

Set up a state for a given model with values for the primary variables defined in the model. Normally all primary variables must be initialized in this way.

Arguments

• name=value: The name of the primary variable together with the value(s) used to initialize the primary variable.

A scalar (or short vector of the right size for VectorVariables) will be repeated over the entire domain, while a vector (or matrix for VectorVariables) with length (number of columns for VectorVariables) equal to the entity count (for example, number of cells for a cell variable) will be used directly.

Note: You likely want to overload [setup_state!]@ref for a custom model instead of setup_state

state, prm = setup_state_and_parameters(model, init)

Simultaneously set up state and parameters from a single init file (typically a Dict containing values that might either be initial values or parameters)

setup_storage!(storage, model::JutulModel; setup_linearized_system = true,
                                                setup_equations = true,
                                                state0 = setup_state(model),
                                                parameters = setup_parameters(model),
                                                tag = nothing,
                                                state0_ad = false,
                                                state_ad = true,
                                                kwarg...)

Allocate storage for a given model. The storage consists of all dynamic quantities used in the simulation. The default implementation allocates properties, equations and linearized system.

setup_storage(model::JutulModel; kwarg...)

Allocate storage for the model. You should overload setup_storage! if you have a custom definition.

si_unit(u::Union{String, Symbol})

Get the multiplicative SI unit conversion factor for a single unit. The return value is given so that x*si_unit(:name) will convert x to the SI representation of the unit with the given name.

Examples

julia> si_unit(:day) # Get days represented as seconds
86400.0

si_units(u1, arg...)

Get multiplicative SI unit conversion factors for multiple units simultaneously. The return value will be a Tuple of values, one for each input argument. Each input arguments can be either a String a Symbol.

Examples

julia> meter, hour = si_units(:meter, :hour)
(1.0, 3600.0)

simulate!(sim::JutulSimulator, timesteps::AbstractVector; forces = nothing,
                                                               config = nothing,
                                                               initialize = true,
                                                               restart = nothing,
                                                               state0 = nothing,
                                                               parameters = nothing,
                                                               kwarg...)

Non-allocating (or perhaps less allocating) version of simulate!.

Arguments

• initialize=true: Perform internal updates as if this is the first time

See also simulate for additional supported input arguments.

simulate(state0, model, timesteps, parameters = setup_parameters(model))
simulate(state0, model, timesteps, info_level = 3)
simulate(state0, model, timesteps; <keyword arguments>)

Simulate a set of timesteps with model for the given initial state0 and optionally specific parameters. Additional keyword arguments are passed onto simulator_config and simulate!. This interface is primarily for convenience, as all storage for the simulator is allocated upon use and discared upon return. If you want to perform multiple simulations with the same model it is advised to instead instantiate Simulator and combine it with simulate!.

Arguments

• state0::Dict: initial state, typically created using setup_state for the model in use.
• model::JutulModel: model that describes the discretized system to solve, for example SimulationModel or MultiModel.
• timesteps::AbstractVector: Vector of desired report steps. The simulator will perform time integration until sum(timesteps) is reached, providing outputs at the end of each report step.
• parameters=setup_parameters(model): Optional overrides the default parameters for the model.
• forces=nothing: Either nothing (for no forces), a single set of forces from setup_forces(model) or a Vector of such forces with equal length to timesteps.
• restart=nothing: If an integer is provided, the simulation will attempt to restart from that step. Requires that output_path is provided here or in the config.
• config=simulator_config(model): Configuration Dict that holds many fine grained settings for output, linear solver, time-steps, outputs etc.

Additional arguments are passed onto simulator_config.

See also simulate!, Simulator, SimulationModel, simulator_config.

simulate(state0, sim::JutulSimulator, timesteps::AbstractVector; parameters = nothing, kwarg...)

Simulate a set of timesteps with simulator for the given initial state0 and optionally specific parameters.

simulator_config(sim; info_level = 3, linear_solver = GenericKrylov())

Set up a simulator configuration object that can be passed onto simulate!.

There are many options available to configure a given simulator. The best way to get an overview of these possible configuration options is to instatiate the config without any arguments and inspecting the resulting table by calling simulator_config(sim) in the REPL.

solve_adjoint_sensitivities!(∇G, storage, states, state0, timesteps, G; forces = setup_forces(model))

Non-allocating version of solve_adjoint_sensitivities.

solve_adjoint_sensitivities(model, states, reports, G; extra_timing = false, state0 = setup_state(model), forces = setup_forces(model), raw_output = false, kwarg...)

Compute sensitivities of model parameter with name target for objective function G.

Solves the adjoint equations: For model equations F the gradient with respect to parameters p is ∇ₚG = Σₙ (∂Fₙ / ∂p)ᵀ λₙ where n ∈ [1, N]. Given Lagrange multipliers λₙ from the adjoint equations (∂Fₙ / ∂xₙ)ᵀ λₙ = - (∂J / ∂xₙ)ᵀ - (∂Fₙ₊₁ / ∂xₙ)ᵀ λₙ₊₁ where the last term is omitted for step n = N and G is the objective function.

solve_numerical_sensitivities(model, states, reports, G, target;
                                            forces = setup_forces(model),
                                            state0 = setup_state(model),
                                            parameters = setup_parameters(model),
                                            epsilon = 1e-8)

Compute sensitivities of model parameter with name target for objective function G.

This method uses numerical perturbation and is primarily intended for testing of solve_adjoint_sensitivities.

solve_timestep!(sim, dT, forces, max_its, config; <keyword arguments>)

Internal function for solving a single time-step with fixed driving forces.

Arguments

• sim: Simulator instance.
• dT: time-step to be solved
• forces: Driving forces for the time-step
• max_its: Maximum number of steps/Newton iterations.
• config: Configuration for solver (typically output from simulator_config).

Note: This function is exported for fine-grained simulation workflows. The general simulate interface is both easier to use and performs additional validation.

subvariable(var, map)

Get subvariable of Jutul variable

successful_reports(old_reports, current_reports, step_index, n = 1)

Get the n last successful solve reports from all previous reports (old_reports) and the current ministep set.

Synchronize backend after allocations.

Some backends may require notification that storage has been allocated.

tpfv_geometry(g)

Generate two-point finite-volume geometry for a given grid, if supported.

See also TwoPointFiniteVolumeGeometry.

Transfer v to the representation expected by a given context.

For the defalt context, the transfer function does nothing. For other context such as the CUDA version, it may convert integers and floats to other types (e.g. Float32) and Arrays to CuArrays.

You will likely have to implement some transfer operators for your own types if you want to simulate with a non-default context.

Two-point potential drop with gravity (generic)

unsafe_reinterpret(Vt, v, n)

Unsafely reinterpret v as a n length vector of value type Vt

Update equation based on currently stored properties

update_equations!(storage, model, dt = nothing)

Update the governing equations using the current set of primary variables, parameters and secondary variables. Does not fill linearized system.

update_equations_and_apply_forces!(storage, model, dt, forces; time = NaN)

Update the model equations and apply boundary conditions and forces. Does not fill linearized system.

update_linearized_system!(storage, model::JutulModel; <keyword arguments>)

Update the linearized system with the current set of equations.

Update a linearized system based on the values and derivatives in the equation.

Update a secondary variable. Normally autogenerated with @jutul_secondary

update_state_dependents!(storage, model, dt, forces; time = NaN, update_secondary = true)

Perform updates of everything that depends on the state: A full linearization for the current primary variables.

This includes properties, governing equations and the linearized system itself.

update_values!(x, dx)

Replace values (for non-Real types, direct assignment)

update_values!(x, dx)

Replace values of x in-place by y, leaving x with the values of y and the partials of x.

value(d::Dict)

Call value on all elements of some Dict.

Take value of AD.

Number of values held by a primary variable. Normally this is equal to the number of degrees of freedom, but some special primary variables are most conveniently defined by having N values and N-1 independent variables.

Define a "typical" numerical value for a variable to scale the linear system entries.

Designate the function as updating a secondary variable.

A generic evaluator is then defined, together with a function for getting the dependencies of that function upon the state. This is most easily documented with an example. If we define the following function annotated with the macro when updating the array containing the values of MyVarType realized for some model:

@jutul_secondary function some_fn!(target, var::MyVarType, model, a, b, c, ix)
    for i in ix
        target[i] = a[i] + b[i] / c[i]
    end
end

The purpose of the macro is to translate this into two functions. The first defines for the dependencies of the function with respect to the fields of the state (primary variables, secondary variables and parameters):

function get_dependencies(var::MyVarType, model)
   return (:a, :b, :c)
end

The second function defines a generic version that takes in state, and automatically expands the set of dependencies into getfield calls.

function update_secondary_variable!(array_target, var::MyVarType, model, state, ix)
    some_fn!(array_target, var, model, state.a, state.b, state.c, ix)
end

Note that the input names of arguments 4 to end-1 matter, as these will be fetched from state, exactly as written.