Skip to content

GeostatInversion.jl

This package provides methods for inverse analysis using parameter fields that are represented using geostatistical (stochastic) methods. Currently, two geostatistical methods are implemented. One is the Principal Component Geostatistical Approach (PCGA) proposed by Kitanidis & Lee. The other utilizes a Randomized Geostatistical Approach (RGA) that builds on PCGA.

Randomized Geostatistical Approach (RGA) references:

Two versions of PCGA are implemented in this package

  • pcgadirect, which uses full matrices and direct solvers during iterations
  • pcgalsqr, which uses low rank representations of the matrices combined with iterative solvers during iterations

The RGA method, rga, can use either of these approaches using the keyword argument. That is, by doing rga(...; pcgafunc=GeostatInversion.pcgadirect) or rga(...; pcgafunc=GeostatInversion.pcgalsqr).

GeostatInversion.jl module functions:

# GeostatInversion.getxisFunction.

Get the parameter subspace that will be explored during the inverse analysis

getxis(samplefield::Function, numfields::Int, numxis::Int, p::Int, q::Int=3, seed=nothing)
getxis(Q::Matrix, numxis::Int, p::Int, q::Int=3, seed=nothing)

Arguments:

  • samplefield : a function that takes no arguments and returns a sample of the field
  • Q : the covariance matrix of the parameter field
  • numfields : the number of fields that will be used to find the subspace
  • numxis : the dimension of the subspace
  • p : oversampling parameter when estimating the range of the covariance matrix (see Halko et al, SIAM Rev., 2011)
  • q : number of power iterations when estimating the range of the covariance matrix (see Halko et al, SIAM Rev., 2011)
  • seed : an optional seed to use when doing the randomized matrix factorization

source

# GeostatInversion.pcgadirectMethod.

Direct principal component geostatistical approach

pcgadirect(forwardmodel::Function, s0::Vector, X::Vector, xis::Array{Array{Float64, 1}, 1}, R, y::Vector; maxiters::Int=5, delta::Float64=sqrt(eps(Float64)), xtol::Float64=1e-6, callback=(s, obs_cal)->nothing)

Arguments:

  • forwardmodel : param to obs map h(s)
  • s0 : initial guess
  • X : mean of parameter prior (replace with B*X drift matrix later for p>1)
  • xis : K columns of Z = randSVDzetas(Q,K,p,q) where Q is the parameter covariance matrix
  • R : covariance of measurement error (data misfit term)
  • y : data vector
  • maxiters : maximum # of PCGA iterations
  • delta : the finite difference step size
  • xtol : convergence tolerence for the parameters
  • callback : a function of the form (params, observations)->... that is called during each iteration

source

# GeostatInversion.pcgalsqrMethod.

Iterative principal component geostatistical approach

pcgalsqr(forwardmodel::Function, s0::Vector, X::Vector, xis::Array{Array{Float64, 1}, 1}, R, y::Vector; maxiters::Int=5, delta::Float64=sqrt(eps(Float64)), xtol::Float64=1e-6)

Arguments:

  • forwardmodel : param to obs map h(s)
  • s0 : initial guess
  • X : mean of parameter prior (replace with B*X drift matrix later for p>1)
  • xis : K columns of Z = randSVDzetas(Q,K,p,q) where Q is the parameter covariance matrix
  • R : covariance of measurement error (data misfit term)
  • y : data vector
  • maxiters : maximum # of PCGA iterations
  • delta : the finite difference step size
  • xtol : convergence tolerence for the parameters

source

# GeostatInversion.rgaMethod.

Randomized (principal component) geostatistical approach

Example:

function rga(forwardmodel::Function, s0::Vector, X::Vector, xis::Array{Array{Float64, 1}, 1}, R, y::Vector, S; maxiters::Int=5, delta::Float64=sqrt(eps(Float64)), xtol::Float64=1e-6, pcgafunc=pcgadirect, callback=(s, obs_cal)->nothing)

Arguments:

  • forwardmodel : param to obs map h(s)
  • s0 : initial guess
  • X : mean of parameter prior (replace with B*X drift matrix later for p>1)
  • xis : K columns of Z = randSVDzetas(Q,K,p,q) where Q is the parameter covariance matrix
  • R : covariance of measurement error (data misfit term)
  • y : data vector
  • S : sketching matrix
  • maxiters : maximum # of PCGA iterations
  • delta : the finite difference step size
  • xtol : convergence tolerance for the parameters
  • callback : a function of the form (params, observations)->... that is called during each iteration

source

Module GeostatInversion.FDDerivatives

GeostatInversion.FDDerivatives module functions:

# GeostatInversion.FDDerivatives.makegradientFunction.

Create Gradient function

source

# GeostatInversion.FDDerivatives.makejacobianFunction.

Create Jacobian function

source

Module GeostatInversion.RandMatFact

GeostatInversion.RandMatFact module functions:

# GeostatInversion.RandMatFact.randsvdMethod.

Random SVD based on algorithm 5.1 from Halko et al.

source

Module GeostatInversion.FFTRF

GeostatInversion.FFTRF module functions:

# GeostatInversion.FFTRF.reducekMethod.

Reduce k

source