diff --git a/Project.toml b/Project.toml index 05038d7f..9981edcc 100644 --- a/Project.toml +++ b/Project.toml @@ -8,6 +8,7 @@ AlgebraicMultigrid = "2169fc97-5a83-5252-b627-83903c6c433c" DataStructures = "864edb3b-99cc-5e75-8d2d-829cb0a9cfe8" Dates = "ade2ca70-3891-5945-98fb-dc099432e06a" DelimitedFiles = "8bb1440f-4735-579b-a4ab-409b98df4dab" +DocStringExtensions = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae" ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210" Jutul = "2b460a1a-8a2b-45b2-b125-b5c536396eb9" Krylov = "ba0b0d4f-ebba-5204-a429-3ac8c609bfb7" @@ -42,22 +43,23 @@ AlgebraicMultigrid = "0.5.1, 0.6.0" DataStructures = "0.18.13" Dates = "1" DelimitedFiles = "1.9.1" +DocStringExtensions = "0.9.3" ForwardDiff = "0.10.30" +GLMakie = "0.9.2" HYPRE = "1.4.0" Jutul = "0.2.18" Krylov = "0.9.1" -LoopVectorization = "0.12.115" LinearAlgebra = "1" +LoopVectorization = "0.12.115" MAT = "0.10.3" MultiComponentFlash = "1.1.3" -GLMakie = "0.9.2" OrderedCollections = "1.6.2" Parsers = "2.7.1" PartitionedArrays = "0.3.2" Polyester = "0.6.11, 0.7.3" PrecompileTools = "1.0.1" -StaticArrays = "1.4.4" SparseArrays = "1" +StaticArrays = "1.4.4" Statistics = "1" TestItems = "0.1.0" TimerOutputs = "0.5.19" diff --git a/README.md b/README.md index 0aa13432..e73f1f02 100644 --- a/README.md +++ b/README.md @@ -3,7 +3,7 @@ [![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) -[![Jutul Darcy logo](docs/src/assets/logo_wide.png)](https://sintefmath.github.io/JutulDarcy.jl/dev/) +[![Jutul Darcy logo](https://github.com/sintefmath/JutulDarcy.jl/raw/main/docs/src/assets/logo_wide.png)](https://sintefmath.github.io/JutulDarcy.jl/dev/) # Reservoir simulation in Julia JutulDarcy.jl: 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 [Applied Computational Science group](https://www.sintef.no/en/digital/departments-new/applied-mathematics/applied-computational-sciences/) at [SINTEF Digital](https://www.sintef.no/en/digital/). diff --git a/docs/Project.toml b/docs/Project.toml index de26df83..87552d64 100644 --- a/docs/Project.toml +++ b/docs/Project.toml @@ -1,9 +1,11 @@ [deps] +CairoMakie = "13f3f980-e62b-5c42-98c6-ff1f3baf88f0" +DocStringExtensions = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae" Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4" -JutulDarcy = "82210473-ab04-4dce-b31b-11573c4f8e0a" +DocumenterCitations = "daee34ce-89f3-4625-b898-19384cb65244" Jutul = "2b460a1a-8a2b-45b2-b125-b5c536396eb9" -Literate = "98b081ad-f1c9-55d3-8b20-4c87d4299306" -CairoMakie = "13f3f980-e62b-5c42-98c6-ff1f3baf88f0" +JutulDarcy = "82210473-ab04-4dce-b31b-11573c4f8e0a" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" -Optim = "429524aa-4258-5aef-a3af-852621145aeb" +Literate = "98b081ad-f1c9-55d3-8b20-4c87d4299306" MultiComponentFlash = "35e5bd01-9722-4017-9deb-64a5d32478ff" +Optim = "429524aa-4258-5aef-a3af-852621145aeb" diff --git a/docs/make.jl b/docs/make.jl index ef1d07fd..9c6244b2 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -3,8 +3,12 @@ using Jutul using Literate using Documenter +using DocumenterCitations +bib = CitationBibliography(joinpath(@__DIR__, "src", "refs.bib")) + function build_jutul_darcy_docs(build_format = nothing; build_examples = true) - DocMeta.setdocmeta!(JutulDarcy, :DocTestSetup, :(using JutulDarcy); recursive=true) + DocMeta.setdocmeta!(JutulDarcy, :DocTestSetup, :(using JutulDarcy; using Jutul); recursive=true) + DocMeta.setdocmeta!(Jutul, :DocTestSetup, :(using Jutul); recursive=true) ## Literate pass # Base directory @@ -42,31 +46,56 @@ function build_jutul_darcy_docs(build_format = nothing; build_examples = true) prettyurls=get(ENV, "CI", "false") == "true", canonical="https://sintefmath.github.io/JutulDarcy.jl", edit_link="main", - assets=String[], + assets=String["assets/citations.css"], ) end makedocs(; - modules=[JutulDarcy], + modules=[JutulDarcy, Jutul], authors="Olav Møyner and contributors", repo="https://github.com/sintefmath/JutulDarcy.jl/blob/{commit}{path}#{line}", sitename="JutulDarcy.jl", warnonly = true, + plugins=[bib], format=build_format, pages=[ "Home" => "index.md", - "Examples" => examples_markdown, - "Usage" => [ - "Supported physical systems" =>"usage/systems.md", - "Solving the equations" => "usage/solution.md" + "Manual" => [ + "Getting started" =>"man/intro.md", + "High-level functions" =>"man/highlevel.md", + "Input files" =>"man/input_files.md", + "Driving forces" =>"man/forces.md", + "Supported physical systems" =>"man/systems.md", + "Wells and controls" =>"man/wells.md", + "Solving the equations" => "man/solution.md", + "Primary variables" => "man/primary.md", + "Secondary variables (properties)" => "man/secondary.md", + "Parameters" => "man/parameters.md", + "Visualization" =>"man/plotting.md", ], - "Internals" => "internals.md" + "Advanced usage" => [ + "Parallel solves with MPI" =>"man/mpi.md", + "JutulDarcy.jl compiled app" =>"man/compiled.md" + ], + "Examples" => examples_markdown, + "Reference" => [ + # "Internals" => "ref/internals.md", + "Jutul functions" => "ref/jutul.md" + ], + "Additional information "=> [ + "References" => "extras/refs.md", + "FAQ" => "extras/faq.md" + ] ], ) deploydocs(; - repo="github.com/sintefmath/JutulDarcy.jl", + repo="github.com/sintefmath/JutulDarcy.jl.git", devbranch="main", ) end ## build_jutul_darcy_docs() + +# ```@autodocs +# Modules = [JutulDarcy] +# ``` diff --git a/docs/src/assets/citations.css b/docs/src/assets/citations.css new file mode 100644 index 00000000..e393fbfd --- /dev/null +++ b/docs/src/assets/citations.css @@ -0,0 +1,18 @@ +.citation dl { + display: grid; + grid-template-columns: max-content auto; } + .citation dt { + grid-column-start: 1; } + .citation dd { + grid-column-start: 2; + margin-bottom: 0.75em; } + .citation ul { + padding: 0 0 2.25em 0; + margin: 0; + list-style: none;} + .citation ul li { + text-indent: -2.25em; + margin: 0.33em 0.5em 0.5em 2.25em;} + .citation ol li { + padding-left:0.75em; +} diff --git a/docs/src/extras/faq.md b/docs/src/extras/faq.md new file mode 100644 index 00000000..f6dca166 --- /dev/null +++ b/docs/src/extras/faq.md @@ -0,0 +1,8 @@ +# Frequently asked questions + +!!! info "Note about units" + JutulDarcy does currently not make us of conversion factors or explicit + units can in principle use any consistent unit system. Some default scaling + of variables assume that the magnitude pressures and velocities roughly + match that of strict SI (e.g. Pascals and cubic meters per second). These + scaling factors are primarily used when iterative linear solvers are used. diff --git a/docs/src/extras/refs.md b/docs/src/extras/refs.md new file mode 100644 index 00000000..78f29bc4 --- /dev/null +++ b/docs/src/extras/refs.md @@ -0,0 +1,4 @@ +# References + +```@bibliography +``` diff --git a/docs/src/index.md b/docs/src/index.md index 4b0eb615..2c8ebeb0 100644 --- a/docs/src/index.md +++ b/docs/src/index.md @@ -1,3 +1,5 @@ +# JutulDarcy.jl - reservoir simulation and porous media flow in Julia + ```@meta CurrentModule = JutulDarcy ``` @@ -9,21 +11,11 @@ DocTestSetup = quote end ``` -# JutulDarcy - Documentation for [JutulDarcy.jl](https://github.com/sintefmath/JutulDarcy.jl). The documentation is currently limited to docstrings and a series of examples. The examples are sorted by complexity. We suggest you start with [Gravity segregation example](Gravity segregation example). -!!! info "Note about units" - JutulDarcy does currently not make us of conversion factors or explicit - units can in principle use any consistent unit system. Some default scaling - of variables assume that the magnitude pressures and velocities roughly - match that of strict SI (e.g. Pascals and cubic meters per second). These - scaling factors are primarily used when iterative linear solvers are used. - JutulDarcy builds upon the general features found in [Jutul.jl](https://github.com/sintefmath/Jutul.jl). You may also find it useful to look at the [Jutul.jl documentation](https://sintefmath.github.io/Jutul.jl/dev/). -## Reading input files -It is also possible to read cases that have been set up in MRST (see [setup_case_from_mrst](@ref) and [simulate_mrst_case](@ref)) or from .DATA files (see [parse_data_file](@ref) and [simulate_data_file](@ref)) +## Package Features ```@index ``` diff --git a/docs/src/internals.md b/docs/src/internals.md deleted file mode 100644 index aeacb664..00000000 --- a/docs/src/internals.md +++ /dev/null @@ -1,3 +0,0 @@ -```@autodocs -Modules = [JutulDarcy] -``` diff --git a/docs/src/man/compiled.md b/docs/src/man/compiled.md new file mode 100644 index 00000000..e69de29b diff --git a/docs/src/man/forces.md b/docs/src/man/forces.md new file mode 100644 index 00000000..a55ea751 --- /dev/null +++ b/docs/src/man/forces.md @@ -0,0 +1,13 @@ +# Driving forces + +## Source terms + +```@docs +SourceTerm +``` + +## Boundary conditions + +```@docs +FlowBoundaryCondition +``` diff --git a/docs/src/man/highlevel.md b/docs/src/man/highlevel.md new file mode 100644 index 00000000..1456b8bc --- /dev/null +++ b/docs/src/man/highlevel.md @@ -0,0 +1,43 @@ +# High-level API + +## Setup + +### Meshes + +### Reservoir + +```@docs +reservoir_domain +``` + +### Wells + +```@docs +setup_well +setup_vertical_well +``` + +### Model + +```@docs +setup_reservoir_model +``` + +### Initial state + +```@docs +setup_reservoir_state +``` + +#### Why is initialization needed? + +#### Simple initialization + +#### Hydrostatic equilibriation + +## Simulation + +```@docs +simulate_reservoir +setup_reservoir_simulator +``` diff --git a/docs/src/man/input_files.md b/docs/src/man/input_files.md new file mode 100644 index 00000000..b834418d --- /dev/null +++ b/docs/src/man/input_files.md @@ -0,0 +1,34 @@ + +# Reading input files + +It is also possible to read cases that have been set up in MRST (see [`setup_case_from_mrst`](@ref) and [`simulate_mrst_case`](@ref)) or from .DATA files (see [`parse_data_file`](@ref) and [`simulate_data_file`](@ref)) + +## MAT-files from the Matlab Reservoir Simulation Toolbox (MRST) + +### Simulation of .MAT files + +```@docs +setup_case_from_mrst +simulate_mrst_case +``` + +### MRST-specific types + +```@docs +Jutul.MRSTWrapMesh +``` + +## DATA-files from commercial reservoir modelling software + +### Parsers + +```@docs +parse_data_file +parse_grdecl_file +``` + +### Simulation of .DATA files + +```@docs +simulate_data_file +``` diff --git a/docs/src/man/intro.md b/docs/src/man/intro.md new file mode 100644 index 00000000..095b0504 --- /dev/null +++ b/docs/src/man/intro.md @@ -0,0 +1,3 @@ +# Getting started + +## A short motivational example diff --git a/docs/src/man/models.md b/docs/src/man/models.md new file mode 100644 index 00000000..e69de29b diff --git a/docs/src/man/mpi.md b/docs/src/man/mpi.md new file mode 100644 index 00000000..3c67beb8 --- /dev/null +++ b/docs/src/man/mpi.md @@ -0,0 +1,50 @@ +# Parallel solution using MPI.jl, PartitionedArrays.jl and HYPRE.jl + +JutulDarcy can use threads by default, but advanced options can improve performance significantly for larger models. + +## Overview of parallel support + +There are two main ways of exploiting multiple cores in Jutul/JutulDarcy: Threads are automatically used for assembly and can be used for parts of the linear solve. If you require the best performance, you have to go to MPI where the linear solvers can use a parallel [BoomerAMG preconditioner](https://hypre.readthedocs.io/en/latest/solvers-boomeramg.html) via [HYPRE.jl](https://github.com/fredrikekre/HYPRE.jl). + +### MPI parallelization + +MPI parallelizes all aspects of the solver using domain decomposition and allows a simulation to be divided between multiple nodes in e.g. a supercomputer. It is significantly more cumbersome to use than standard simulations as the program must be launched in MPI mode. This is typically a non-interactive process where you launch your MPI processes and once they complete the simulation the result is available on disk. + +### Thread parallelization + +JutulDarcy also supports threads. By defualt, this only parallelizes property evaluations and assembly of the linear system. For many problems, the linear solve is the limiting factor for performance.Using threads is automatic if you start Julia with multiple threads. + +An experimental thread-parallel backend for matrices and linear algebra can be enabled by setting `backend=:csr` in the call to [`setup_reservoir_model`](@ref). This backend provides additional features such as a parallel zero-overlap ILU(0) implementation and parallel apply for AMG, but these features are still work in progress. + +### Mixed-mode parallelism + +You can mix the two approaches: Adding multiple threads to each MPI process can use threads to speed up assembly and property evaluations. + +## Solving with MPI in practice + +### Setting up the environment + +Tou will have to set up an environment with the following packages under Julia 1.9+: +`PartitionedArrays`, `MPI`, `JutulDarcy` and `HYPRE`. This is generally the best performing solver setup available, even if you are working in a shared memory environment. + +### Writing the script + +Write your script as usual. The following options must then be set: + +- [`setup_reservoir_model`](@ref) should have the extra keyword argument `split_wells=true`. We also recommend `backend=:csr` for the best performance. +- [`simulate_reservoir`](@ref) or [`setup_reservoir_simulator`](@ref) should get the optional argument `mode = :mpi` + +You must then run the file using the approprioate `mpiexec` as described in the MPI.jl documentation. Specialized functions will be called by `simulate_reservoir` when this is the case. We document them here, even if we recommend using the high level version of this interface: + +```@docs +setup_reservoir_simulator_parray +simulate_reservoir_parray +``` + +### Limitations for running in MPI + +!!! note + You should be familiar with the MPI programming model to use this feature. See [MPI.jl](https://juliaparallel.org/MPI.jl/stable/) for more details, and how MPI is handled in Julia specifically. + +!!! note + MPI consolidates results by writing files to disk. Unless you have a plan to work with the distributed states in-memory returned by the `simulate!` call, it is best to specify a `output_path` optional argument to [`setup_reservoir_simulator`](@ref). After the simulation, that folder will contain output just as if you had run the case in serial. diff --git a/docs/src/man/parameters.md b/docs/src/man/parameters.md new file mode 100644 index 00000000..3ae9d0e3 --- /dev/null +++ b/docs/src/man/parameters.md @@ -0,0 +1,34 @@ +# Parameters + +## General + +```@docs +FluidVolume +``` + +## Reservoir parameters + +### Transmissibility + +```@docs +JutulDarcy.Transmissibilities +``` + +```@docs +reservoir_transmissibility +``` + +### Other + +```@docs +TwoPointGravityDifference +ConnateWater +EndPointScalingCoefficients +``` + +## Well parameters + +```@docs +WellIndices +PerforationGravityDifference +``` diff --git a/docs/src/man/plotting.md b/docs/src/man/plotting.md new file mode 100644 index 00000000..4206fc36 --- /dev/null +++ b/docs/src/man/plotting.md @@ -0,0 +1,9 @@ + +# Plotting and visualization + +```@docs +plot_reservoir +plot_well_results +JutulDarcy.plot_reservoir_simulation_result +Jutul.plot_cell_data +``` diff --git a/docs/src/man/primary.md b/docs/src/man/primary.md new file mode 100644 index 00000000..2284be09 --- /dev/null +++ b/docs/src/man/primary.md @@ -0,0 +1,41 @@ +# Primary variables + +## Fluid systems + +### General + +```@docs +Pressure +ImmiscibleSaturation +``` + +### Immiscible flow + +```@docs +Saturations +``` + +### Black-oil flow + +```@docs +BlackOilUnknown +BlackOilX +``` + +### Compositional flow + +```@docs +OverallMoleFractions +``` + +## Wells + +```@docs +TotalMassFlux +``` + +## WellGroup / Facility + +```@docs +TotalSurfaceMassRate +``` diff --git a/docs/src/man/secondary.md b/docs/src/man/secondary.md new file mode 100644 index 00000000..6da29c2e --- /dev/null +++ b/docs/src/man/secondary.md @@ -0,0 +1,62 @@ +# Secondary variables (properties) + +## Fluid systems + +### General + +#### Relative permeabilities + +```@docs +BrooksCoreyRelativePermeabilities +RelativePermeabilities +JutulDarcy.ReservoirRelativePermeabilities +``` + +```@docs +PhaseRelativePermeability +``` + +```@docs +EndPointScalingCoefficients +``` + +#### Phase viscosities + +```@docs +DeckPhaseViscosities +``` + +```@docs +PhaseMassDensities +``` + +### Immiscible flow + +#### Phase densities + +#### Shrinkage factors + +```@docs +DeckShrinkageFactors +``` + +```@docs +ConstantCompressibilityDensities +``` + +### Black-oil flow + +```@docs +``` + +### Compositional flow + +```@docs +PhaseMassFractions +``` + +## Wells + +```@docs +TotalMass +``` diff --git a/docs/src/usage/solution.md b/docs/src/man/solution.md similarity index 71% rename from docs/src/usage/solution.md rename to docs/src/man/solution.md index 23ed5d70..f1497e19 100644 --- a/docs/src/usage/solution.md +++ b/docs/src/man/solution.md @@ -1,10 +1,15 @@ # Solving the equations + By default, Jutul solves a system as a fully-coupled implicit system of equations discretized with a two-point flux approximation with single-point upwind. + ## Newton's method + The standard way of solving a system of non-linear equations is by Newton's method (also known as Newton-Raphson's method). A quick recap: For a vector valued residual ``\mathbf{r}(x)`` of the primary variable vector ``\mathbf{x}`` we can defined a Newton update: + ```math \mathbf{x}^{k+1} = \mathbf{r}^k - J^{-1} \mathbf{r}(\mathbf{x}^k), \quad J_{ij} = \frac{\partial \mathbf{r}_i}{\partial \mathbf{x}_j}. ``` + JutulDarcy solves systems that generally have both non-smooth behavior and physical constraints on the values for ``\textbf{x}``. For that reason, we modify Newton's method slightly: ```math @@ -18,17 +23,28 @@ Here, ``\omega`` is a function that limits the variables so that they do not cha ``` Starting with ``\mathbf{x}^0``as some initial guess taken from the previous time-step, we can solve the system by iterating upon this loop. + ## Linear solvers and linear systems + For most practical applications it is not feasible or efficient to invert the Jacobian. JutulDarcy uses preconditioned iterative solvers by default, but it is possible to use direct solvers as well when working with smaller models. The high level interface for setting up a reservoir model [`setup_reservoir_model`](@ref) has an optional `block_backend=true` keyword argument that determines the matrix format, and consequently the linear solver type to be used. ### Direct solvers + If `block_backend` is set to `false`, Jutul will assemble into the standard Julia CSC sparse matrix with `Float64` elements and Julia's default direct solver will be used. It is also possible to use other Julia solvers on this system, but the default preconditioners assume that block backend is enabled. -### Iterative solver +### Iterative solver + If `block_backend` is set to `true`, Jutul will by default use a constrained-pressure residual (CPR) preconditioner for BiCGStab. Jutul relies on [Krylov.jl](https://github.com/JuliaSmoothOptimizers/Krylov.jl) for iterative solvers. The main function that selects the linear solver is [`reservoir_linsolve`](@ref) that allows for the selection of different preconditioners and linear solvers. +```@docs +reservoir_linsolve +Jutul.GenericKrylov +``` + #### Single model (only porous medium) + If the model is a single model (e.g. only a reservoir) the matrix format is a block-CSC matrix that combines Julia's builtin sparse matrix format with statically sized elements from the [StaticArrays.jl](https://github.com/JuliaArrays/StaticArrays.jl) package. If we consider the two-phase immiscible system from [Multi-phase, immiscible flow](@ref) we have a pair of equations ``R_n, R_w`` together with the corresponding primary variables pressure and first saturation ``p, S_n`` defined for all ``N_c`` cells. Let us simplify the notation a bit so that the subscripts of the primary variables are ``p, s`` and define a ``N_c \times N_c`` block Jacobian linear system where the entires are given by: + ```math J_{ij} = \begin{bmatrix} \left(\frac{\partial r_{n}}{\partial p}\right)_{ij} & \left(\frac{\partial r_{n}}{\partial s}\right)_{ij} \\ @@ -37,22 +53,37 @@ J_{ij} = \begin{bmatrix} J_{wp} & J_{ws} \end{bmatrix}_{ij} ``` -This block system has several advantages: - - We immediately get access to more powerful version of standard Julia preconditioners provided that all operations used are applicable for matrices and are applied in the right commutative order. For example, JutulDarcy uses the [ILUZero.jl](https://github.com/mcovalt/ILUZero.jl) package when a CSC linear system is preconditioned with incomplete LU factorization with zero fill-in. - - Sparse matrix vector products are much more efficient as less indicies need to be looked up for each element wise multiplication. - - Performing local reductions over variables is much easier when they are located in a local matrix. +This block system has several advantages: +- We immediately get access to more powerful version of standard Julia + preconditioners provided that all operations used are applicable for matrices + and are applied in the right commutative order. For example, JutulDarcy uses + the [ILUZero.jl](https://github.com/mcovalt/ILUZero.jl) package when a CSC + linear system is preconditioned with incomplete LU factorization with zero + fill-in. +- Sparse matrix vector products are much more efficient as less indicies need to + be looked up for each element wise multiplication. +- Performing local reductions over variables is much easier when they are + located in a local matrix. ##### Constrained Pressure Residual -The CPR preconditioner [`CPRPreconditioner`](@ref) is a multi-stage physics-informed preconditioner that seeks to decouple the global pressure part of the system from the local transport part. In the limits of incompressible flow without gravity it can be thought of as an elliptic / hyperbolic splitting. + +The CPR preconditioner [wallis-cpr, cao-cpr](@cite) [`CPRPreconditioner`](@ref) is a multi-stage physics-informed preconditioner that seeks to decouple the global pressure part of the system from the local transport part. In the limits of incompressible flow without gravity it can be thought of as an elliptic / hyperbolic splitting. + +```@docs +CPRPreconditioner +``` The short version of the CPR preconditioner can be motivated by our test system: + ```math r_n = \frac{\partial}{\partial t} ((1 - S_w) \rho_n \phi) + \nabla \cdot (\rho_n \vec{v}_n) - \rho_n q_n = 0,\\ r_w = \frac{\partial}{\partial t} (S_w \rho_w \phi) + \nabla \cdot (\rho_w \vec{v}_w) - \rho_w q_w = 0. ``` + For simplicity, we assume that there is no gravity, source terms, or compressibility. Each equation can then be divided by their respective densities and summed up to produce a pressure equation: + ```math r_p = \frac{\partial}{\partial t} ((1 - S_w) \phi) + \nabla \cdot \vec{v}_n + \frac{\partial}{\partial t} (S_w \phi) + \nabla \cdot \vec{v}_w \\ = \frac{\partial}{\partial t} ((S_w - S_w) \phi) + \nabla \cdot (\vec{v}_n + \vec{v}_w) \\ @@ -60,6 +91,7 @@ r_p = \frac{\partial}{\partial t} ((1 - S_w) \phi) + \nabla \cdot \vec{v}_n + \f = - \nabla \mathbf{K}(k_{rw}/\mu_w + k_{rn}/\mu_n) \nabla p \\ = - \nabla \mathbf{K}\lambda_t \nabla p = 0 ``` + The final equation is the variable coefficient Poisson equation and is referred to as the incompressible pressure equation for a porous media. We know that algebraic multigrid preconditioners (AMG) are highly efficient for linear systems made by discretizing this equation. The idea in CPR is to exploit this by constructing an approximate pressure equation that is suited for AMG inside the preconditioner. Constructing the preconditioner is done in two stages: @@ -69,15 +101,15 @@ Constructing the preconditioner is done in two stages: During the linear solve, the preconditioner is then made up of two broad stages: First, a preconditioner is applied to the pressure part (typically AMG), then the full system is preconditioned (typically ILU(0)) after the residual has been corrected by the pressure estimate: -1. Form weighted pressure residual ``r_p = \sum_i w_i r_i ``. +1. Form weighted pressure residual ``r_p = \sum_i w_i r_i``. 2. Apply pressure preconditioer ``M_p``: ``\Delta p = M_p^{-1} r_p``. 3. Correct global residual ``r^* = r - J P(\Delta p)`` where ``P`` expands the pressure update to the full system vector, with zero entries outside the pressure indices. 4. Precondition the full system ``\Delta x^* = M^{-1}r^*`` 5. Correct the global update with the pressure to obtain the final update: ``\Delta x = \Delta x^* + P(\Delta p)`` #### Multi model (porous medium with wells) -If a model is a porous medium with wells, the same preconditioners can be used, but an additional step is required to incorporate the well system. In practical terms, this means that our linearized system is expanded to multiple linear systems: +If a model is a porous medium with wells, the same preconditioners can be used, but an additional step is required to incorporate the well system. In practical terms, this means that our linearized system is expanded to multiple linear systems: ```math J \Delta \mathbf{x} = \begin{bmatrix} @@ -94,9 +126,11 @@ J \Delta \mathbf{x} = \begin{bmatrix} \mathbf{r}_w \end{bmatrix} ``` + Here, ``J_{rr}`` is the reservoir equations differentiated with respect to the reservoir primary variables, i.e. the Jacobian from the previous section. ``J_{ww}`` is the well system differentiated with respect to the well primary variables. The cross terms, ``J_{rw}``and ``J_{wr}``, are the same equations differentiated with respect to the primary variables of the other system. The well system is generally much smaller than the reservoir system and can be solved by a direct solver. We would like to reuse the block preconditioners defined for the base system. The approach we use is a [Schur complement](https://en.wikipedia.org/wiki/Schur_complement) approach to solve the full system. If we linearly eliminate the dependence of the reservoir equations on the well primary variables, we obtain the reduced system: + ```math J \Delta \mathbf{x} = \begin{bmatrix} J_{rr} - J_{rw}J_{ww}^{-1}J_{wr} & 0 \\ @@ -112,41 +146,18 @@ J \Delta \mathbf{x} = \begin{bmatrix} \mathbf{r}_w \end{bmatrix} ``` + We can then solve the system in terms of the reservoir degrees of freedom where the system is a block linear system and we already have a working preconditioner: + ```math \left(J_{rr} - J_{rw}J_{ww}^{-1}J_{wr}\right)\mathbf{x}_r = \mathbf{r}_r - J_{rw}J_{ww}^{-1}\mathbf{r}_w ``` + Once that system is solved for ``\mathbf{x}_r``, we can recover the well degrees of freedom ``\mathbf{r}_w`` directly: + ```math \mathbf{r}_w = J_{ww}^{-1}(\mathbf{r}_w - J_{wr}\mathbf{x}_r) ``` - !!! note "Efficiency of Schur complement" - Explicitly forming the matrix ``J_{rr} - J_{rw}J_{ww}^{-1}J_{wr}`` will generally lead to a lot of fill-in in the linear system. JutulDarcy instead uses the action of ``J_{rr} - J_{rw}J_{ww}^{-1}J_{wr}`` as a linear operator from [LinearOperators.jl](https://github.com/JuliaSmoothOptimizers/LinearOperators.jl). This means that we must apply the inverse of the well system every time we need to compute the residual or action of the system matrix, but fortunately performing the action of the Schur complement is inexpensive as long as ``J_{ww}`` is small and the factorization can be stored. - -## References -[Wallis, J.R. "Incomplete Gaussian Elimination as a Preconditioning for Generalized Conjugate Gradient Acceleration." Paper presented at the SPE Reservoir Simulation Symposium, San Francisco, California, November 1983](https://doi.org/10.2118/12265-MS) - -[Cao, H., Tchelepi, H. A., Wallis, J., and H. Yardumian. "Parallel Scalable Unstructured CPR-Type Linear Solver for Reservoir Simulation." Paper presented at the SPE Annual Technical Conference and Exhibition, Dallas, Texas, October 2005](https://doi.org/10.2118/96809-MS) -# Parallel simulations - -JutulDarcy can use threads by default, but advanced options can improve performance significantly for larger models. - -## Overview of parallel support -There are two main ways of exploiting multiple cores in Jutul/JutulDarcy: Threads are automatically used for assembly and can be used for parts of the linear solve. If you require the best performance, you have to go to MPI where the linear solvers can use a parallel [BoomerAMG preconditioner](https://hypre.readthedocs.io/en/latest/solvers-boomeramg.html) via [HYPRE.jl](https://github.com/fredrikekre/HYPRE.jl). -#### Shared memory -An experimental thread-parallel backend for matrices and linear algebra can be enabled by setting `backend=:csr` in the call to [`setup_reservoir_model`](@ref). This backend provides additional features such as a parallel zero-overlap ILU(0) implementation and parallel apply for AMG, but these features are still work in progress. - -#### MPI support for distributed memory -It is possible to run cases using MPI. You will have to set up an environment with the following packages under Julia 1.9+: -`PartitionedArrays`, `MPI`, `JutulDarcy` and `HYPRE`. This is generally the best performing solver setup available, even if you are working in a shared memory environment. - -Write your script as usual, but in your call to `setup_reservoir_simulator`, pass the optional argument `mode = :mpi`. You must then run the file using `mpiexec`. - -!!! note - You should be familiar with the MPI programming model to use this feature. See [MPI.jl](https://juliaparallel.org/MPI.jl/stable/) and [MPIClusterManagers.jl](https://github.com/JuliaParallel/MPIClusterManagers.jl) for more details, and how MPI is handled in Julia specifically. - -!!! note - MPI consolidates results by writing files to disk. Unless you have a plan to work with the distributed states in-memory returned by the `simulate!` call, it is best to specify a `output_path` optional argument to `setup_reservoir_simulator`. After the simulation, that folder will contain output just as if you had run the case in serial. - + Explicitly forming the matrix ``J_{rr} - J_{rw}J_{ww}^{-1}J_{wr}`` will generally lead to a lot of fill-in in the linear system. JutulDarcy instead uses the action of ``J_{rr} - J_{rw}J_{ww}^{-1}J_{wr}`` as a linear operator from [LinearOperators.jl](https://github.com/JuliaSmoothOptimizers/LinearOperators.jl). This means that we must apply the inverse of the well system every time we need to compute the residual or action of the system matrix, but fortunately performing the action of the Schur complement is inexpensive as long as ``J_{ww}`` is small and the factorization can be stored. diff --git a/docs/src/usage/systems.md b/docs/src/man/systems.md similarity index 87% rename from docs/src/usage/systems.md rename to docs/src/man/systems.md index b37ffb50..5c96bf82 100644 --- a/docs/src/usage/systems.md +++ b/docs/src/man/systems.md @@ -1,7 +1,9 @@ # Supported physical systems + JutulDarcy supports a number of different systems. These are [`JutulSystem`](https://sintefmath.github.io/Jutul.jl/dev/usage/#Jutul.JutulSystem) instances that describe a particular type of physics for porous media flow. We describe these in roughly the order of complexity that they can model. ## Summary + The general form of the flow systems we will discuss is a conservation law for ``N`` components on residual form: ```math @@ -14,17 +16,38 @@ The following table gives an overview of the available features that are describ | System | Number of phases | Number of components | ``M`` | ``V`` | |---|---|---|---|---| -| [`SinglePhaseSystem`](@ref) | 1 | 1 | ``\rho \phi `` | ``\rho \vec{v}`` | +| [`SinglePhaseSystem`](@ref) | 1 | 1 | ``\rho \phi`` | ``\rho \vec{v}`` | | [`ImmiscibleSystem`](@ref) | Any | (Any) | ``S_\alpha \rho_\alpha \phi`` | ``\rho_\alpha \vec{v}_\alpha`` | -| [`StandardBlackOilSystem`](@ref) | 2-3 | (2-3) | ``\rho_o^s(b_g S_g + R_s b_o S_o)`` | `` b_g \vec{v}_g + R_s b_g \vec{v}_o`` | +| [`StandardBlackOilSystem`](@ref) | 2-3 | (2-3) | ``\rho_o^s(b_g S_g + R_s b_o S_o)`` | ``b_g \vec{v}_g + R_s b_g \vec{v}_o`` | | [`MultiPhaseCompositionalSystemLV`](@ref) | 2-3 | Any | ``\rho_l X_i S_l + \rho_v Y_i S_v`` | ``\rho_l X_i \vec{v}_l + \rho_v Y_i \vec{v}_v`` | +### Phases +Phases are defined using specific types. Some constructors take a list of phases present in the model. Phases do not contain any data themselves and the distinction between different phases applies primarily for well controls. + +```@docs +LiquidPhase +AqueousPhase +VaporPhase +``` ### Implementation details -In the above the discrete version of ``M_i`` is implemented in the update function for [`TotalMasses`](@ref) that should by convention be named [`update_total_masses!`](@ref). The discrete component fluxes are implemented by [`component_mass_fluxes!`](@ref). The source terms are implemented by [`apply_forces_to_equation!`](@ref) for boundary conditions and sources, and [`update_cross_term_in_entity!`](@ref) for wells. We use Julia's multiple dispatch to pair the right implementation with the right physics system. + +In the above the discrete version of ``M_i`` is implemented in the update function for [`TotalMasses`](@ref) that should by convention be named [`JutulDarcy.update_total_masses!`](@ref). The discrete component fluxes are implemented by [`JutulDarcy.component_mass_fluxes!`](@ref). + +```@docs +JutulDarcy.update_total_masses! +JutulDarcy.component_mass_fluxes! +``` + +The source terms are implemented by [`Jutul.apply_forces_to_equation!`](@ref) for boundary conditions and sources, and [`Jutul.update_cross_term_in_entity!`](@ref) for wells. We use Julia's multiple dispatch to pair the right implementation with the right physics system. ## Single-phase flow + +```@docs +SinglePhaseSystem +``` + The simplest form of porous media flow is the single-phase system. ```math @@ -37,44 +60,62 @@ r(p) = \frac{\partial}{\partial t}( \rho \phi) + \nabla \cdot (\rho \vec{v}) - \ \vec{v} = - \frac{\mathbf{K}}{\mu} (\nabla p + \rho g \nabla z) ``` -Here, ``\mathbf{K}`` is a positive-definite permeability tensor, ``\mu`` the fluid viscosity, ``g`` the magnitude of gravity oriented down and ``z`` the depth. +Here, ``\mathbf{K}`` is a positive-definite permeability tensor, ``\mu`` the fluid viscosity, ``g`` the magnitude of gravity oriented down and ``z`` the depth. !!! note "Single-phase implementation" The [`SinglePhaseSystem`](@ref) is a dedicated single phase system. This is mathematically equivalent to an [`ImmiscibleSystem`](@ref) when set up with a single phase. For single phase flow, the fluid [`Pressure`](@ref) is the primary variable in each cell. The equation supports two types of compressibility: That of the fluid where density is a function ``\rho(p)`` of pressure and that of the pores where the porosity ``\phi(p)`` changes with pressure. - !!! tip JutulDarcy uses the notion of depth rather than coordinate when defining buoyancy forces. This is consistent with the convention in the literature on subsurface flow. - ## Multi-phase, immiscible flow + +```@docs +ImmiscibleSystem +``` + The flow systems immediately become more interesting if we add more phases. We can extend the above single-phase system by introducing the phase saturation of phase with label ``\alpha`` as ``S_\alpha``. The phase saturation represents the volumetric fraction of the rock void space occupied by the phase. If we consider a pair of phases ``\{n, w\}`` non-wetting and wetting we can write the system as ```math r_\alpha = \frac{\partial}{\partial t} (S_\alpha \rho_\alpha \phi) + \nabla \cdot (\rho_\alpha \vec{v}_\alpha) - \rho_\alpha q_\alpha = 0, \quad \alpha \in \{n, w\} ``` + This requires an additional closure such that the amount of saturation of all phases exactly fills the available fluid volume: + ```math S_w + S_n = 1, \quad 1 \ge S_\alpha \ge 0 \quad \alpha \in \{n, w\} ``` + This equation is local and linear in the saturations and can be eliminated to produce the classical two-equation system for two-phase flow, + ```math r_n = \frac{\partial}{\partial t} ((1 - S_w) \rho_n \phi) + \nabla \cdot (\rho_n \vec{v}_n) - \rho_n q_n = 0,\\ r_w = \frac{\partial}{\partial t} (S_w \rho_w \phi) + \nabla \cdot (\rho_w \vec{v}_w) - \rho_w q_w = 0. ``` + To complete this description we also need expressions for the phase fluxes. We use the standard multiphase extension of Darcy's law, + ```math \vec{v}_\alpha = - \mathbf{K} \frac{k_{r\alpha}}{\mu_\alpha} (\nabla p\alpha + \rho_\alpha g \nabla z) ``` + Here, we have introduced the relative permeability of the phase ``k_{r\alpha}``, an empirical relationship between the saturation and the flow rate. Relative permeability is a complex topic with many different relationships and functional forms, but we limit the discussion to monotone, non-negative functions of their respective saturations, for example a simple Brooks-Corey type of ``k_{r\alpha}(S_\alpha) = S_\alpha^2``. We have also introduced separate phase pressures ``p_\alpha`` that account for capillary pressure, e.g. ``p_w = p_n + p_c(S_w)``. +### Primary variables + !!! note "Immiscible implementation" The [`ImmiscibleSystem`](@ref) implements this system for any number of phases. The primary variables for this system is a single reference [`Pressure`](@ref) and phase [`Saturations`](@ref). As we do not solve for the volume closure equation, there is one less degree of freedom associated with the saturations than there are number of phases. ## Black-oil: Multi-phase, pseudo-compositional flow + +```@docs +StandardBlackOilSystem +``` + The black-oil equations is an extension of the immiscible description to handle limited miscibility between the phases. Originally developed for certain types of oil and gas simulation, these equations are useful when the number of components is low and tabulated values for dissolution and vaporization are available. The assumptions of the black-oil model is that the "oil" and "gas" pseudo-components have uniform composition throughout the domain. JutulDarcy supports two- and three-phase black oil flow. The difference between two and three phases amounts to an additional immiscible aqueous phase that is identical to that of the previous section. For that reason, we focus on the miscible pseudo-components: + ```math r_o = \rho_o^s \left( \frac{\partial}{\partial t}( (b_o S_o + R_v b_g (1 - S_o)) \phi) + \nabla \cdot ( b_o \vec{v}_o + R_v b_o \vec{v}_g) - q_o^s \right ) \\ r_g = \rho_g^s \left( \frac{\partial}{\partial t}( (b_g S_g + R_s b_o S_o) \phi) + \nabla \cdot ( b_g \vec{v}_g + R_s b_g \vec{v}_o) - q_g^s \right ) @@ -85,8 +126,14 @@ The model uses the notion of surface (or reference densities) ``\rho_o^s, \rho_g !!! note "Blackoil implementation" The [`StandardBlackOilSystem`](@ref) implements the black-oil equations. It is possible to run cases with and without water, with and without ``R_s`` and with and without ``R_v``. The primary variables for the most general case is the reference [`Pressure`](@ref), an [`ImmiscibleSaturation`](@ref) for the aqueous phase and the special [`BlackOilUnknown`](@ref) that will represent either ``S_o``, ``R_s`` or ``R_v`` on a cell-by-cell basis depending on what phases are present and saturated. -A full description of the black-oil equations is outside the scope of this documentation. Please see [Lie, An Introduction to Reservoir Simulation Using MATLAB/GNU Octave, Cambridge University Press, 2019](https://doi.org/10.1017/9781108591416) for more details. +A full description of the black-oil equations is outside the scope of this documentation. Please see [mrst-book-i](@cite) for more details. + ## Compositional: Multi-phase, multi-component flow + +```@docs +MultiPhaseCompositionalSystemLV +``` + The more general case of multi-component flow is often referred to as a compositional model. The typical version of this model describes the fluid as a system of ``N`` components where the phases present and fluid properties are determined by an equation-of-state. This can be highly accurate if the equation-of-state is tuned for the mixtures that are encountered, but comes at a significant computational cost as the equation-of-state must be evaluated many times. JutulDarcy implements a standard compositional model that assumes local instantaneous equilibrium and that the components are present in up to two phases with an optional immiscible phase added. This is sometimes referred to as a "simple water" or "dead water" description. By default the solvers use [MultiComponentFlash.jl](https://github.com/moyner/MultiComponentFlash.jl) to solve thermodynamic equilibrium. This package implements the generalized cubic approach and defaults to Peng-Robinson. @@ -97,10 +144,11 @@ Assume that we have two phases liquid and vapor referred to as ``l`` and ``v`` w r_i = \frac{\partial}{\partial t} \left( (\rho_l X_i S_l + \rho_v Y_i S_v) \phi \right) + \nabla \cdot (\rho_l X_i \vec{v}_l + \rho_v Y_i \vec{v}_v) - Q_i, \quad M \in \{1, \dots, M\} ``` -For additional details, please see [Chapter 8 - Compositional Simulation with the AD-OO Framework Advanced Modeling with the MATLAB Reservoir Simulation Toolbox, Møyner, 2021](https://doi.org/10.1017/9781009019781). +For additional details, please see Chapter 8 - Compositional Simulation with the AD-OO Framework in [mrst-book-ii](@cite). !!! note "Compositional implementation" The [`MultiPhaseCompositionalSystemLV`](@ref) implements the compositional model. The primary variables for the most general case is the reference [`Pressure`](@ref), an [`ImmiscibleSaturation`](@ref) for the optional immiscible phase and ``M-1`` [`OverallMoleFractions`](@ref). ## Thermal flow + Currently experimental and undocumented. See [`ThermalSystem`](@ref) if you are feeling brave. diff --git a/docs/src/man/wells.md b/docs/src/man/wells.md new file mode 100644 index 00000000..5a5b797c --- /dev/null +++ b/docs/src/man/wells.md @@ -0,0 +1,65 @@ +# Wells and controls + +## Well setup routines + +Wells can be set up using the convenience functions [`setup_well`](@ref) and [`setup_vertical_well`](@ref). These routines act on the output from [`reservoir_domain`](@ref) and can set up both types of wells. We recommend that you use these functions instead of manually calling the well constructors. + +## Types of wells + +### Simple wells + +```@docs +JutulDarcy.SimpleWell +``` + +#### Equations + +### Multisegment wells + +```@docs +MultiSegmentWell +``` + +## Well controls and limits + +### Types of well controls + +```@docs +InjectorControl +ProducerControl +DisabledControl +``` + +```@docs +JutulDarcy.replace_target +JutulDarcy.default_limits +``` + +### Types of well targets + +```@docs +BottomHolePressureTarget +SinglePhaseRateTarget +SurfaceLiquidRateTarget +SurfaceOilRateTarget +SurfaceGasRateTarget +SurfaceWaterRateTarget +TotalRateTarget +HistoricalReservoirVoidageTarget +ReservoirVoidageTarget +DisabledTarget +``` + +### Imposing limits on wells (multiple constraints) + +## Well forces + +### Changing perforations + +```@docs +PerforationMask +``` + +### Other forces + +Can use [`SourceTerm`](@ref) or [`FlowBoundaryCondition`](@ref) diff --git a/docs/src/ref/internals.md b/docs/src/ref/internals.md new file mode 100644 index 00000000..e69de29b diff --git a/docs/src/ref/jutul.md b/docs/src/ref/jutul.md new file mode 100644 index 00000000..febabc42 --- /dev/null +++ b/docs/src/ref/jutul.md @@ -0,0 +1,7 @@ +# Documentation from Jutul.jl + +JutulDarcy.jl builds upon Jutul.jl. You can use JutulDarcy.jl without knowing the inner workings of Jutul.jl, but if you want to dive under the hood these docstrings may be helpful. + +```@autodocs +Modules = [Jutul] +``` diff --git a/docs/src/refs.bib b/docs/src/refs.bib new file mode 100644 index 00000000..1e697d95 --- /dev/null +++ b/docs/src/refs.bib @@ -0,0 +1,35 @@ +@book{mrst-book-i, + title={An introduction to reservoir simulation using MATLAB/GNU Octave: User guide for the MATLAB Reservoir Simulation Toolbox (MRST)}, + author={Lie, Knut-Andreas}, + year={2019}, + doi={10.1017/9781108591416}, + url={https://www.cambridge.org/core/books/an-introduction-to-reservoir-simulation-using-matlabgnu-octave/F48C3D8C88A3F67E4D97D4E16970F894}, + publisher={Cambridge University Press} +} + +@book{mrst-book-ii, + title={Advanced modelling with the MATLAB reservoir simulation toolbox}, + author={Lie, Knut-Andreas and M{\o}yner, Olav}, + year={2021}, + doi={10.1017/9781009019781}, + url={https://www.cambridge.org/core/books/advanced-modeling-with-the-matlab-reservoir-simulation-toolbox/7AC2425C73F6F729DB88DB1A504FA1E7}, + publisher={Cambridge University Press} +} + +@inproceedings{cao-cpr, + title={Parallel scalable unstructured CPR-type linear solver for reservoir simulation}, + author={Cao, Hui and Tchelepi, Hamdi A and Wallis, J and Yardumian, H}, + booktitle={SPE Annual Technical Conference and Exhibition}, + year={2005}, + url={https://doi.org/10.2118/12265-MS}, + organization={Society of Petroleum Engineers} +} + +@inproceedings{wallis-cpr, + title={Incomplete Gaussian elimination as a preconditioning for generalized conjugate gradient acceleration}, + author={Wallis, John R}, + booktitle={SPE Reservoir Simulation Conference}, + url={https://doi.org/10.2118/12265-MS}, + year={1983}, + organization={Society of Petroleum Engineers} +} diff --git a/examples/co2_brine_2d_vertical.jl b/examples/co2_brine_2d_vertical.jl index ca883619..90d60c0d 100644 --- a/examples/co2_brine_2d_vertical.jl +++ b/examples/co2_brine_2d_vertical.jl @@ -45,7 +45,7 @@ L, V = LiquidPhase(), VaporPhase() sys = MultiPhaseCompositionalSystemLV(eos, (L, V)) model, parameters = setup_reservoir_model(res, sys, wells = [inj, prod], reference_densities = rhoS); push!(model[:Reservoir].output_variables, :Saturations) -kr = BrooksCoreyRelPerm(sys, 2.0, 0.0, 1.0) +kr = BrooksCoreyRelativePermeabilities(sys, 2.0, 0.0, 1.0) model = replace_variables!(model, RelativePermeabilities = kr) T0 = repeat([303.15*Kelvin], 1, nc) parameters[:Reservoir][:Temperature] = T0 diff --git a/examples/compositional_5components.jl b/examples/compositional_5components.jl index 2eb7812a..bf63411a 100644 --- a/examples/compositional_5components.jl +++ b/examples/compositional_5components.jl @@ -42,13 +42,12 @@ L, V = LiquidPhase(), VaporPhase() # Define system and realize on grid sys = MultiPhaseCompositionalSystemLV(eos, (L, V)) model, parameters = setup_reservoir_model(res, sys, wells = [inj, prod], reference_densities = rhoS, block_backend = true); -kr = BrooksCoreyRelPerm(sys, 2.0, 0.0, 1.0) +kr = BrooksCoreyRelativePermeabilities(sys, 2.0, 0.0, 1.0) model = replace_variables!(model, RelativePermeabilities = kr) push!(model[:Reservoir].output_variables, :Saturations) -T0 = repeat([387.45*Kelvin], 1, nc) -parameters[:Reservoir][:Temperature] = T0 +parameters[:Reservoir][:Temperature] = 387.45*Kelvin state0 = setup_reservoir_state(model, Pressure = 225*bar, OverallMoleFractions = [0.463, 0.01640, 0.20520, 0.19108, 0.12432]); dt = repeat([2.0]*day, 365) diff --git a/examples/five_spot_ensemble.jl b/examples/five_spot_ensemble.jl index 9cad3979..b155452a 100644 --- a/examples/five_spot_ensemble.jl +++ b/examples/five_spot_ensemble.jl @@ -37,7 +37,7 @@ function simulate_qfs(porosity = 0.2) model, parameters = setup_reservoir_model(domain, sys, wells = [Inj, Prod]) c = [1e-6/bar, 1e-6/bar] ρ = ConstantCompressibilityDensities(p_ref = 150*bar, density_ref = rhoS, compressibility = c) - kr = BrooksCoreyRelPerm(sys, [2.0, 2.0]) + kr = BrooksCoreyRelativePermeabilities(sys, [2.0, 2.0]) replace_variables!(model, PhaseMassDensities = ρ, RelativePermeabilities = kr); state0 = setup_reservoir_state(model, Pressure = 150*bar, Saturations = [1.0, 0.0]) diff --git a/examples/optimize_simple_bl.jl b/examples/optimize_simple_bl.jl index 37c873fb..f727ab17 100644 --- a/examples/optimize_simple_bl.jl +++ b/examples/optimize_simple_bl.jl @@ -18,7 +18,7 @@ function setup_bl(;nc = 100, time = 1.0, nstep = 100, poro = 0.1, perm = 9.8692e sys = ImmiscibleSystem((LiquidPhase(), VaporPhase())) model = SimulationModel(G, sys) model.primary_variables[:Pressure] = Pressure(minimum = -Inf, max_rel = nothing) - kr = BrooksCoreyRelPerm(sys, [2.0, 2.0]) + kr = BrooksCoreyRelativePermeabilities(sys, [2.0, 2.0]) replace_variables!(model, RelativePermeabilities = kr) tot_time = sum(tstep) diff --git a/examples/two_phase_buckley_leverett.jl b/examples/two_phase_buckley_leverett.jl index 3a8bbc16..4ae11e7a 100644 --- a/examples/two_phase_buckley_leverett.jl +++ b/examples/two_phase_buckley_leverett.jl @@ -32,7 +32,7 @@ function solve_bl(;nc = 100, time = 1.0, nstep = nc) p0 = 100*bar sys = ImmiscibleSystem((LiquidPhase(), VaporPhase())) model = SimulationModel(domain, sys) - kr = BrooksCoreyRelPerm(sys, [2.0, 2.0], [0.2, 0.2]) + kr = BrooksCoreyRelativePermeabilities(sys, [2.0, 2.0], [0.2, 0.2]) replace_variables!(model, RelativePermeabilities = kr) tot_time = sum(timesteps) pv = pore_volume(domain) diff --git a/examples/two_phase_unstable_gravity.jl b/examples/two_phase_unstable_gravity.jl index d2a6d32d..12883f35 100644 --- a/examples/two_phase_unstable_gravity.jl +++ b/examples/two_phase_unstable_gravity.jl @@ -29,7 +29,7 @@ sys = ImmiscibleSystem([L, V]) model, parameters = setup_reservoir_model(domain, sys) density = ConstantCompressibilityDensities(sys, p0, [rhoLS, rhoVS], [cl, cv]) # Replace density with a lighter pair replace_variables!(model, PhaseMassDensities = density); -kr = BrooksCoreyRelPerm(sys, [2.0, 3.0]) +kr = BrooksCoreyRelativePermeabilities(sys, [2.0, 3.0]) replace_variables!(model, RelativePermeabilities = kr) # ### Define initial saturation diff --git a/ext/JutulDarcyGLMakieExt/well_plots.jl b/ext/JutulDarcyGLMakieExt/well_plots.jl index d17a91f1..e76a7ad0 100644 --- a/ext/JutulDarcyGLMakieExt/well_plots.jl +++ b/ext/JutulDarcyGLMakieExt/well_plots.jl @@ -55,7 +55,6 @@ function well_cells_for_plot(w) return w.perforations.reservoir end -# export plot_well_results function JutulDarcy.plot_well_results(well_data::Dict, arg...; name = "Data", kwarg...) JutulDarcy.plot_well_results([well_data], arg...; names = [name], kwarg...) end diff --git a/src/InputParser/deckinput/keywords/grid.jl b/src/InputParser/deckinput/keywords/grid.jl index 9984c84e..ce03a9d3 100644 --- a/src/InputParser/deckinput/keywords/grid.jl +++ b/src/InputParser/deckinput/keywords/grid.jl @@ -97,6 +97,21 @@ function parse_keyword!(data, outer_data, units, cfg, f, v::Union{Val{:PERMX}, V parse_and_set_grid_data!(data, outer_data, units, cfg, f, k, unit = :permeability) end +function parse_keyword!(data, outer_data, units, cfg, f, v::Val{:MULTREGT}) + parser_message(cfg, outer_data, "MULTREGT", PARSER_MISSING_SUPPORT) + skip_record(f) +end + +function parse_keyword!(data, outer_data, units, cfg, f, v::Union{ + Val{:MULTX}, Val{:MULTY}, Val{:MULTZ}, + Val{Symbol("MULTX-")}, Val{Symbol("MULTY-")}, Val{Symbol("MULTZ-")} + } + ) + k = unpack_val(v) + parser_message(cfg, outer_data, "$k", PARSER_JUTULDARCY_MISSING_SUPPORT) + parse_and_set_grid_data!(data, outer_data, units, cfg, f, k, unit = :id) +end + function parse_keyword!(data, outer_data, units, cfg, f, v::Val{:MULTPV}) k = unpack_val(v) parse_and_set_grid_data!(data, outer_data, units, cfg, f, k) diff --git a/src/InputParser/deckinput/keywords/props.jl b/src/InputParser/deckinput/keywords/props.jl index ab99c492..c08cd589 100644 --- a/src/InputParser/deckinput/keywords/props.jl +++ b/src/InputParser/deckinput/keywords/props.jl @@ -179,12 +179,12 @@ function parse_keyword!(data, outer_data, units, cfg, f, ::Val{:PVTG}) end function parse_keyword!(data, outer_data, units, cfg, f, ::Val{:PVTW}) - rec = read_record(f) tdims = [NaN, NaN, NaN, NaN, NaN] utypes = (:pressure, :liquid_formation_volume_factor, :compressibility, :viscosity, :compressibility) nreg = number_of_tables(outer_data, :pvt) out = [] for i = 1:nreg + rec = read_record(f) t = parse_defaulted_line(rec, tdims) swap_unit_system_axes!(t, units, utypes) @assert all(isfinite, t) "PVTW cannot be defaulted, found defaulted record in region $i" @@ -194,12 +194,12 @@ function parse_keyword!(data, outer_data, units, cfg, f, ::Val{:PVTW}) end function parse_keyword!(data, outer_data, units, cfg, f, ::Val{:PVCDO}) - rec = read_record(f) tdims = [NaN, NaN, NaN, NaN, NaN] utypes = (:pressure, :liquid_formation_volume_factor, :compressibility, :viscosity, :compressibility) nreg = number_of_tables(outer_data, :pvt) out = [] for i = 1:nreg + rec = read_record(f) t = parse_defaulted_line(rec, tdims) swap_unit_system_axes!(t, units, utypes) @assert all(isfinite, t) "PVCDO cannot be defaulted, found defaulted record in region $i" diff --git a/src/InputParser/deckinput/keywords/schedule.jl b/src/InputParser/deckinput/keywords/schedule.jl index 68f6eaee..d630135b 100644 --- a/src/InputParser/deckinput/keywords/schedule.jl +++ b/src/InputParser/deckinput/keywords/schedule.jl @@ -81,7 +81,12 @@ end function parse_keyword!(data, outer_data, units, cfg, f, ::Val{:WCONPROD}) d = "Default" - defaults = [d, "OPEN", d, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, 0] + defaults = [ + d, "OPEN", d, -1.0, -1.0, + -1.0, -1.0, -1.0, -1.0, -1.0, + 0, 0.0, NaN, NaN, NaN, + NaN, NaN, NaN, NaN, NaN + ] wells = get_wells(outer_data) wconprod = parse_defaulted_group_well(f, defaults, wells, 1) for kw in wconprod @@ -272,6 +277,10 @@ function parse_keyword!(data, outer_data, units, cfg, f, ::Val{:DATES}) out = Vector{DateTime}() sizehint!(out, length(dt)) for t in dt + if length(t[end]) > 8 + # TODO: Fix dirty hack for skipping ms + t[end] = t[end][1:8] + end push!(out, convert_date_kw(t)) end data["DATES"] = out diff --git a/src/InputParser/deckinput/parser.jl b/src/InputParser/deckinput/parser.jl index dd18bd4f..3f14f2e8 100644 --- a/src/InputParser/deckinput/parser.jl +++ b/src/InputParser/deckinput/parser.jl @@ -64,15 +64,33 @@ end """ parse_data_file(filename; unit = :si) + data = parse_data_file("MY_MODEL.DATA") -Parse a .DATA file (industry standard input file) into a Dict. Units will be -converted to strict SI unless you pass something else like `units = :field`. -Setting `units = nothing` will skip unit conversion. Note that JutulDarcy -assumes that the unit system is internally consistent. It is highly recommended -to parse to the SI units if you want to perform simulations. +Parse a .DATA file given by the `filename` (industry standard input file) into a +Dict. Units will be converted to strict SI unless you pass something else like +`units = :field`. Setting `units = nothing` will skip unit conversion. Note that +JutulDarcy assumes that the unit system is internally consistent. It is highly +recommended to parse to the SI units if you want to perform simulations. -NOTE: This function is experimental and only covers a small portion of the -keywords that exist for various simulators. +The best publicly available documentation on this format is available from the +Open Porous Media (OPM) project's webpages: [OPM Flow manual +](https://opm-project.org/?page_id=955). + +# Keyword arguments +- `warn_parsing=true`: Produce a warning when keywords are not supported (or + partially supported) by the parser. +- `warn_feature`=true`: Produce a warning when keywords are supported, but have + limited or missing support in the numerical solvers. +- `units=:si`: Symbol that indicates the unit system to be used in the output. + Setting this to `nothing` will return values without conversion, i.e. exactly + what is in the input files. `:si` will use strict SI. Other alternatives are + `:field` and `:metric`. `:lab` is currently unsupported. + +# Note +This function is experimental and only covers a small portion of the keywords +that exist for various simulators. You will get warnings that indicate the level +of support for keywords in both the parser and the numerical solvers when known +keywords with limited support. Pull requests for new keywords are welcome! """ function parse_data_file(filename; kwarg...) outer_data = Dict{String, Any}() @@ -201,6 +219,19 @@ function parse_data_file!(outer_data, filename, data = outer_data; return outer_data end +""" + parse_grdecl_file("mygrid.grdecl"; actnum_path = missing, kwarg...) + +Parse a GRDECL file separately from the full input file. Note that the GRID +section does not contain units - passing the `input_units` keyword is therefore +highly recommended. + +# Keyword arguments + - `units=:si`: Units to use for return values. Requires `input_units` to be set. + - `input_units=nothing`: The units the file is given in. + - `verbose=false`: Toggle verbosity. + +""" function parse_grdecl_file(filename; actnum_path = missing, kwarg...) outer_data = Dict{String, Any}() data = new_section(outer_data, :GRID) diff --git a/src/InputParser/deckinput/utils.jl b/src/InputParser/deckinput/utils.jl index b25b80a8..4369d67b 100644 --- a/src/InputParser/deckinput/utils.jl +++ b/src/InputParser/deckinput/utils.jl @@ -15,8 +15,10 @@ function read_record(f; fix = true) end line = strip(line) if !startswith(line, "--") - if endswith(line, '/') - line = strip(rstrip(line, '/')) + if contains(line, '/') + # TODO: Think this is OK for parsing ASCII. + ix = findfirst('/', line) + line = line[1:ix-1] active = false end if length(line) > 0 @@ -44,7 +46,9 @@ function parse_defaulted_group_well(f, defaults, wells, namepos = 1) out = [] line = read_record(f) while length(line) > 0 - parsed = parse_defaulted_line(line, defaults) + allow_wildcard = fill(true, length(defaults)) + allow_wildcard[1] = false + parsed = parse_defaulted_line(line, defaults, allow_wildcard = allow_wildcard) name = parsed[namepos] if occursin('*', name) || occursin('?', name) re = Regex(replace(name, "*" => ".*", "?" => ".")) @@ -78,7 +82,7 @@ function parse_defaulted_line(lines::String, defaults; kwarg...) return parse_defaulted_line([lines], defaults; kwarg...) end -function parse_defaulted_line(lines, defaults; required_num = 0, keyword = "") +function parse_defaulted_line(lines, defaults; required_num = 0, keyword = "", allow_wildcard = missing) out = similar(defaults, 0) sizehint!(out, length(defaults)) pos = 1 @@ -89,25 +93,33 @@ function parse_defaulted_line(lines, defaults; required_num = 0, keyword = "") if length(s) == 0 continue end - if occursin('*', s) && !startswith(s, '\'') # Could be inside a string for wildcard matching + default = defaults[pos] + if ismissing(allow_wildcard) + allow_star = true + else + allow_star = allow_wildcard[pos] + end + if allow_star && occursin('*', s) && !startswith(s, '\'') # Could be inside a string for wildcard matching if s == "*" num_defaulted = 1 else - num_defaulted = Parsers.parse(Int, match(r"\d+\*", s).match[1:end-1]) + parse_wildcard = match(r"\d+\*", s) + if isnothing(parse_wildcard) + error("Unable to parse string for * expansion: $s") + end + num_defaulted = Parsers.parse(Int, parse_wildcard.match[1:end-1]) end for i in 1:num_defaulted push!(out, defaults[pos]) pos += 1 end else - default = defaults[pos] if default isa String converted = strip(s, [' ', '\'']) else T = typeof(default) converted = Parsers.tryparse(T, s) if isnothing(converted) - @info s converted = T.(Parsers.tryparse(Float64, s)) end end @@ -349,36 +361,81 @@ end function parse_keyword!(data, outer_data, units, cfg, f, v::Val{T}) where T # Keywords where we read a single record and don't do anything proper - skip_kw = [ - :PETOPTS, - :PARALLEL, - :VECTABLE, - :MULTSAVE - ] - # Keywords that are a single record where we should warn - skip_kw_with_warn = Symbol[ + skip_kw_with_warn = Symbol[ + :SATOPTS, + :EQLOPTS, + :TRACERS, + :PIMTDIMS, + :FLUXNUM, + :OPTIONS ] - # Single word keywords are trivial to parse, just set a true flag. - single_word_kw = [ - :MULTOUT, - :NOSIM, - :NONNC, - :NEWTRAN - ] - single_word_kw_with_warn = Symbol[] - if T in skip_kw - data["$T"] = read_record(f) - elseif T in skip_kw_with_warn - parser_message(cfg, outer_data, "$T", PARSER_MISSING_SUPPORT) - data["$T"] = read_record(f) - elseif T in single_word_kw - data["$T"] = true - elseif T in single_word_kw_with_warn - parser_message(cfg, outer_data, "$T", PARSER_JUTULDARCY_MISSING_SUPPORT) - data["$T"] = true - else - error("Unhandled keyword $T encountered.") + + skip_list = [] + function skip_kw!(kw, num, msg = nothing) + push!(skip_list, (kw, num, msg)) + end + skip_kw!(:PETOPTS, 1) + skip_kw!(:PARALLEL, 1) + skip_kw!(:MULTSAVE, 1) + skip_kw!(:VECTABLE, 1) + skip_kw!(:MULTSAVE, 1) + + skip_kw!(:MULTOUT, 0) + skip_kw!(:NOSIM, 0) + skip_kw!(:NONNC, 0) + skip_kw!(:NEWTRAN, 0) + + skip_kw!(:SATOPTS, 1, PARSER_MISSING_SUPPORT) + skip_kw!(:EQLOPTS, 1, PARSER_MISSING_SUPPORT) + skip_kw!(:TRACERS, 1, PARSER_MISSING_SUPPORT) + skip_kw!(:PIMTDIMS, 1, PARSER_MISSING_SUPPORT) + skip_kw!(:FLUXNUM, 1, PARSER_MISSING_SUPPORT) + skip_kw!(:OPTIONS, 1, PARSER_MISSING_SUPPORT) + skip_kw!(:EHYSTR, 1, PARSER_MISSING_SUPPORT) + skip_kw!(:SWATINIT, 1, PARSER_MISSING_SUPPORT) + skip_kw!(:ZIPPY2, 1, PARSER_MISSING_SUPPORT) + skip_kw!(:DRSDT, 1, PARSER_MISSING_SUPPORT) + skip_kw!(:WPAVE, 1, PARSER_MISSING_SUPPORT) + skip_kw!(:VAPPARS, 1, PARSER_MISSING_SUPPORT) + skip_kw!(:NETBALAN, 1, PARSER_MISSING_SUPPORT) + + skip_kw!(:TRACER, Inf, PARSER_MISSING_SUPPORT) + skip_kw!(:THPRES, Inf, PARSER_MISSING_SUPPORT) + skip_kw!(:PIMULTAB, Inf, PARSER_MISSING_SUPPORT) + skip_kw!(:VFPPROD, Inf, PARSER_MISSING_SUPPORT) + skip_kw!(:VFPINJ, Inf, PARSER_MISSING_SUPPORT) + skip_kw!(:WTRACER, Inf, PARSER_MISSING_SUPPORT) + skip_kw!(:GCONINJE, Inf, PARSER_MISSING_SUPPORT) + skip_kw!(:WTEST, Inf, PARSER_MISSING_SUPPORT) + + found = false + for (kw, num, msg) in skip_list + if kw != T + continue + end + if !isnothing(msg) + parser_message(cfg, outer_data, "$kw", msg) + end + if num == 0 + # Single word keywords are trivial to parse, just set a true flag. + data["$T"] = true + elseif num == 1 + data["$T"] = read_record(f) + else + skip_record(f) + end + found = true + break + end + + if !found + if startswith("$T", "TVDP") + parser_message(cfg, outer_data, "$T", PARSER_MISSING_SUPPORT) + read_record(f) + else + error("Unhandled keyword $T encountered.") + end end end diff --git a/src/JutulDarcy.jl b/src/JutulDarcy.jl index 1e059813..04db83fd 100644 --- a/src/JutulDarcy.jl +++ b/src/JutulDarcy.jl @@ -1,31 +1,124 @@ + +""" +$(README) + +--- +## Module exports: + +$(EXPORTS) + +""" module JutulDarcy - import Jutul: number_of_cells, number_of_faces, - degrees_of_freedom_per_entity, - values_per_entity, - absolute_increment_limit, relative_increment_limit, maximum_value, minimum_value, - select_primary_variables!, - select_secondary_variables!, - select_equations!, - select_parameters!, - select_minimum_output_variables!, - initialize_primary_variable_ad!, - update_primary_variable!, - update_secondary_variable!, - default_value, - initialize_variable_value!, - initialize_variable_value, - initialize_variable_ad!, - update_half_face_flux!, - update_accumulation!, - update_equation!, - setup_parameters, - count_entities, - count_active_entities, - associated_entity, - active_entities, - number_of_entities, - declare_entities, - get_neighborship + export MultiPhaseSystem, ImmiscibleSystem, SinglePhaseSystem + export reservoir_linsolve + export parse_data_file, parse_grdecl_file + export get_1d_reservoir + export DeckPhaseViscosities + export DeckShrinkageFactors + export CPRPreconditioner + export MuBTable + export ConstMuBTable + export DeckPhaseMassDensities, RelativePermeabilities + export ThreePhaseCompositionalDensitiesLV + export PhaseMassFractions + export PhaseMassFractions + export ThreePhaseLBCViscositiesLV + export plot_well! + export plot_well_results + export plot_reservoir_simulation_result + export plot_reservoir + export simulate_reservoir_parray + export setup_reservoir_simulator_parray + export component_mass_fluxes!, update_total_masses! + export table_to_relperm + export AqueousPhase, LiquidPhase, VaporPhase + export number_of_phases + export SourceTerm + export Pressure, Saturations, TotalMasses, TotalMass + export fluid_volume, pore_volume + export MinimalTPFAGrid + export compute_peaceman_index + export discretized_domain_tpfv_flow + export discretized_domain_well + export StandardBlackOilSystem + export PhaseRelativePermeability + export FlowBoundaryCondition + export ReservoirSimResult + export reservoir_domain + export reservoir_model + export setup_reservoir_model + export setup_reservoir_simulator + export simulate_reservoir + export setup_reservoir_state + export setup_reservoir_forces + export full_well_outputs + export well_output + export well_symbols + export wellgroup_symbols + export available_well_targets + export BlackOilUnknown + export BlackOilX + export TotalSurfaceMassRate + export WellGroup + export DisabledControl + export Wells + export TotalMassVelocityMassFractionsFlow + export BottomHolePressureTarget, TotalRateTarget + export SinglePhaseRateTarget, DisabledTarget + export SurfaceLiquidRateTarget, SurfaceOilRateTarget + export SurfaceWaterRateTarget, SurfaceGasRateTarget + export HistoricalReservoirVoidageTarget, ReservoirVoidageTarget + export SinglePhaseRateTarget, BottomHolePressureTarget + export PerforationMask + export WellDomain, MultiSegmentWell + export TotalMassFlux, PotentialDropBalanceWell + export SegmentWellBoreFrictionHB + export InjectorControl, ProducerControl + export Perforations + export MixedWellSegmentFlow + export segment_pressure_drop + export setup_well, setup_vertical_well + export well_mismatch + export simulate_data_file, setup_case_from_data_file + export get_test_setup, get_well_from_mrst_data + export setup_case_from_mrst + export simulate_mrst_case + export MultiPhaseCompositionalSystemLV + export StandardVolumeSource, VolumeSource, MassSource + export OverallMoleFractions + export ImmiscibleSaturation + export ThermalSystem + export PhaseMassDensities, ConstantCompressibilityDensities + export BrooksCoreyRelativePermeabilities, TabulatedSimpleRelativePermeabilities + + import Jutul: + number_of_cells, number_of_faces, + degrees_of_freedom_per_entity, + values_per_entity, + absolute_increment_limit, relative_increment_limit, maximum_value, minimum_value, + select_primary_variables!, + select_secondary_variables!, + select_equations!, + select_parameters!, + select_minimum_output_variables!, + initialize_primary_variable_ad!, + update_primary_variable!, + update_secondary_variable!, + default_value, + initialize_variable_value!, + initialize_variable_value, + initialize_variable_ad!, + update_half_face_flux!, + update_accumulation!, + update_equation!, + setup_parameters, + count_entities, + count_active_entities, + associated_entity, + active_entities, + number_of_entities, + declare_entities, + get_neighborship import Jutul: update_preconditioner!, partial_update_preconditioner! @@ -50,8 +143,9 @@ module JutulDarcy using PrecompileTools using Dates import DataStructures: OrderedDict + using DocStringExtensions + - export reservoir_linsolve, get_1d_reservoir include("types.jl") include("deck_types.jl") include("porousmedia_grids.jl") @@ -99,7 +193,7 @@ module JutulDarcy include("InputParser/InputParser.jl") using .InputParser import .InputParser: parse_data_file - export parse_data_file + @compile_workload begin precompile_darcy_multimodels() diff --git a/src/blackoil/blackoil.jl b/src/blackoil/blackoil.jl index 64e8a482..64eb787f 100644 --- a/src/blackoil/blackoil.jl +++ b/src/blackoil/blackoil.jl @@ -30,13 +30,13 @@ function select_secondary_variables!(S, system::BlackOilSystem, model) S[:PhaseState] = BlackOilPhaseState() spe1_data = blackoil_bench_pvt(:spe1) pvt = spe1_data[:pvt] - S[:PhaseMassDensities] = DeckDensity(pvt) + S[:PhaseMassDensities] = DeckPhaseMassDensities(pvt) S[:ShrinkageFactors] = DeckShrinkageFactors(pvt) g = physical_representation(model.domain) if !(g isa WellDomain) S[:SurfaceVolumeMobilities] = SurfaceVolumeMobilities() end - S[:PhaseViscosities] = DeckViscosity(pvt) + S[:PhaseViscosities] = DeckPhaseViscosities(pvt) if has_disgas(system) S[:Rs] = Rs() end diff --git a/src/blackoil/variables/density.jl b/src/blackoil/variables/density.jl index ab9c91af..6994aa1b 100644 --- a/src/blackoil/variables/density.jl +++ b/src/blackoil/variables/density.jl @@ -1,5 +1,5 @@ # Density for all three cases -@jutul_secondary function update_deck_density!(rho, m::DeckDensity, model::StandardBlackOilModel, Rs, Rv, ShrinkageFactors, ix) +@jutul_secondary function update_deck_density!(rho, m::DeckPhaseMassDensities, model::StandardBlackOilModel, Rs, Rv, ShrinkageFactors, ix) b = ShrinkageFactors sys = model.system w, o, g = phase_indices(sys) @@ -11,7 +11,7 @@ end end -@jutul_secondary function update_deck_density!(rho, m::DeckDensity, model::VapoilBlackOilModel, Rv, ShrinkageFactors, ix) +@jutul_secondary function update_deck_density!(rho, m::DeckPhaseMassDensities, model::VapoilBlackOilModel, Rv, ShrinkageFactors, ix) b = ShrinkageFactors sys = model.system has_water = has_other_phase(sys) @@ -31,7 +31,7 @@ end end end -@jutul_secondary function update_deck_density!(rho, m::DeckDensity, model::DisgasBlackOilModel, Rs, ShrinkageFactors, ix) +@jutul_secondary function update_deck_density!(rho, m::DeckPhaseMassDensities, model::DisgasBlackOilModel, Rs, ShrinkageFactors, ix) b = ShrinkageFactors sys = model.system has_water = has_other_phase(sys) diff --git a/src/blackoil/variables/variables.jl b/src/blackoil/variables/variables.jl index bad218c3..38c31209 100644 --- a/src/blackoil/variables/variables.jl +++ b/src/blackoil/variables/variables.jl @@ -32,12 +32,23 @@ function Jutul.line_plot_data(model::SimulationModel, ::Rv) return JutulLinePlotData(X[2:end]./1e5, F[2:end], title = "Saturated oil-in-gas ratio", xlabel = "Pressure [bar]", ylabel = "Rv") end +""" + BlackOilUnknown(dr_max = Inf, ds_max = Inf) + +Variable defining the variable switching black-oil variable. The `ds_max` +argument limits the maximum saturation change over a single Newton iteration +when both a oileic and gaseous phase is present. The `dr_max` limits the maximum +change in the undersaturated variable, taken relative to the maximum value of +the undersaturated variable. + +During simulation, this variable can take on the following interpretations: Gas +saturation, Rs or Rv, depending on the phase conditions and miscibility model. +""" Base.@kwdef struct BlackOilUnknown{R} <: ScalarVariable dr_max::R = Inf ds_max::R = 0.2 end -export BlackOilX struct BlackOilX{T} val::T phases_present::PresentPhasesBlackOil diff --git a/src/blackoil/variables/viscosity.jl b/src/blackoil/variables/viscosity.jl index 6173373a..73fb68bb 100644 --- a/src/blackoil/variables/viscosity.jl +++ b/src/blackoil/variables/viscosity.jl @@ -1,4 +1,4 @@ -@jutul_secondary function update_deck_viscosity!(μ, ρ::DeckViscosity, model::StandardBlackOilModel, Pressure, Rs, Rv, ix) +@jutul_secondary function update_deck_viscosity!(μ, ρ::DeckPhaseViscosities, model::StandardBlackOilModel, Pressure, Rs, Rv, ix) pvt, reg = ρ.pvt, ρ.regions w, o, g = phase_indices(model.system) muW = pvt[w] @@ -15,7 +15,7 @@ end end -@jutul_secondary function update_deck_viscosity!(μ, ρ::DeckViscosity, model::VapoilBlackOilModel, Pressure, Rv, ix) +@jutul_secondary function update_deck_viscosity!(μ, ρ::DeckPhaseViscosities, model::VapoilBlackOilModel, Pressure, Rv, ix) pvt, reg = ρ.pvt, ρ.regions w, o, g = phase_indices(model.system) muW = pvt[w] @@ -31,7 +31,7 @@ end end end -@jutul_secondary function update_deck_viscosity!(μ, ρ::DeckViscosity, model::DisgasBlackOilModel, Pressure, Rs, ix) +@jutul_secondary function update_deck_viscosity!(μ, ρ::DeckPhaseViscosities, model::DisgasBlackOilModel, Pressure, Rs, ix) pvt, reg = ρ.pvt, ρ.regions has_wat = has_other_phase(model.system) if has_wat diff --git a/src/cpgrid/cpgrid.jl b/src/cpgrid/cpgrid.jl index 16bbf6f0..53469263 100644 --- a/src/cpgrid/cpgrid.jl +++ b/src/cpgrid/cpgrid.jl @@ -398,6 +398,9 @@ function grid_from_primitives(primitives) nodes, cartdims ) = primitives + if sum(active) == 0 + error("Grid has no active cells.") + end # Faces mapping to nodes faces = Vector{Int}() face_pos = [1] diff --git a/src/cpr.jl b/src/cpr.jl index 3a05a34f..7ca1e470 100644 --- a/src/cpr.jl +++ b/src/cpr.jl @@ -1,4 +1,3 @@ -export CPRPreconditioner struct CPRStorage{P, R, S, F, W, V} A_p::P # pressure system @@ -41,7 +40,12 @@ end """ -Constrained pressure residual + CPRPreconditioner(p = default_psolve(), s = ILUZeroPreconditioner(); strategy = :quasi_impes, weight_scaling = :unit, update_frequency = 1, update_interval = :iteration, partial_update = true) + +Construct a constrained pressure residual (CPR) preconditioner. + +By default, this is a AMG-BILU(0) version (algebraic multigrid for pressure, +block-ILU(0) for the global system). """ mutable struct CPRPreconditioner{P, S} <: JutulPreconditioner storage::Union{CPRStorage, Nothing} @@ -61,14 +65,6 @@ mutable struct CPRPreconditioner{P, S} <: JutulPreconditioner psolver end - -""" - CPRPreconditioner(p = default_psolve(), s = ILUZeroPreconditioner(); strategy = :quasi_impes, weight_scaling = :unit, update_frequency = 1, update_interval = :iteration, partial_update = true) - -Construct a constrained pressure residual (CPR) preconditioner. - -By default, this is a AMG-BILU(0) version (algebraic multigrid for pressure, block-ILU(0) for the global system). -""" function CPRPreconditioner(p = default_psolve(), s = ILUZeroPreconditioner(); strategy = :true_impes, weight_scaling = :unit, diff --git a/src/deck_support.jl b/src/deck_support.jl index a3959cbb..3db88991 100644 --- a/src/deck_support.jl +++ b/src/deck_support.jl @@ -1,6 +1,4 @@ -export DeckViscosity, DeckShrinkage - -@jutul_secondary function update_deck_viscosity!(mu, μ::DeckViscosity, model, Pressure, ix) +@jutul_secondary function update_deck_viscosity!(mu, μ::DeckPhaseViscosities, model, Pressure, ix) pvt, reg = μ.pvt, μ.regions @inbounds for ph in axes(mu, 1) pvt_ph = pvt[ph] @@ -11,7 +9,7 @@ export DeckViscosity, DeckShrinkage end end -@jutul_secondary function update_deck_density!(rho, ρ::DeckDensity, model, Pressure, ix) +@jutul_secondary function update_deck_density!(rho, ρ::DeckPhaseMassDensities, model, Pressure, ix) rhos = reference_densities(model.system) pvt, reg = ρ.pvt, ρ.regions # Note immiscible assumption @@ -60,8 +58,8 @@ end function Jutul.line_plot_data(model::SimulationModel, k::DeckPhaseVariables) deck_pvt_type(::DeckShrinkageFactors) = :shrinkage - deck_pvt_type(::DeckDensity) = :density - deck_pvt_type(::DeckViscosity) = :viscosity + deck_pvt_type(::DeckPhaseMassDensities) = :density + deck_pvt_type(::DeckPhaseViscosities) = :viscosity pvt_type = deck_pvt_type(k) has_reg = !isnothing(k.regions) diff --git a/src/deck_types.jl b/src/deck_types.jl index 925c5902..bbf0383f 100644 --- a/src/deck_types.jl +++ b/src/deck_types.jl @@ -1,33 +1,45 @@ abstract type DeckPhaseVariables <: PhaseVariables end abstract type AbstractReservoirDeckTable end - -export MuBTable, ConstMuBTable - abstract type AbstractTablePVT <: AbstractReservoirDeckTable end +""" + DeckPhaseViscosities(pvt, regions = nothing) -struct DeckViscosity{T, R} <: DeckPhaseVariables +Secondary variable used to evaluate viscosities when a case is generated from a +input file. Typically not instantiated in user scripts. +""" +struct DeckPhaseViscosities{T, R} <: DeckPhaseVariables pvt::T regions::R - function DeckViscosity(pvt; regions = nothing) + function DeckPhaseViscosities(pvt; regions = nothing) check_regions(regions) pvt_t = Tuple(pvt) new{typeof(pvt_t), typeof(regions)}(pvt_t, regions) end end -export DeckDensity, RelativePermeabilities, ThreePhaseCompositionalDensitiesLV, PhaseMassFractions, PhaseMassFractions -export ThreePhaseLBCViscositiesLV -struct DeckDensity{T, R} <: DeckPhaseVariables +""" + DeckPhaseMassDensities(pvt, regions = nothing) + +Secondary variable used to evaluate densities when a case is generated from a +input file. Typically not instantiated in user scripts. +""" +struct DeckPhaseMassDensities{T, R} <: DeckPhaseVariables pvt::T regions::R - function DeckDensity(pvt; regions = nothing) + function DeckPhaseMassDensities(pvt; regions = nothing) check_regions(regions) pvt_t = Tuple(pvt) new{typeof(pvt_t), typeof(regions)}(pvt_t, regions) end end +""" +DeckShrinkageFactors(pvt, regions = nothing) + +Secondary variable used to evaluate shrinkage factors when a case is generated +from a input file. Typically not instantiated in user scripts. +""" struct DeckShrinkageFactors{T, R} <: DeckPhaseVariables pvt::T regions::R @@ -38,7 +50,13 @@ struct DeckShrinkageFactors{T, R} <: DeckPhaseVariables end end +""" + MuBTable(pvt, regions = nothing) +Table used to evaluate viscosities and shrinkage factors when a case is +generated from a input file. Typically used to wrap tables (e.g. PVDG, PVDO) for +use in simulation. +""" struct MuBTable{V, I} pressure::V shrinkage::V @@ -81,6 +99,12 @@ struct ConstMuBTable{R} mu_c::R end +""" + ConstMuBTable(pvtw::M) where M<:AbstractVector + +Create a constant viscosity and formation-volume-factor table from a vector. +Typical usage is to wrap a PVTW type table generated from external software. +""" function ConstMuBTable(pvtw::M) where M<:AbstractVector pvtw = flat_region_expand(pvtw) # Only one region supported atm diff --git a/src/ext.jl b/src/ext.jl index aee0141a..ae121802 100644 --- a/src/ext.jl +++ b/src/ext.jl @@ -1,4 +1,3 @@ -export plot_well!, plot_well_results, plot_reservoir_simulation_result function plot_well! @@ -40,7 +39,6 @@ function plot_reservoir_simulation_result(model::MultiModel, res::ReservoirSimRe return fig end -export plot_reservoir function plot_reservoir(model, arg...; well_fontsize = 18, well_linewidth = 5, kwarg...) rmodel = reservoir_model(model) fig = plot_interactive(rmodel.data_domain, arg...; kwarg...) @@ -55,13 +53,11 @@ function plot_reservoir(model, arg...; well_fontsize = 18, well_linewidth = 5, k return fig end -export simulate_reservoir_parray function simulate_reservoir_parray(case, mode = :mpi; kwarg...) sim, cfg = setup_reservoir_simulator(case; mode = mode, kwarg...) return simulate!(sim, case.dt, forces = case.forces, config = cfg) end -export setup_reservoir_simulator_parray function setup_reservoir_simulator_parray end diff --git a/src/facility/cross_terms.jl b/src/facility/cross_terms.jl index c405089a..d5e33768 100644 --- a/src/facility/cross_terms.jl +++ b/src/facility/cross_terms.jl @@ -1,4 +1,3 @@ -export ReservoirFromWellFlowCT, FacilityFromWellFlowCT, WellFromFacilityFlowCT abstract type AbstractReservoirFromWellCT <: Jutul.AdditiveCrossTerm end struct ReservoirFromWellFlowCT{T<:AbstractVector, I<:AbstractVector} <: AbstractReservoirFromWellCT diff --git a/src/facility/facility.jl b/src/facility/facility.jl index 60abba49..f37c66d8 100644 --- a/src/facility/facility.jl +++ b/src/facility/facility.jl @@ -1,6 +1,3 @@ -export TotalSurfaceMassRate, WellGroup, DisabledControl -export HistoryMode, PredictionMode, Wells - """ (Absolute) Minimum well rate for a well that is not disabled. """ diff --git a/src/facility/types.jl b/src/facility/types.jl index 4360c2f3..e13e65e3 100644 --- a/src/facility/types.jl +++ b/src/facility/types.jl @@ -1,4 +1,3 @@ -export TotalMassVelocityMassFractionsFlow abstract type FacilitySystem <: JutulSystem end struct PredictionMode <: FacilitySystem end @@ -18,6 +17,14 @@ end const WellGroupModel = SimulationModel{WellGroup, <:Any, <:Any, <:Any} struct Wells <: JutulEntity end + +""" + TotalSurfaceMassRate(max_absolute_change = nothing, max_relative_change = nothing) + +Variable, typically representing the primary variable for a [`WellGroup`](@ref). +The variable is a single entry per well and solves for the total surface mass +rate from a well to the facility model. +""" Base.@kwdef struct TotalSurfaceMassRate <: ScalarVariable "Maximum absolute change betweeen two Newton updates (nominally kg/s)" max_absolute_change::Union{Float64, Nothing} = nothing @@ -37,7 +44,10 @@ abstract type WellTarget end abstract type SurfaceVolumeTarget <: WellTarget end """ -Perforations are connections from well cells to reservoir vcells + Perforations() + +Entity that defines perforations: Connections from well cells to reservoir +cells. """ struct Perforations <: JutulEntity end @@ -61,14 +71,26 @@ end Base.show(io::IO, t::SurfaceVolumeTarget) = print(io, "$(typeof(t)) with value $(t.value) [m^3/s] for $(join([typeof(p) for p in lumped_phases(t)], ", "))") -# Basics -export BottomHolePressureTarget, TotalRateTarget, SinglePhaseRateTarget, DisabledTarget -# Phase mixtures -export SurfaceLiquidRateTarget, SurfaceOilRateTarget, SurfaceWaterRateTarget, SurfaceGasRateTarget +""" + BottomHolePressureTarget(q, phase) + +Bottom-hole pressure (bhp) target with target pressure value `bhp`. A well +operating under a bhp constraint will keep the well pressure at the bottom hole +(typically the top of the perforations) fixed at this value unless doing so +would violate other constraints, like the well switching from injection to +production when declared as an injector. -struct BottomHolePressureTarget <: WellTarget - value::AbstractFloat +# Examples +```julia-repl +julia> BottomHolePressureTarget(100e5) +BottomHolePressureTarget with value 100.0 [bar] +``` + +""" +struct BottomHolePressureTarget{T} <: WellTarget + value::T end +Base.show(io::IO, t::BottomHolePressureTarget) = print(io, "BottomHolePressureTarget with value $(convert_from_si(t.value, :bar)) [bar]") """ SinglePhaseRateTarget(q, phase) @@ -82,17 +104,25 @@ SinglePhaseRateTarget of 0.001 [m^3/s] for LiquidPhase() ``` """ -struct SinglePhaseRateTarget <: SurfaceVolumeTarget - value::AbstractFloat - phase::AbstractPhase +struct SinglePhaseRateTarget{T, P} <: SurfaceVolumeTarget + value::T + phase::P end lumped_phases(t::SinglePhaseRateTarget) = (t.phase, ) """ SurfaceLiquidRateTarget(q) -Well target of specified liquid rate with value `q` (liquid/oil and water, but not gas) -at surface conditions. +Well target of specified liquid rate at surface conditions with value `q`. +Typically used for a [`ProducerControl`](@ref) as you have full control over the +mixture composition during injection. + +Liquid rate, sometimes abbreviated LRAT, is made up of the phases that remain +liquid at surface conditions. Typically, this will be water and oil if present +in the model, but never different types of gas. If a producing becomes nearly or +completely flooded by gas the well can go to very high or even infinite flows. +It is therefore important to combine this control with a limit such as a +bottom-hole-pressure constraint. """ struct SurfaceLiquidRateTarget{T} <: SurfaceVolumeTarget where T<:AbstractFloat value::T @@ -104,6 +134,8 @@ lumped_phases(::SurfaceLiquidRateTarget) = (AqueousPhase(), LiquidPhase()) SurfaceOilRateTarget(q) Well target of specified oil rate with value `q` at surface conditions. +Typically used for a [`ProducerControl`](@ref) as oil, for economic reasons, is +rarely injected into the subsurface. Abbreviated as ORAT in some settings. """ struct SurfaceOilRateTarget{T} <: SurfaceVolumeTarget where T<:AbstractFloat value::T @@ -115,6 +147,11 @@ lumped_phases(::SurfaceOilRateTarget) = (LiquidPhase(), ) SurfaceGasRateTarget(q) Well target of specified gas rate with value `q` at surface conditions. + +Often used for both [`InjectorControl`](@ref) [`ProducerControl`](@ref). +Abbreviated as GRAT in some settings. If used for production it is important to +also impose limits, as the well rate may become very high if there is little gas +present. """ struct SurfaceGasRateTarget{T} <: SurfaceVolumeTarget where T<:AbstractFloat value::T @@ -126,6 +163,10 @@ lumped_phases(::SurfaceGasRateTarget) = (VaporPhase(), ) SurfaceWaterRateTarget(q) Well target of specified water rate with value `q` at surface conditions. + +Often used for both [`InjectorControl`](@ref) [`ProducerControl`](@ref). If used +for production it is important to also impose limits, as the well rate may +become very high if there is little water present. """ struct SurfaceWaterRateTarget{T} <: SurfaceVolumeTarget where T<:AbstractFloat value::T @@ -136,7 +177,10 @@ lumped_phases(::SurfaceWaterRateTarget) = (AqueousPhase(), ) """ TotalRateTarget(q) -Well target of specified total rate of all phases with value `q` at surface conditions. +Well target of specified total rate (sum of all phases) with value `q` at surface +conditions. + +Often used for both [`InjectorControl`](@ref) [`ProducerControl`](@ref). """ struct TotalRateTarget{T} <: SurfaceVolumeTarget where T<:AbstractFloat value::T @@ -146,7 +190,10 @@ Base.show(io::IO, t::TotalRateTarget) = print(io, "TotalRateTarget with value $( """ HistoricalReservoirVoidageTarget(q, weights) -Historical RESV target for history matching cases. +Historical RESV target for history matching cases. See +[`ReservoirVoidageTarget`](@ref). For historical rates, the weights described in +that target are computed based on the reservoir pressure and conditions at the +previous time-step. """ struct HistoricalReservoirVoidageTarget{T, K} <: WellTarget where {T<:AbstractFloat, K<:Tuple} value::T @@ -157,7 +204,16 @@ Base.show(io::IO, t::HistoricalReservoirVoidageTarget) = print(io, "HistoricalRe """ ReservoirVoidageTarget(q, weights) -RESV targets with weights for each pseudo-component +RESV target for history matching cases. The `weights` input should +have one entry per phase (or pseudocomponent) in the system. The well control +equation is then: + +``|q_{ctrl} - \\sum_i w_i q_i^s|`` + +where ``q_i^s`` is the surface rate of phase ``i`` and ``w_i`` the weight of +component stream ``i``. + +This constraint is typically set up from .DATA files for black-oil and immiscible cases. """ struct ReservoirVoidageTarget{T, K} <: WellTarget where {T<:AbstractFloat, K<:Tuple} value::T @@ -167,7 +223,8 @@ end """ DisabledTarget(q) -Disabled target used when a well is under `DisabledControl()` only. +Disabled target used when a well is under `DisabledControl()` only. The well +will be disconnected from the surface. """ struct DisabledTarget <: WellTarget end abstract type WellForce <: JutulForce end @@ -176,10 +233,13 @@ abstract type WellControlForce <: WellForce end """ default_limits(ctrl) -Create reasonable default limits for well control `ctrl`, for example to avoid BHP injectors turning into producers. +Create reasonable default limits for well control `ctrl`, for example to avoid +BHP injectors turning into producers. """ +function default_limits(ctrl) + as_limit(ctrl.target) +end -default_limits(ctrl) = as_limit(ctrl.target) as_limit(target) = NamedTuple([Pair(translate_target_to_symbol(target, shortname = true), target.value)]) as_limit(T::DisabledTarget) = nothing as_limit(T::HistoricalReservoirVoidageTarget) = nothing @@ -188,8 +248,10 @@ as_limit(target::ReservoirVoidageTarget) = NamedTuple([Pair(translate_target_to_ """ DisabledControl() -Control that disables a well. If a well is disabled, it is disconnected from the surface network and no flow occurs -between the well and the top side. Mass transfer can still occur inside the well, and between the well and the reservoir. +Control that disables a well. If a well is disabled, it is disconnected from the +surface network and no flow occurs between the well and the top side. Mass +transfer can still occur inside the well, and between the well and the reservoir +unless perforations are also closed by a [`PerforationMask`](@ref). See also [`ProducerControl`](@ref), [`InjectorControl`](@ref). """ @@ -350,7 +412,21 @@ function (D::WellSegmentFlow)(i, ::Cells) return (faces = cd.faces[loc], signs = signs, cells = cells) end -export PerforationMask +""" + mask = PerforationMask(mask::Vector) + +Create a perforation mask. This can be passed to [`setup_forces`](@ref) for a +well under the `mask` argument. The mask should equal the number of perforations +in the well and is applied to the reference well indices in a multiplicative +fashion. For example, if a well named `:Injector` has two perforations, the +following mask would disable the first perforation and decrease the connection +strength for the second perforation by 50%: +```julia +mask = PerforationMask([0.0, 0.5]) +iforces = setup_forces(W, mask = mask) +forces = setup_reservoir_forces(model, control = controls, Injector = iforces) +``` +""" struct PerforationMask{V} <: JutulForce where V<:AbstractVector values::V function PerforationMask(v::T) where T<:AbstractVecOrMat diff --git a/src/facility/wells/mswells_equations.jl b/src/facility/wells/mswells_equations.jl index 776333ad..a90e752b 100644 --- a/src/facility/wells/mswells_equations.jl +++ b/src/facility/wells/mswells_equations.jl @@ -56,7 +56,14 @@ function segment_pressure_drop(f::SegmentWellBoreFrictionHB, v, ρ, μ) return Δp end +""" + PotentialDropBalanceWell(discretization) +Equation for the pressure drop equation in a multisegment well. This equation +lives on the segment between each node and balances the pressure difference +across the segment with the hydrostatic difference and well friction present in +the current flow regime. +""" struct PotentialDropBalanceWell{T} <: JutulEquation flow_discretization::T end diff --git a/src/facility/wells/wells.jl b/src/facility/wells/wells.jl index a25d853e..d642015f 100644 --- a/src/facility/wells/wells.jl +++ b/src/facility/wells/wells.jl @@ -1,11 +1,3 @@ -export WellDomain, MultiSegmentWell -export TotalMassFlux, PotentialDropBalanceWell, SegmentWellBoreFrictionHB - -export InjectorControl, ProducerControl, SinglePhaseRateTarget, BottomHolePressureTarget - -export Perforations -export MixedWellSegmentFlow -export segment_pressure_drop abstract type WellPotentialFlowDiscretization <: PotentialFlowDiscretization end @@ -17,12 +9,29 @@ struct MixedWellSegmentFlow <: WellPotentialFlowDiscretization end include("separator.jl") -# Total velocity in each well segment +""" + TotalMassFlux(scale = si_unit(:day), max_abs = nothing, max_rel = nothing) + +Variable normally used as primary variable. Represents the total mass flux going +through a face. The typical usage is the mass flow through a segment of a +[`MultiSegmentWell`](@ref). + +Note that the flow direction can often switch signs over a segment during a +complex simulation. Setting `max_rel` to something other than `nothing` can +therefore lead to severe convergence issues in the case of flow reversal. + +# Fields (as keyword arguments) + +$FIELDS +""" struct TotalMassFlux{R} <: ScalarVariable + "Scaling for variable" scale::Union{Nothing, R} + "Max absolute change between Newton iterations" max_abs::Union{Nothing, R} + "Maximum relative change between Newton iterations" max_rel::Union{Nothing, R} - function TotalMassFlux(;scale = 3600*24, max_abs = nothing, max_rel = nothing) + function TotalMassFlux(; scale = si_unit(:day), max_abs = nothing, max_rel = nothing) new{Float64}(scale, max_abs, max_rel) end end @@ -52,7 +61,6 @@ function common_well_setup(nr; dz = nothing, WI = nothing, gravity = gravity_con return (WI, gdz) end -export setup_well, setup_vertical_well """ setup_well(D::DataDomain, reservoir_cells; skin = 0.0, Kh = nothing, radius = 0.1, dir = :z) @@ -71,6 +79,7 @@ function setup_well(g, K, reservoir_cells::AbstractVector; skin = 0.0, Kh = nothing, radius = 0.1, + accumulator_volume = missing, simple_well = false, simple_well_regularization = 1.0, WI = missing, @@ -126,7 +135,9 @@ function setup_well(g, K, reservoir_cells::AbstractVector; volumes[i] = h*π*r_i^2 end if simple_well - accumulator_volume = simple_well_regularization*maximum(volumes) + if ismissing(accumulator_volume) + accumulator_volume = simple_well_regularization*maximum(volumes) + end W = SimpleWell(reservoir_cells; WI = WI_computed, volume = accumulator_volume, dz = dz, kwarg...) else # Depth differences are taken care of via centers. diff --git a/src/flux.jl b/src/flux.jl index 6010915b..bf18e551 100644 --- a/src/flux.jl +++ b/src/flux.jl @@ -140,11 +140,14 @@ end end end -export component_mass_fluxes!, update_total_masses! """ - component_mass_fluxes!(q, face, state, model, kgrad, upw) + component_mass_fluxes!(q, face, state, model, flux_type, kgrad, upw) + +Implementation of component fluxes for a given system for a given face. Should +return a `StaticVector` with one entry per component. + +$(SIGNATURES) -Implementation of component fluxes for a given system for a given face. """ function component_mass_fluxes! @@ -155,6 +158,9 @@ end Update total masses for a given system. Number of input arguments varies based on physical system under consideration. + +$(SIGNATURES) + """ function update_total_masses! diff --git a/src/forces/bc.jl b/src/forces/bc.jl index 1b4bef46..2064f34b 100644 --- a/src/forces/bc.jl +++ b/src/forces/bc.jl @@ -1,7 +1,7 @@ """ FlowBoundaryCondition( cell, - pressure = DEFAULT_MINIMUM_PRESSURE, + pressure = DEFAULT_MINIMUM_PRESSURE, temperature = 298.15; fractional_flow = nothing, density = nothing, @@ -9,7 +9,7 @@ trans_thermal = 1e-6 ) -Boundary condition for constant values (pressure/temperature) +Dirchlet boundary condition for constant values (pressure/temperature) at some inflow boundary """ function FlowBoundaryCondition( cell, @@ -36,7 +36,6 @@ end function Jutul.subforce(s::AbstractVector{S}, model) where S<:FlowBoundaryCondition s = deepcopy(s) m = global_map(model.domain) - n = length(s) keep = repeat([false], n) for (i, bc) in enumerate(s) diff --git a/src/forces/sources.jl b/src/forces/sources.jl index 69e88bdd..3373dc54 100644 --- a/src/forces/sources.jl +++ b/src/forces/sources.jl @@ -1,9 +1,19 @@ """ SourceTerm(cell, value; fractional_flow = [1.0], type = MassSource) -Create source term in given `cell` with given total `value`. The optional -`fractional_flow` argument controls how this term is divided over components if -used for inflow. +Create source term in given `cell` with given total `value`. + +The optional `fractional_flow` argument controls how this term is divided over +components if used for inflow and should contain one entry per component in the +system: (`number_of_components(system)`). `fractional_flow` should sum up to +1.0. The `type` argument should be an instance of the `FlowSourceType` enum, +with interpretations as follows: + +- `MassSource`: Source is directly interpreted as component masses. +- `StandardVolumeSource`: Source is volume at standard/surface conditions. + References densities are used to convert into mass sources. +- `VolumeSource`: Source is volume at in-situ / reservoir conditions. + """ function SourceTerm(cell, value; fractional_flow = [1.0], type = MassSource) @assert sum(fractional_flow) == 1.0 "Fractional flow for source term in cell $cell must sum to 1." diff --git a/src/gradients/objectives.jl b/src/gradients/objectives.jl index 88c20e00..0398697e 100644 --- a/src/gradients/objectives.jl +++ b/src/gradients/objectives.jl @@ -1,5 +1,3 @@ -export well_mismatch - function compute_well_qoi(well_model, state, well::Symbol, pos, rhoS, control) well_state = state[well] well_state = convert_to_immutable_storage(well_state) diff --git a/src/input_simulation/data_input.jl b/src/input_simulation/data_input.jl index 92ce3938..3ffd9ea5 100644 --- a/src/input_simulation/data_input.jl +++ b/src/input_simulation/data_input.jl @@ -1,5 +1,3 @@ -export simulate_data_file, setup_case_from_data_file - """ simulate_data_file(inp; parse_arg = NamedTuple(), kwarg...) @@ -16,6 +14,21 @@ function simulate_data_file(data; setup_arg = NamedTuple(), kwarg...) return result end +""" + case = setup_case_from_data_file( + filename; + parse_arg = NamedTuple(), + kwarg... + ) + + case, data = setup_case_from_data_file(filename; include_data = true) + +Set up a [`JutulCase`](@ref) from a standard input file (with extension .DATA). +Additional arguments to the parser can be passed as key/values in a `NamedTuple` +given as `parse_arg`. The optional input `include_data=true` will make the +function return the parsed data in addition the case itself. Additional keyword +arguments are passed on to [`setup_case_from_parsed_data`](@ref). +""" function setup_case_from_data_file( fn; parse_arg = NamedTuple(), @@ -52,11 +65,11 @@ function setup_case_from_parsed_data(datafile; simple_well = true, use_ijk_trans svar = submodel.secondary_variables # PVT pvt = tuple(pvt...) - rho = DeckDensity(pvt) + rho = DeckPhaseMassDensities(pvt) if sys isa StandardBlackOilSystem set_secondary_variables!(submodel, ShrinkageFactors = JutulDarcy.DeckShrinkageFactors(pvt)) end - mu = DeckViscosity(pvt) + mu = DeckPhaseViscosities(pvt) set_secondary_variables!(submodel, PhaseViscosities = mu, PhaseMassDensities = rho) end if k == :Reservoir @@ -90,7 +103,18 @@ function parse_well_from_compdat(domain, wname, cdat, wspecs, compord; simple_we if isnan(rd) rd = nothing end - W = setup_well(domain, wc, name = Symbol(wname), WI = WI, reference_depth = rd, simple_well = simple_well) + if simple_well + avol = 0.1*mean(domain[:volumes][wc])*mean(domain[:porosity][wc]) + else + avol = missing + end + W = setup_well(domain, wc, + name = Symbol(wname), + accumulator_volume = avol, + WI = WI, + reference_depth = rd, + simple_well = simple_well + ) return (W, wc, WI, open) end @@ -404,10 +428,20 @@ function parse_state0_direct_assignment(model, datafile) return init end -function parse_reservoir(data_file) - grid = data_file["GRID"] +function mesh_from_grid_section(f, actnum = missing) + if f isa String + f = parse_grdecl_file(f) + end + f::AbstractDict + if haskey(f, "GRID") + grid = f["GRID"] + else + grid = f + end + if ismissing(actnum) + actnum = get_effective_actnum(grid) + end cartdims = grid["cartDims"] - actnum = get_effective_actnum(grid) if haskey(grid, "COORD") coord = grid["COORD"] zcorn = grid["ZCORN"] @@ -428,6 +462,13 @@ function parse_reservoir(data_file) # We always want to return an unstructured mesh. G = UnstructuredMesh(G) end + return G +end + +function parse_reservoir(data_file) + grid = data_file["GRID"] + cartdims = grid["cartDims"] + G = mesh_from_grid_section(grid) active_ix = G.cell_map nc = number_of_cells(G) # TODO: PERMYY etc for full tensor perm @@ -951,11 +992,13 @@ function producer_control(sys, flag, ctrl, orat, wrat, grat, lrat, bhp; is_hist t = BottomHolePressureTarget(self_val) is_rate = false elseif ctrl == "RESV" - selv_val = -(wrat + orat + grat) + self_val = -(wrat + orat + grat) + w = [wrat, orat, grat] + w = w./sum(w) if is_hist - t = HistoricalReservoirVoidageTarget(selv_val) + t = HistoricalReservoirVoidageTarget(self_val, w) else - t = ReservoirVoidageTarget(selv_val) + t = ReservoirVoidageTarget(self_val, w) end else error("$ctype control not supported") @@ -992,11 +1035,15 @@ function injector_limits(; bhp = Inf, surface_rate = Inf, reservoir_rate = Inf) end function injector_control(sys, streams, name, flag, type, ctype, surf_rate, res_rate, bhp; is_hist = false) + if occursin('*', flag) + # This is a bit of a hack. + flag = "OPEN" + end if flag == "SHUT" || flag == "STOP" ctrl = DisabledControl() lims = nothing else - @assert flag == "OPEN" + @assert flag == "OPEN" "Unsupported well flag: $flag" if ctype == "RATE" is_rate = true t = TotalRateTarget(surf_rate) diff --git a/src/input_simulation/mrst_input.jl b/src/input_simulation/mrst_input.jl index f34fc1fd..0ff1fbfa 100644 --- a/src/input_simulation/mrst_input.jl +++ b/src/input_simulation/mrst_input.jl @@ -1,6 +1,3 @@ -export get_minimal_tpfa_grid_from_mrst, plot_interactive, get_test_setup, get_well_from_mrst_data -export setup_case_from_mrst - function get_mrst_input_path(name) function valid_mat_path(S) base, ext = splitext(S) @@ -320,7 +317,7 @@ end function deck_relperm(props; oil, water, gas, satnum = nothing) if haskey(props, "SCALECRS") - if length(props["SCALECRS"]) == 0 || lowercase(only(props["SCALECRS"])) == "no" + if length(props["SCALECRS"]) == 0 || lowercase(first(props["SCALECRS"])) == "no" @info "Found two-point rel. perm. scaling" scaling = TwoPointKrScale else @@ -352,17 +349,17 @@ function deck_relperm(props; oil, water, gas, satnum = nothing) @assert haskey(props, "SGFN") for (sof3, swfn, sgfn) in zip(props["SOF3"], props["SWFN"], props["SGFN"]) # Water - krw = PhaseRelPerm(swfn[:, 1], swfn[:, 2], label = :w) + krw = PhaseRelativePermeability(swfn[:, 1], swfn[:, 2], label = :w) # Oil pairs so = sof3[:, 1] krow_t = sof3[:, 2] krog_t = sof3[:, 3] - krow = PhaseRelPerm(so, krow_t, label = :ow) - krog = PhaseRelPerm(so, krog_t, label = :og) + krow = PhaseRelativePermeability(so, krow_t, label = :ow) + krog = PhaseRelativePermeability(so, krog_t, label = :og) # Gas - krg = PhaseRelPerm(sgfn[:, 1], sgfn[:, 2], label = :g) + krg = PhaseRelativePermeability(sgfn[:, 1], sgfn[:, 2], label = :g) push!(KRW, krw) push!(KRG, krg) @@ -402,7 +399,7 @@ function deck_relperm(props; oil, water, gas, satnum = nothing) krarg = (g = kr_1, og = kr_2) end end - return ReservoirRelativePermeability(; krarg..., regions = satnum, scaling = scaling) + return ReservoirRelativePermeabilities(; krarg..., regions = satnum, scaling = scaling) end function flat_region_expand(x::AbstractMatrix, n = nothing) @@ -632,11 +629,11 @@ function model_from_mat_deck(G, data_domain, mrst_data, res_context) svar = model.secondary_variables # PVT pvt = tuple(pvt...) - rho = DeckDensity(pvt) + rho = DeckPhaseMassDensities(pvt) if !is_immiscible set_secondary_variables!(model, ShrinkageFactors = DeckShrinkageFactors(pvt)) end - mu = DeckViscosity(pvt) + mu = DeckPhaseViscosities(pvt) set_secondary_variables!(model, PhaseViscosities = mu, PhaseMassDensities = rho) set_deck_specialization!(model, props, satnum, has_oil, has_wat, has_gas) param = setup_parameters(model) @@ -685,10 +682,10 @@ function set_deck_relperm!(vars, param, props; kwarg...) if scaling_type(kr) != NoKrScale ph = kr.phases @assert ph == :wog - param[:RelPermScalingW] = RelPermScalingCoefficients(:w) - param[:RelPermScalingOW] = RelPermScalingCoefficients(:ow) - param[:RelPermScalingOG] = RelPermScalingCoefficients(:og) - param[:RelPermScalingG] = RelPermScalingCoefficients(:g) + param[:RelPermScalingW] = EndPointScalingCoefficients(:w) + param[:RelPermScalingOW] = EndPointScalingCoefficients(:ow) + param[:RelPermScalingOG] = EndPointScalingCoefficients(:og) + param[:RelPermScalingG] = EndPointScalingCoefficients(:g) end end @@ -799,12 +796,17 @@ function init_from_mat(mrst_data, model, param) return init end +""" + setup_case_from_mrst("filename.mat"; kwarg...) + +Set up a [`Jutul.SimulationCase`](@ref) from a MRST-exported .mat file. +""" function setup_case_from_mrst(casename; wells = :ms, backend = :csc, block_backend = true, split_wells = false, use_well_lengths = false, - facility_grouping = :onegroup, + facility_grouping = missing, minbatch = 1000, steps = :full, nthreads = Threads.nthreads(), @@ -820,6 +822,13 @@ function setup_case_from_mrst(casename; wells = :ms, kwarg...) data_domain, mrst_data = reservoir_domain_from_mrst(casename, extraout = true, convert_grid = convert_grid) G = discretized_domain_tpfv_flow(data_domain; kwarg...) + if ismissing(facility_grouping) + if split_wells + facility_grouping = :perwell + else + facility_grouping = :onegroup + end + end # Set up initializers models = OrderedDict() initializer = Dict() @@ -1063,7 +1072,7 @@ function setup_case_from_mrst(casename; wells = :ms, if legacy_output return (models, parameters, initializer, timesteps, forces, mrst_data) else - model = reservoir_multimodel(models, split_wells = split_wells) + model = reservoir_multimodel(models) # Replace various variables - if they are available replace_variables!(model, OverallMoleFractions = OverallMoleFractions(dz_max = dz_max), throw = false) replace_variables!(model, Saturations = Saturations(ds_max = ds_max), throw = false) @@ -1177,26 +1186,36 @@ function mrst_well_ctrl(model, wdata, is_comp, rhoS) return ctrl end -export simulate_mrst_case - """ - simulate_mrst_case(file_name; kwarg...) + simulate_mrst_case(file_name) + simulate_mrst_case(file_name; ) Simulate a MRST case from `file_name` as exported by `writeJutulInput` in MRST. # Arguments - `file_name::String`: The path to a `.mat` file that is to be simulated. -- `extra_outputs::Vector{Symbol} = [:Saturations]`: Additional variables to output from the simulation. + +# Keyword arguments +- `extra_outputs::Vector{Symbol} = [:Saturations]`: Additional variables to + output from the simulation. - `write_output::Bool = true`: Write output (in the default JLD2 format) -- `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 +- `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. +- `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_case_from_mrst`](@ref), +[`setup_reservoir_simulator`](@ref) and [`simulator_config`](@ref) if +applicable. """ function simulate_mrst_case(fn; extra_outputs::Vector{Symbol} = [:Saturations], output_path = nothing, diff --git a/src/io.jl b/src/io.jl index c5da1d6f..30b1c7f2 100644 --- a/src/io.jl +++ b/src/io.jl @@ -1,14 +1,12 @@ -export table_to_relperm - function table_to_relperm(swof; swcon = 0.0, first_label = :w, second_label = :ow) sw = vec(swof[:, 1]) krw = vec(swof[:, 2]) - krw = PhaseRelPerm(sw, krw, label = first_label) + krw = PhaseRelativePermeability(sw, krw, label = first_label) kro = vec(swof[end:-1:1, 3]) so = 1 .- sw so = vec(so[end:-1:1]) @. so = so - swcon - krow = PhaseRelPerm(so, kro, label = second_label) + krow = PhaseRelativePermeability(so, kro, label = second_label) return (krw, krow) end diff --git a/src/multicomponent/multicomponent.jl b/src/multicomponent/multicomponent.jl index d8217252..a69c0cc1 100644 --- a/src/multicomponent/multicomponent.jl +++ b/src/multicomponent/multicomponent.jl @@ -1,7 +1,4 @@ using MultiComponentFlash -export MultiPhaseCompositionalSystemLV -export StandardVolumeSource, VolumeSource, MassSource - const MINIMUM_COMPOSITIONAL_SATURATION = 1e-10 @inline function is_pure_single_phase(s_immiscible) diff --git a/src/multicomponent/variables/others.jl b/src/multicomponent/variables/others.jl index 62d74fb8..5d026d1d 100644 --- a/src/multicomponent/variables/others.jl +++ b/src/multicomponent/variables/others.jl @@ -1,3 +1,8 @@ +""" + PhaseMassFractions(:liquid) + +Variable that defines the component mass fractions in a specific phase. +""" struct PhaseMassFractions{T} <: CompositionalFractions phase::T end diff --git a/src/multicomponent/variables/primary.jl b/src/multicomponent/variables/primary.jl index 53d74701..a99b5863 100644 --- a/src/multicomponent/variables/primary.jl +++ b/src/multicomponent/variables/primary.jl @@ -16,7 +16,6 @@ function update_primary_variable!(state, p::CompositionalFractions, state_symbol Jutul.unit_sum_update!(s, p, model, dx, w) end -export OverallMoleFractions struct OverallMoleFractions <: CompositionalFractions dz_max::Float64 end @@ -68,13 +67,13 @@ function Jutul.increment_norm(dX, state, model, X, pvar::OverallMoleFractions) end """ -A single saturation that represents the "other" phase in a -three phase compositional system where two phases are predicted by an EoS +A single saturation variable that represents the "other" phase in a three phase +compositional system where two phases are predicted by an EoS """ Base.@kwdef struct ImmiscibleSaturation <: ScalarVariable ds_max::Float64 = 0.2 end -maximum_value(::ImmiscibleSaturation) = 1.0# - MINIMUM_COMPOSITIONAL_SATURATION +maximum_value(::ImmiscibleSaturation) = 1.0 minimum_value(::ImmiscibleSaturation) = 0.0 absolute_increment_limit(s::ImmiscibleSaturation) = s.ds_max diff --git a/src/multiphase.jl b/src/multiphase.jl index a6194eff..15d6434d 100644 --- a/src/multiphase.jl +++ b/src/multiphase.jl @@ -1,10 +1,3 @@ -export MultiPhaseSystem, ImmiscibleSystem, SinglePhaseSystem -export AqueousPhase, LiquidPhase, VaporPhase -export number_of_phases, get_short_name, get_name, subscript -export update_linearized_system! -export SourceTerm -export setup_storage, update_equations! -export Pressure, Saturations, TotalMasses, TotalMass # Abstract multiphase system @@ -47,15 +40,28 @@ end function subscript(prefix::String, phase::AbstractPhase) return string(prefix, "_", get_short_name(phase)) end -# Aqueous phase + +""" + AqueousPhase() + +`AbstractPhase` subtype for water-like phases. +""" struct AqueousPhase <: AbstractPhase end get_name(::AqueousPhase) = "Aqueous" -# Liquid phase +""" + LiquidPhase() + +`AbstractPhase` subtype for liquid-like phases. +""" struct LiquidPhase <: AbstractPhase end get_name(::LiquidPhase) = "Liquid" -# Vapor phases +""" + VaporPhase() + +`AbstractPhase` subtype for vapor or gaseous phases. +""" struct VaporPhase <: AbstractPhase end get_name(::VaporPhase) = "Vapor" @@ -84,7 +90,7 @@ struct Pressure <: ScalarVariable end """ - Pressure(; max_abs = nothing, max_rel = 0.2, scale = 1e8, maximum = Inf, minimum = DEFAULT_MINIMUM_PRESSURE) + Pressure(; max_abs = nothing, max_rel = 0.2, scale = si_unit(:bar), maximum = Inf, minimum = DEFAULT_MINIMUM_PRESSURE) Pressure variable definition. `max_abs`/`max_rel` maximum allowable absolute/relative change over a Newton iteration, `scale` is a "typical" value @@ -106,10 +112,10 @@ struct Saturations <: FractionVariables end """ - Saturations(;ds_max = 0.2) + Saturations(; ds_max = 0.2) Saturations as primary variable. `ds_max` controls maximum allowable saturation -change between two Newton iterations. +change between two successive Newton iterations. """ function Saturations(;ds_max = 0.2) Saturations(ds_max) @@ -144,7 +150,11 @@ function Jutul.default_values(model, ::ConnateWater) return swcon end -# Total component masses +""" + TotalMasses() + +Variable that defines total component masses in each cell of the domain. +""" struct TotalMasses <: VectorVariables end function degrees_of_freedom_per_entity(model::SimulationModel{G, S}, v::TotalMasses) where {G<:Any, S<:MultiPhaseSystem} @@ -154,11 +164,26 @@ end struct PhaseMassMobilities <: PhaseVariables end struct PhaseMobilities <: PhaseVariables end -@inline function minimum_value(::TotalMasses) 0 end +@inline function minimum_value(::TotalMasses) + 0.0 +end + +""" + TotalMasses() +Variable that defines total mass of all components in each cell of the domain. +""" struct TotalMass <: ScalarVariable end -@inline function minimum_value(::TotalMass) 0 end +@inline function minimum_value(::TotalMass) + return 0.0 +end +""" + Transmissibilities() + +Variable/parameter used to define the cell-to-cell transmissibilities when using +a two-point flux approximation scheme. +""" struct Transmissibilities <: ScalarVariable end Jutul.variable_scale(::Transmissibilities) = 1e-10 Jutul.minimum_value(::Transmissibilities) = 0.0 @@ -291,7 +316,6 @@ end number_of_equations_per_entity(system::MultiPhaseSystem, e::ConservationLaw) = number_of_components(system) number_of_equations_per_entity(system::SinglePhaseSystem, e::ConservationLaw) = 1 -export fluid_volume, pore_volume function pore_volume(data_domain::DataDomain; throw = true) if haskey(data_domain, :pore_volume, Cells()) pv = data_domain[:pore_volume] diff --git a/src/porousmedia.jl b/src/porousmedia.jl index 4c45c968..b62c4ab1 100644 --- a/src/porousmedia.jl +++ b/src/porousmedia.jl @@ -1,6 +1,5 @@ import Jutul: compute_half_face_trans, compute_face_trans -export compute_peaceman_index function compute_peaceman_index(g::T, K, r, pos; kwarg...) where T<:Jutul.JutulMesh Δ = Jutul.cell_dims(g, pos) K = Jutul.expand_perm(K, dim(g)) @@ -60,42 +59,6 @@ function Jutul.discretize_domain(d::DataDomain, system::MultiPhaseSystem, ::Val{ return discretized_domain_tpfv_flow(d; kwarg...) end -export discretized_domain_tpfv_flow -function discretized_domain_tpfv_flow(geometry; porosity = 0.1, - permeability = 9.869232667160131e-14, # 100 mD - T = nothing, - gdz = nothing, - gravity = true, - pore_volume = nothing, - general_ad = false) - N = geometry.neighbors - error("Deprecated.") - if isnothing(pore_volume) - pore_volume = porosity.*geometry.volumes - end - nc = length(pore_volume) - if isnothing(T) - T = compute_face_trans(geometry, permeability) - end - cc = geometry.cell_centroids - if gravity && size(cc, 1) == 3 - z = vec(cc[3, :]) - else - z = zeros(nc) - end - if isnothing(gdz) - gdz = compute_face_gdz(N, z) - end - G = MinimalTPFAGrid(pore_volume, N, trans = T, gdz = gdz) - if general_ad - flow = PotentialFlow(N, nc) - else - flow = TwoPointPotentialFlowHardCoded(G, ncells = nc) - end - disc = (mass_flow = flow, heat_flow = flow) - return DiscretizedDomain(G, disc) -end - function discretized_domain_tpfv_flow(domain::Jutul.DataDomain; general_ad = false) N = domain[:neighbors] nc = number_of_cells(physical_representation(domain)) @@ -113,7 +76,6 @@ function Jutul.discretize_domain(d::DataDomain{W}, system::MultiPhaseSystem, ::V return discretized_domain_well(physical_representation(d); kwarg...) end -export discretized_domain_well function discretized_domain_well(W::MultiSegmentWell; z = nothing, kwarg...) if isnothing(z) z = vec(W.centers[3, :]) diff --git a/src/porousmedia_grids.jl b/src/porousmedia_grids.jl index 90eb74de..cb8d2fb5 100644 --- a/src/porousmedia_grids.jl +++ b/src/porousmedia_grids.jl @@ -1,7 +1,3 @@ -export MinimalTPFAGrid -export subgrid - -export transfer, get_1d_reservoir import Base.eltype diff --git a/src/test_utils/setup_reservoir.jl b/src/test_utils/setup_reservoir.jl index 932abf44..c0a49ea7 100644 --- a/src/test_utils/setup_reservoir.jl +++ b/src/test_utils/setup_reservoir.jl @@ -75,7 +75,7 @@ function get_test_setup(mesh_or_casename; case_name = "single_phase_simple", con sys = ImmiscibleSystem([L, V]) model = SimulationModel(G, sys, context = context) - kr = BrooksCoreyRelPerm(sys, [2, 3]) + kr = BrooksCoreyRelativePermeabilities(sys, [2, 3]) s = model.secondary_variables s[:RelativePermeabilities] = kr s[:PhaseMassDensities] = ConstantCompressibilityDensities(sys, pRef, rhoLS, cl) @@ -113,7 +113,7 @@ function get_test_setup(mesh_or_casename; case_name = "single_phase_simple", con sys = ImmiscibleSystem([L, V]) model = SimulationModel(G, sys, context = context) - kr = BrooksCoreyRelPerm(sys, [2, 3]) + kr = BrooksCoreyRelativePermeabilities(sys, [2, 3]) s = model.secondary_variables s[:RelativePermeabilities] = kr s[:PhaseMassDensities] = ConstantCompressibilityDensities(sys, pRef, rhoLS, cl) @@ -151,7 +151,7 @@ function get_test_setup(mesh_or_casename; case_name = "single_phase_simple", con sys = ImmiscibleSystem([A, L, V]) model = SimulationModel(G, sys, context = context) - kr = BrooksCoreyRelPerm(sys, [2, 2, 2]) + kr = BrooksCoreyRelativePermeabilities(sys, [2, 2, 2]) s = model.secondary_variables s[:RelativePermeabilities] = kr s[:PhaseMassDensities] = ConstantCompressibilityDensities(sys, pRef, rhoLS, cl) @@ -228,7 +228,7 @@ function get_test_setup(mesh_or_casename; case_name = "single_phase_simple", con sys = MultiPhaseCompositionalSystemLV(eos, (L, V)) model = SimulationModel(G, sys, context = context) - kr = BrooksCoreyRelPerm(sys, [2, 3]) + kr = BrooksCoreyRelativePermeabilities(sys, [2, 3]) s = model.secondary_variables s[:RelativePermeabilities] = kr parameters = setup_parameters(model, Temperature = T0) @@ -269,7 +269,7 @@ function get_test_setup(mesh_or_casename; case_name = "single_phase_simple", con sys = MultiPhaseCompositionalSystemLV(eos, (A, L, V)) model = SimulationModel(G, sys, context = context) - kr = BrooksCoreyRelPerm(sys, [2, 3, 2]) + kr = BrooksCoreyRelativePermeabilities(sys, [2, 3, 2]) s = model.secondary_variables s[:RelativePermeabilities] = kr parameters = setup_parameters(model, Temperature = T0) diff --git a/src/thermal/thermal.jl b/src/thermal/thermal.jl index 20aa6f54..30ee24a4 100644 --- a/src/thermal/thermal.jl +++ b/src/thermal/thermal.jl @@ -1,4 +1,11 @@ -export ThermalSystem + + +""" + ThermalSystem(num_phases = 1, formulation = :Temperature) + +Geothermal system that defines heat transfer through fluid advection and through +the rock itself. Can be combined with a multiphase system using [`Jutul.CompositeSystem`](@ref). +""" struct ThermalSystem{T} <: JutulSystem nph::Int64 function ThermalSystem(; nphases = 1, formulation = :Temperature) diff --git a/src/types.jl b/src/types.jl index 70278123..37a8ff6a 100644 --- a/src/types.jl +++ b/src/types.jl @@ -1,6 +1,13 @@ abstract type AbstractPhase end +""" +Abstract supertype for all multiphase flow systems. +""" abstract type MultiPhaseSystem <: JutulSystem end +""" +Abstract supertype for multicomponent systems, i.e. flow systems where the +number of components is decoupled from the number of phases. +""" abstract type MultiComponentSystem <: MultiPhaseSystem end const DarcyFlowModel = SimulationModel{<:Any, <:MultiPhaseSystem, <:Any, <:Any} @@ -64,7 +71,6 @@ function Base.show(io::IO, sys::MultiPhaseCompositionalSystemLV) print(io, "MultiPhaseCompositionalSystemLV $name with $(MultiComponentFlash.eostype(eos)) EOS with $n components: $cnames") end -export StandardBlackOilSystem struct StandardBlackOilSystem{D, V, W, R, F, T, P, Num} <: BlackOilSystem rs_max::D rv_max::V @@ -166,7 +172,10 @@ end """ ImmiscibleSystem(phases; reference_densities = ones(length(phases))) - ImmiscibleSystem((LiquidPhase(), VaporPhase()), (1000.0, 700.0)) + ImmiscibleSystem((LiquidPhase(), VaporPhase()), reference_densities = (1000.0, 700.0)) + +Immiscible flow system: Each component exists only in a single phase, and the +number of components equal the number of phases. Set up an immiscible system for the given phases with optional reference densitites. This system is easy to specify with [Pressure](@ref) and @@ -199,8 +208,7 @@ end number_of_components(sys::SinglePhaseSystem) = 1 - -struct PhaseRelPerm{T, N} +struct PhaseRelativePermeability{T, N} k::T label::Symbol connate::N @@ -209,9 +217,20 @@ struct PhaseRelPerm{T, N} k_max::N end -export PhaseRelPerm -function PhaseRelPerm(s, k; label = :w, connate = s[1], epsilon = 1e-16) +""" + PhaseRelativePermeability(s, k; label = :w, connate = s[1], epsilon = 1e-16) + +Type that stores a sorted phase relative permeability table (given as vectors of +equal length `s` and `k`): + +``K_r = K(S)`` + +Optionally, a label for the phase, the connate +saturation and a small epsilon value used to avoid extrapolation can be +specified. +""" +function PhaseRelativePermeability(s, k; label = :w, connate = s[1], epsilon = 1e-16) for i in eachindex(s) if i == 1 if s[1] == 0.0 @@ -234,20 +253,24 @@ function PhaseRelPerm(s, k; label = :w, connate = s[1], epsilon = 1e-16) s, k = JutulDarcy.add_missing_endpoints(s, k) JutulDarcy.ensure_endpoints!(s, k, epsilon) kr = get_1d_interpolator(s, k, cap_endpoints = false) - return PhaseRelPerm(kr, label, connate, crit, s_max, k_max) + return PhaseRelativePermeability(kr, label, connate, crit, s_max, k_max) end -(kr::PhaseRelPerm)(S) = kr.k(S) +(kr::PhaseRelativePermeability)(S) = kr.k(S) -function Base.show(io::IO, t::MIME"text/plain", kr::PhaseRelPerm) - println(io, "PhaseRelPerm for $(kr.label):") +function Base.show(io::IO, t::MIME"text/plain", kr::PhaseRelativePermeability) + println(io, "PhaseRelativePermeability for $(kr.label):") println(io, " .k: Internal representation: $(kr.k)") println(io, " Connate saturation = $(kr.connate)") println(io, " Critical saturation = $(kr.critical)") println(io, " Maximum rel. perm = $(kr.k_max) at $(kr.s_max)") end - +""" +MassSource: Source is directly interpreted as component masses. +StandardVolumeSource: Source is volume at standard/surface conditions. References densities are used to convert into mass sources. +VolumeSource: Source is volume at in-situ / reservoir conditions. +""" @enum FlowSourceType begin MassSource StandardVolumeSource @@ -261,7 +284,6 @@ struct SourceTerm{I, F, T} <: JutulForce type::FlowSourceType end -export FlowBoundaryCondition struct FlowBoundaryCondition{I, F, T} <: JutulForce cell::I pressure::F @@ -301,20 +323,27 @@ struct SimpleWell{SC, P, V} <: WellDomain where {SC, P} end """ - SimpleWell(reservoir_cells) + SimpleWell(reservoir_cells; ) Set up a simple well. -NOTE: `setup_vertical_well` or `setup_well` are the recommended way of setting -up wells. +# Note + +[`setup_vertical_well`](@ref) or [`setup_well`](@ref) are the recommended +way of setting up wells. + +# Fields + +$FIELDS + """ function SimpleWell( - reservoir_cells; - name = :Well, - explicit_dp = true, - surface_conditions = default_surface_cond(), - volume = 1000.0, # Regularization volume for well, not a real volume - kwarg... + reservoir_cells; + name = :Well, + explicit_dp = true, + surface_conditions = default_surface_cond(), + volume = 1000.0, # Regularization volume for well, not a real volume + kwarg... ) nr = length(reservoir_cells) WI, gdz = common_well_setup(nr; kwarg...) @@ -322,14 +351,22 @@ function SimpleWell( return SimpleWell(volume, perf, surface_conditions, name, explicit_dp) end struct MultiSegmentWell{V, P, N, A, C, SC, S} <: WellDomain - volumes::V # One per cell - perforations::P # (self -> local cells, reservoir -> reservoir cells, WI -> connection factor) - neighborship::N # Well cell connectivity - top::A # "Top" node where scalar well quantities live - centers::C # Coordinate centers of nodes - surface::SC # p, T at surface - name::Symbol # Symbol that names the well - segment_models::S # Segment pressure drop model for each segment + "One of volumes per node (cell)" + volumes::V + "(self -> local cells, reservoir -> reservoir cells, WI -> connection factor)" + perforations::P + "Well cell connectivity (connections between nodes)" + neighborship::N + "Top node where scalar well quantities live" + top::A + "Coordinate centers of nodes" + centers::C + "pressure and temperature conditions at surface" + surface::SC + "Name of the well as a Symbol" + name::Symbol + "Pressure drop model for seg well segment" + segment_models::S end """ @@ -348,8 +385,15 @@ end Create well perforated in a vector of `reservoir_cells` with corresponding `volumes` and cell `centers`. -NOTE: `setup_vertical_well` or `setup_well` are the recommended way of setting -up wells. +# Note + +[`setup_vertical_well`](@ref) or [`setup_well`](@ref) are the recommended +way of setting up wells. + +# Fields + +$FIELDS + """ function MultiSegmentWell(reservoir_cells, volumes::AbstractVector, centers; N = nothing, @@ -420,15 +464,34 @@ function MultiSegmentWell(reservoir_cells, volumes::AbstractVector, centers; end -export ReservoirSimResult struct ReservoirSimResult + "Well results as a Dict (output from [`full_well_outputs`](@ref))" wells::AbstractDict + "Reservoir states for each time-step" states::AbstractVector + "The time the states and well solutions are given at" time::AbstractVector + "Raw simulation results with more detailed well results and reports of solution progress" result::Jutul.SimResult + "Dict for holding additional useful data connected to the simulation" extra::AbstractDict end +""" + ReservoirSimResult(model, result::Jutul.SimResult, forces, extra = Dict(); kwarg...) + +Create a specific reservoir simulation results that contains well curves, +reservoir states, and so on. This is the return type from `simulate_reservoir`. + +A `ReservoirSimResult` can be unpacked into wells and states: +```julia +ws, states = res_result +``` + +# Fields + +$FIELDS +""" function ReservoirSimResult(model, result::Jutul.SimResult, forces, extra = Dict(); kwarg...) for (k, v) in kwarg extra[k] = v diff --git a/src/utils.jl b/src/utils.jl index ab1d850d..8fc875f8 100644 --- a/src/utils.jl +++ b/src/utils.jl @@ -1,4 +1,3 @@ -export reservoir_model reservoir_model(model) = model reservoir_model(model::MultiModel) = model.models.Reservoir @@ -6,7 +5,6 @@ reservoir_storage(model, storage) = storage reservoir_storage(model::MultiModel, storage) = storage.Reservoir -export reservoir_domain """ reservoir_domain(g; permeability = convert_to_si(0.1, :darcy), porosity = 0.1, kwarg...) @@ -51,7 +49,6 @@ function reservoir_domain(case::JutulCase) return reservoir_domain(case.model) end -export setup_reservoir_model """ setup_reservoir_model(reservoir, system; wells = [], ) setup_reservoir_model(reservoir, system; wells = [], context = DefaultContext(), reservoir_context = nothing, backend = :csc, ) @@ -183,7 +180,6 @@ function set_reservoir_variable_defaults!(model; p_min, p_max, dp_max_abs, dp_ma return model end -export setup_reservoir_simulator """ setup_reservoir_simulator(models, initializer, parameters = nothing; ) @@ -332,8 +328,29 @@ function setup_reservoir_simulator(case::JutulCase; return (sim, cfg) end -export simulate_reservoir +""" + simulate_reservoir(state0, model, dt; + parameters = setup_parameters(model), + restart = false, + forces = setup_forces(model), + kwarg... + ) + simulate_reservoir(case; + kwarg... + ) + +Convenience function for simulating a reservoir model. This function internally +calls [`setup_reservoir_simulator`](@ref), simulates the problem and returns a +[`ReservoirSimResult`](@ref). + +You can optionally unpack this result into the most typical desired outputs: + +`wellsols, states = simulate_reservoir(...)` + +where `wellsols` contains the well results and `states` the reservoir results +(pressure, saturations and so on, in each cell of the reservoir domain). +""" function simulate_reservoir(state0, model, dt; parameters = setup_parameters(model), restart = false, @@ -521,7 +538,6 @@ function reservoir_multimodel(models::AbstractDict; specialize = false, split_we return model end -export setup_reservoir_state """ setup_reservoir_state(model, ) # Ex: For immiscible two-phase @@ -611,7 +627,6 @@ function handle_alternate_primary_variable_spec!(res_init, found, system) return res_init end -export setup_reservoir_forces """ setup_reservoir_forces(model; control = nothing, limits = nothing, set_default_limits = true, ) @@ -651,7 +666,6 @@ function setup_reservoir_forces(model::MultiModel; control = nothing, limits = n end 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) @@ -742,6 +756,11 @@ function well_output(model::MultiModel, states, well_symbol, forces, target = Bo return d end +""" + well_symbols(model::MultiModel) + +Get the keys of a `MultiModel` models that correspond to well models. +""" function well_symbols(model::MultiModel) models = model.models symbols = Vector{Symbol}() diff --git a/src/variables/pvt.jl b/src/variables/pvt.jl index c7edab88..cf45c18f 100644 --- a/src/variables/pvt.jl +++ b/src/variables/pvt.jl @@ -1,11 +1,37 @@ """ -Mass density of each phase +Abstract type representing the evaluation of mass density of each phase (i.e. +units of mass per units of volume, for each cell in the model domain.) """ abstract type PhaseMassDensities <: PhaseVariables end + +""" + ConstantCompressibilityDensities( + sys_or_nph::Union{MultiPhaseSystem, Integer}, + reference_pressure = 1.0, + reference_density = 0.0, + compressibility = 1.0 + ) + +Secondary variable that implements a constant compressibility relationship for +density. Given the reference pressure, compressibility and density at the +reference pressure, each phase density can be computed as: + +``ρ(S) = ρ_{ref} e^{(p - p_{ref})c}`` + +The constructor can take in either one value per phase or a single value for all +phases for the reference pressure, compressibility and density at reference +conditions. + +## Fields +$FIELDS +""" struct ConstantCompressibilityDensities{T} <: PhaseMassDensities + "Reference pressure for each phase (where the reference densities are given)" reference_pressure::T + "Densities at the reference point" reference_densities::T + "Compressibility factor used when expanding around reference pressure, typically between 1e-3 and 1e-10" compressibility::T function ConstantCompressibilityDensities(sys_or_nph::Union{MultiPhaseSystem, Integer}, reference_pressure = 101325.0, reference_density = 1000.0, compressibility = 1e-10) if isa(sys_or_nph, Integer) @@ -96,6 +122,15 @@ end end end +""" + FluidVolume() + +Variable typically taken to be a parameter. Represents the per-cell volume that +where multiphase flow can occur. For a well, this is the volume inside the +well-bore free flow can occur. For a porous medium, this is the void space +inside the pores that is, to some extent, connected and open to flow (effective +pore-volume). +""" struct FluidVolume <: ScalarVariable end Jutul.minimum_value(::FluidVolume) = eps() diff --git a/src/variables/relperm.jl b/src/variables/relperm.jl index 303b1746..0d2478c8 100644 --- a/src/variables/relperm.jl +++ b/src/variables/relperm.jl @@ -8,6 +8,38 @@ function Jutul.default_value(model, v::AbstractRelativePermeabilities) end struct RelativePermeabilitiesParameter <: AbstractRelativePermeabilities end + +""" + RelativePermeabilities((kr1, kr2, ...)) + +A simple relative permeability implementation. Assumes that each phase has a +relative permeability on the form: + +``K_{r,phase} = F(S_{phase})`` + +Supports multiple fluid regions through the `regions` keyword. + +# Examples +Single region: +``` +kr1 = S -> S^2 +kr2 = S -> S^3 + +kr = RelativePermeabilities((kr1, kr2)) +``` +Two regions: +``` +kr1_reg1 = S -> S^2 +kr2_reg1 = S -> S^3 + +kr1_reg2 = S -> S^3 +kr2_reg2 = S -> S^4 + +regions # should be a vector with one entry that is 1 or 2 for each cell in the domain + +kr = RelativePermeabilities(((kr1_reg1, kr2_reg1), (kr1_reg2, kr2_reg2)), regions = regions) +``` +""" struct RelativePermeabilities{K, R} <: AbstractRelativePermeabilities relperms::K regions::R @@ -65,12 +97,33 @@ function Jutul.line_plot_data(model::SimulationModel, k::RelativePermeabilities) return data end -struct BrooksCoreyRelPerm{V, T} <: AbstractRelativePermeabilities +""" + BrooksCoreyRelativePermeabilities( + sys_or_nph::Union{MultiPhaseSystem, Integer}, + exponents = 1.0, + residuals = 0.0, + endpoints = 1.0 + ) + +Secondary variable that implements the family of Brooks-Corey relative +permeability functions. This is a simple analytical expression for relative +permeabilities that has a limited number of parameters: + +``K(S) = K_{max} * ((S - S_r)/(1 - S_r^{tot}))^N`` + +## Fields +$FIELDS +""" +struct BrooksCoreyRelativePermeabilities{V, T} <: AbstractRelativePermeabilities + "Exponents for each phase" exponents::V + "Residual saturations for each phase" residuals::V + "Maximum relative permeability for each phase" endpoints::V + "Total residual saturation over all phases" residual_total::T - function BrooksCoreyRelPerm(sys_or_nph::Union{MultiPhaseSystem, Integer}, exponents = 1.0, residuals = 0.0, endpoints = 1.0) + function BrooksCoreyRelativePermeabilities(sys_or_nph::Union{MultiPhaseSystem, Integer}, exponents = 1.0, residuals = 0.0, endpoints = 1.0) if isa(sys_or_nph, Integer) nph = sys_or_nph else @@ -86,13 +139,16 @@ struct BrooksCoreyRelPerm{V, T} <: AbstractRelativePermeabilities end """ -Interpolated multiphase rel. perm. that is simple (single region, no magic for more than two phases) + TabulatedSimpleRelativePermeabilities(s::AbstractVector, kr::AbstractVector; regions::Union{AbstractVector, Nothing} = nothing, kwarg...) + +Interpolated multiphase relative permeabilities that assumes that the relative +permeability of each phase depends only on the phase saturation of that phase. """ -struct TabulatedRelPermSimple{V, M, I} <: AbstractRelativePermeabilities +struct TabulatedSimpleRelativePermeabilities{V, M, I} <: AbstractRelativePermeabilities s::V kr::M interpolators::I - function TabulatedRelPermSimple(s::AbstractVector, kr::AbstractVector; regions::Union{AbstractVector, Nothing} = nothing, kwarg...) + function TabulatedSimpleRelativePermeabilities(s::AbstractVector, kr::AbstractVector; regions::Union{AbstractVector, Nothing} = nothing, kwarg...) nph = length(kr) n = length(kr[1]) @assert nph > 0 @@ -116,28 +172,31 @@ struct TabulatedRelPermSimple{V, M, I} <: AbstractRelativePermeabilities end end -""" -Interpolated multiphase rel. perm. that is simple (single region, no magic for more than two phases) -""" -struct ReservoirRelativePermeability{Scaling, ph, O, OW, OG, G, R} <: AbstractRelativePermeabilities +struct ReservoirRelativePermeabilities{Scaling, ph, O, OW, OG, G, R} <: AbstractRelativePermeabilities + "Water relative permeability as a function of water saturation: ``k_{rw}(S_w)``" krw::O + "Oil relative permeability (in the presence of water) as a function of oil saturation: ``k_{row}(S_o)``" krow::OW + "Gas relative permeability (in the presence of water) as a function of gas saturation: ``k_{row}(S_g)``" krog::OG + "Gas relative permeability as a function of gas saturation: ``k_{rg}(S_g)``" krg::G + "Regions to use for each cell of the domain. Can be `nothing` if a single region is used throughout the domain." regions::R + "Symbol designating the type of system, :wog for three-phase, :og for oil-gas, :wg for water-gas, etc." phases::Symbol end -struct RelPermScalingCoefficients{phases} <: VectorVariables +struct EndPointScalingCoefficients{phases} <: VectorVariables end -function RelPermScalingCoefficients(phases::Symbol) - return RelPermScalingCoefficients{phases}() +function EndPointScalingCoefficients(phases::Symbol) + return EndPointScalingCoefficients{phases}() end -Jutul.degrees_of_freedom_per_entity(model, ::RelPermScalingCoefficients) = 4 +Jutul.degrees_of_freedom_per_entity(model, ::EndPointScalingCoefficients) = 4 -function Jutul.default_values(model, scalers::RelPermScalingCoefficients{P}) where P +function Jutul.default_values(model, scalers::EndPointScalingCoefficients{P}) where P nc = number_of_cells(model.domain) relperm = Jutul.get_variable(model, :RelativePermeabilities) kr = relperm[P] @@ -155,7 +214,20 @@ function Jutul.default_values(model, scalers::RelPermScalingCoefficients{P}) whe return kscale end -function ReservoirRelativePermeability(; w = nothing, g = nothing, ow = nothing, og = nothing, scaling = NoKrScale, regions = nothing) +""" + ReservoirRelativePermeabilities( + w = nothing, g = nothing, ow = nothing, og = nothing, + scaling = NoKrScale, regions = nothing) + +Relative permeability with advanced features for reservoir simulation. Includes +features like rel. perm. endpoint scaling, connate water adjustment and separate +phase pair relative permeabilites for the oil phase. + +# Fields + +$FIELDS +""" +function ReservoirRelativePermeabilities(; w = nothing, g = nothing, ow = nothing, og = nothing, scaling = NoKrScale, regions = nothing) has_w = !isnothing(w) has_g = !isnothing(g) has_og = !isnothing(og) @@ -179,7 +251,7 @@ function ReservoirRelativePermeability(; w = nothing, g = nothing, ow = nothing, phases = :og end else - error("ReservoirRelativePermeability only implements two-phase (WO, OG, WG) or three-phase (WOG)") + error("ReservoirRelativePermeabilities only implements two-phase (WO, OG, WG) or three-phase (WOG)") end F = x -> region_wrap(x, regions) @@ -188,10 +260,10 @@ function ReservoirRelativePermeability(; w = nothing, g = nothing, ow = nothing, krog = F(og) krg = F(g) - return ReservoirRelativePermeability{scaling, phases, typeof(krw), typeof(krow), typeof(krog), typeof(krg), typeof(regions)}(krw, krow, krog, krg, regions, phases) + return ReservoirRelativePermeabilities{scaling, phases, typeof(krw), typeof(krow), typeof(krog), typeof(krg), typeof(regions)}(krw, krow, krog, krg, regions, phases) end -function Base.getindex(m::ReservoirRelativePermeability, s::Symbol) +function Base.getindex(m::ReservoirRelativePermeabilities, s::Symbol) if s == :w return m.krw elseif s == :g @@ -205,9 +277,9 @@ function Base.getindex(m::ReservoirRelativePermeability, s::Symbol) end end -scaling_type(::ReservoirRelativePermeability{T}) where T = T +scaling_type(::ReservoirRelativePermeabilities{T}) where T = T -function Jutul.line_plot_data(model::SimulationModel, k::ReservoirRelativePermeability) +function Jutul.line_plot_data(model::SimulationModel, k::ReservoirRelativePermeabilities) s = collect(0:0.01:1) has_reg = !isnothing(k.regions) if has_reg @@ -251,15 +323,15 @@ function Jutul.line_plot_data(model::SimulationModel, k::ReservoirRelativePermea return data end -function Jutul.subvariable(k::ReservoirRelativePermeability, map::FiniteVolumeGlobalMap) +function Jutul.subvariable(k::ReservoirRelativePermeabilities, map::FiniteVolumeGlobalMap) c = map.cells regions = Jutul.partition_variable_slice(k.regions, c) scaling = scaling_type(k) - return ReservoirRelativePermeability(; w = k.krw, ow = k.krow, og = k.krog, g = k.krg, regions = regions, scaling = scaling) + return ReservoirRelativePermeabilities(; w = k.krw, ow = k.krow, og = k.krog, g = k.krg, regions = regions, scaling = scaling) end -function Base.show(io::IO, t::MIME"text/plain", kr::ReservoirRelativePermeability) - println(io, "ReservoirRelativePermeability") +function Base.show(io::IO, t::MIME"text/plain", kr::ReservoirRelativePermeabilities) + println(io, "ReservoirRelativePermeabilities") println(io, " functions:") for f in [:krw, :krow, :krog, :krg] k = getfield(kr, f) @@ -280,7 +352,7 @@ function Base.show(io::IO, t::MIME"text/plain", kr::ReservoirRelativePermeabilit end -@jutul_secondary function update_kr!(kr, kr_def::BrooksCoreyRelPerm, model, Saturations, ix) +@jutul_secondary function update_kr!(kr, kr_def::BrooksCoreyRelativePermeabilities, model, Saturations, ix) n, sr, kwm, sr_tot = kr_def.exponents, kr_def.residuals, kr_def.endpoints, kr_def.residual_total for i in ix for ph in axes(kr, 1) @@ -297,7 +369,7 @@ function brooks_corey_relperm(s::T, n::Real, sr::Real, kwm::Real, sr_tot::Real) return kwm*sat^n end -@jutul_secondary function update_kr!(kr, kr_def::TabulatedRelPermSimple, model, Saturations, ix) +@jutul_secondary function update_kr!(kr, kr_def::TabulatedSimpleRelativePermeabilities, model, Saturations, ix) I = kr_def.interpolators for c in ix @inbounds for j in eachindex(I) @@ -308,7 +380,7 @@ end return kr end -@jutul_secondary function update_kr!(kr, relperm::ReservoirRelativePermeability{NoKrScale, :wog}, model, Saturations, ConnateWater, ix) +@jutul_secondary function update_kr!(kr, relperm::ReservoirRelativePermeabilities{NoKrScale, :wog}, model, Saturations, ConnateWater, ix) s = Saturations phases = phase_indices(model.system) for c in ix @@ -317,7 +389,7 @@ end return kr end -@jutul_secondary function update_kr!(kr, relperm::ReservoirRelativePermeability{NoKrScale, :wo}, model, Saturations, ix) +@jutul_secondary function update_kr!(kr, relperm::ReservoirRelativePermeabilities{NoKrScale, :wo}, model, Saturations, ix) s = Saturations regions = relperm.regions indices = phase_indices(model.system) @@ -327,7 +399,7 @@ end return kr end -@jutul_secondary function update_kr!(kr, relperm::ReservoirRelativePermeability{NoKrScale, :og}, model, Saturations, ix) +@jutul_secondary function update_kr!(kr, relperm::ReservoirRelativePermeabilities{NoKrScale, :og}, model, Saturations, ix) s = Saturations regions = relperm.regions indices = phase_indices(model.system) @@ -337,7 +409,7 @@ end return kr end -@jutul_secondary function update_kr!(kr, relperm::ReservoirRelativePermeability{NoKrScale, :wg}, model, Saturations, ix) +@jutul_secondary function update_kr!(kr, relperm::ReservoirRelativePermeabilities{NoKrScale, :wg}, model, Saturations, ix) s = Saturations regions = relperm.regions indices = phase_indices(model.system) @@ -366,7 +438,7 @@ Base.@propagate_inbounds function two_phase_relperm!(kr, s, regions, Kr_1, Kr_2, kr[i2, c] = kr2(sg) end -@jutul_secondary function update_kr_with_scaling!(kr, relperm::ReservoirRelativePermeability{<:Any, :wog}, model, Saturations, RelPermScalingW, RelPermScalingOW, RelPermScalingOG, RelPermScalingG, ConnateWater, ix) +@jutul_secondary function update_kr_with_scaling!(kr, relperm::ReservoirRelativePermeabilities{<:Any, :wog}, model, Saturations, RelPermScalingW, RelPermScalingOW, RelPermScalingOG, RelPermScalingG, ConnateWater, ix) s = Saturations phases = phase_indices(model.system) scalers = (w = RelPermScalingW, ow = RelPermScalingOW, og = RelPermScalingOG, g = RelPermScalingG) @@ -396,7 +468,7 @@ Base.@propagate_inbounds @inline function update_three_phase_relperm!(kr, relper kr[g, c] = Krg end -function get_kr_scalers(kr::PhaseRelPerm) +function get_kr_scalers(kr::PhaseRelativePermeability) return (kr.connate, kr.critical, kr.s_max, kr.k_max) end function get_kr_scalers(scaler::AbstractMatrix, c) diff --git a/src/variables/variables.jl b/src/variables/variables.jl index 0531bcd0..0ecc1a1c 100644 --- a/src/variables/variables.jl +++ b/src/variables/variables.jl @@ -1,6 +1,3 @@ -export PhaseMassDensities, ConstantCompressibilityDensities -export BrooksCoreyRelPerm, TabulatedRelPermSimple - include("capillary.jl") include("relperm.jl") include("pvt.jl") @@ -31,7 +28,7 @@ function select_default_darcy_secondary_variables!(S, domain, system, formulatio is_multiphase = !isa(system, SinglePhaseSystem) is_bo = system isa BlackOilSystem if is_multiphase && !is_well - S[:RelativePermeabilities] = BrooksCoreyRelPerm(system) + S[:RelativePermeabilities] = BrooksCoreyRelativePermeabilities(system) end if !is_well S[:PhaseMobilities] = PhaseMobilities() @@ -44,7 +41,7 @@ end function select_default_darcy_parameters!(prm, domain, system::SinglePhaseSystem, formulation) prm[:PhaseViscosities] = PhaseViscosities() prm[:FluidVolume] = FluidVolume() - prm[:RelativePermeabilities] = BrooksCoreyRelPerm(system) + prm[:RelativePermeabilities] = BrooksCoreyRelativePermeabilities(system) prm[:Saturations] = Saturations() end diff --git a/test/parray.jl b/test/parray.jl index 6e11881d..5e33f71b 100644 --- a/test/parray.jl +++ b/test/parray.jl @@ -61,7 +61,7 @@ function setup_bl_case(nc, backend = :csr; nstep = nc) ctx = DefaultContext(matrix_layout = BlockMajorLayout()) end model = SimulationModel(domain, sys, context = ctx) - kr = BrooksCoreyRelPerm(sys, [2.0, 2.0], [0.2, 0.2]) + kr = BrooksCoreyRelativePermeabilities(sys, [2.0, 2.0], [0.2, 0.2]) replace_variables!(model, RelativePermeabilities = kr) tot_time = sum(timesteps) pv = pore_volume(domain) diff --git a/test/relperm.jl b/test/relperm.jl index f86502b1..7184a1eb 100644 --- a/test/relperm.jl +++ b/test/relperm.jl @@ -33,7 +33,7 @@ function kr_test_sat(n) end function test_brooks_corey_kr() - bc = BrooksCoreyRelPerm(2, [2.0, 3.0], [0.2, 0.3], [0.9, 1.0]) + bc = BrooksCoreyRelativePermeabilities(2, [2.0, 3.0], [0.2, 0.3], [0.9, 1.0]) S = kr_test_sat(2) kr = similar(S) @test JutulDarcy.update_kr!(kr, bc, nothing, S, entity_eachindex(kr)) ≈ [0.9 0.324 0.0; 0.0 0.064 1.0] @@ -65,7 +65,7 @@ end function test_rel_perm_wrapper() s = [0.1, 0.15, 0.2, 0.8, 1.0] kr = [0.0, 0.0, 0.4, 0.9, 0.9] - relperm = PhaseRelPerm(s, kr) + relperm = PhaseRelativePermeability(s, kr) @testset "Detection of points" begin @test relperm.connate == 0.1 @test relperm.k_max == 0.9 diff --git a/test/sens_bl.jl b/test/sens_bl.jl index 2e5d9eb0..6c58ae75 100644 --- a/test/sens_bl.jl +++ b/test/sens_bl.jl @@ -39,7 +39,7 @@ function setup_bl(;nc = 100, time = 1.0, nstep = 100) # Define system and realize on grid sys = ImmiscibleSystem((LiquidPhase(), VaporPhase())) model = SimulationModel(G, sys) - kr = BrooksCoreyRelPerm(sys, [2.0, 2.0], [0.2, 0.2]) + kr = BrooksCoreyRelativePermeabilities(sys, [2.0, 2.0], [0.2, 0.2]) replace_variables!(model, RelativePermeabilities = kr) tot_time = sum(timesteps) pv = pore_volume(G)