WENO support, improved CPR support for non-standard discretizations, more tests

Project.toml

@@ -62,7 +62,7 @@ ForwardDiff = "0.10.30"
 GLMakie = "0.10.13"
 GeoEnergyIO = "1.1.12"
 HYPRE = "1.6.0, 1.7"
-Jutul = "0.3"
+Jutul = "0.3.2"
 Krylov = "0.9.1"
 LazyArtifacts = "1"
 LinearAlgebra = "1"

docs/make.jl

@@ -30,13 +30,15 @@ function build_jutul_darcy_docs(build_format = nothing;
         "Intro: Simulating Eclipse/DATA input files" => "data_input_file",
         "Intro: Sensitivities in JutulDarcy" => "intro_sensitivities",
         "Intro: Compositional flow" => "co2_brine_2d_vertical",
+        "CO2 injection in saline aquifer" => "co2_sloped",
+        "Consistent discretizations: Average MPFA and nonlinear TPFA" => "consistent_avgmpfa",
+        "High resolution and consistency: WENO, NTPFA and AvgMPFA" => "mpfa_weno_discretizations",
         "Quarter-five-spot with variation" => "five_spot_ensemble",
         "Gravity circulation with CPR preconditioner" => "two_phase_unstable_gravity",
-        "CO2 property calculations" => "co2_props",
-        "Relative permeability" => "relperms",
+        "Properties: CO2-brine correlations with salinity" => "co2_props",
+        "Properties: Relative Permeabilities in JutulDarcy" => "relperms",
         "Model coarsening" => "model_coarsening",
         "History matching a coarse model - CGNet" => "cgnet_egg",
-        "CO2 injection in saline aquifer" => "co2_sloped",
         "Compositional with five components" => "compositional_5components",
         "Parameter matching of Buckley-Leverett" => "optimize_simple_bl",
         "Hybrid simulation with relative permeability" => "hybrid_simulation_relperm",

docs/src/refs.bib

@@ -136,3 +136,35 @@ @article{cgnet1
   year={2024},
   publisher={Springer}
 }
+
+@article{zhang_nonlinear,
+  title={Cell-centered nonlinear finite-volume methods with improved robustness},
+  author={Zhang, Wenjuan and Al Kobaisi, Mohammed},
+  journal={SPE Journal},
+  volume={25},
+  number={01},
+  pages={288--309},
+  year={2020},
+  publisher={SPE}
+}
+
+@article{schneider_nonlinear,
+  title={Monotone nonlinear finite-volume method for challenging grids},
+  author={Schneider, Martin and Flemisch, Bernd and Helmig, Rainer and Terekhov, K and Tchelepi, H},
+  journal={Computational Geosciences},
+  volume={22},
+  pages={565--586},
+  year={2018},
+  publisher={Springer}
+}
+
+@article{raynaud_discretization,
+  title={Toward Accurate Reservoir Simulations on Unstructured Grids: Design of Simple Error Estimators and Critical Benchmarking of Consistent Discretization Methods for Practical Implementation},
+  author={Raynaud, X and Pizzolato, A and Johansson, A and Caresani, F and Ferrari, A and M{\o}yner, O and Nilsen, HM and Cominelli, A and Lie, K-A},
+  journal={SPE Journal},
+  volume={26},
+  number={06},
+  pages={4109--4127},
+  year={2021},
+  publisher={OnePetro}
+}

examples/co2_props.jl

@@ -1,4 +1,4 @@
-# # Plotting CO2 properties used in compositional model
+# # CO2-brine correlations with salinity
 # JutulDarcy includes a set of pre-generated tables for simulation of CO2
 # storage in saline aquifers, as well as functions for calculating PVT and
 # solubility properties of CO2-brine mixtures with varying degrees of salinity.

examples/consistent_avgmpfa.jl

@@ -0,0 +1,129 @@
+# # Consistent discretizations: Average MPFA and nonlinear TPFA
+# This example demonstrates how to use alternative discretizations for the
+# pressure gradient term in the Darcy equation, i.e. the approximation of the
+# Darcy flux:
+#
+# ``\mathbf{}{K}(\nabla p + \rho g \Delta z)``.
+#
+# It is well-known that for certain combinations of grid geometry and
+# permeability fields, the classical two-point flux approximation scheme can
+# give incorrect results. This is due to the fact that the TPFA scheme is not a
+# formally consistent method when the product of the permeability tensor and the
+# normal vector does not align with the cell-to-cell vectors over a face (lack
+# of K-orthogonality).
+#
+# In such cases, it is often beneficial to use a consistent discretization.
+# JutulDarcy includes a class of linear and nonlinear schemes that are designed
+# to be accurate even for challenging grids.
+# 
+# For further details on this class of methods, which differ a bit from the
+# classical MPFA-O type method often seen in the literature, see
+# [schneider_nonlinear](@cite), [zhang_nonlinear](@cite) and
+# [raynaud_discretization](@cite).
+
+# ## Define a mesh and twist the nodes
+# This makes the mesh non K-orthogonal and will lead to wrong solutions for the
+# default TPFA scheme.
+using Jutul
+using JutulDarcy
+using LinearAlgebra
+using GLMakie
+using Test #src
+
+sys = SinglePhaseSystem()
+nx = nz = 100
+pdims = (1.0, 1.0)
+g = CartesianMesh((nx, nz), pdims)
+
+g = UnstructuredMesh(g)
+D = dim(g)
+
+v = 0.1
+for i in eachindex(g.node_points)
+    x, y = g.node_points[i]
+    shiftx =  v*sin(π*x)*sin(3*(-π/2 + π*y))
+    shifty =  v*sin(π*y)*sin(3*(-π/2 + π*x))
+    g.node_points[i] += [shiftx, shifty]
+end
+
+nc = number_of_cells(g)
+domain = reservoir_domain(g, permeability = 0.1*si_unit(:darcy))
+
+fig = Figure()
+Jutul.plot_mesh_edges!(Axis(fig[1, 1]), g)
+fig
+# ## Create a test problem function
+# We set up a problem for our given domain with left and right boundary boundary
+# conditions that correspond to a linear pressure drop. We can expect the
+# steady-state pressure solution to be linear between the two faces as there is
+# no variation in permeability or significant compressibility. The function will
+# return the pressure solution at the end of the simulation for a given scheme.
+function solve_test_problem(scheme)
+    model, parameters = setup_reservoir_model(domain, sys,
+        general_ad = true,
+        kgrad = scheme,
+        block_backend = false
+    )
+    state0 = setup_reservoir_state(model, Pressure = 1e5)
+    nc = number_of_cells(g)
+    bcells = Int64[]
+    bpres = Float64[]
+    for k in 1:nz
+        bnd_l = Jutul.cell_index(g, (1, k))
+        bnd_r = Jutul.cell_index(g, (nx, k))
+        push!(bcells, bnd_l)
+        push!(bpres, 1e5)
+        push!(bcells, bnd_r)
+        push!(bpres, 2e5)
+    end
+    bc = flow_boundary_condition(bcells, domain, bpres)
+    forces = setup_reservoir_forces(model, bc = bc)
+
+    dt = [si_unit(:day)]
+    _, states = simulate_reservoir(state0, model, dt,
+        forces = forces, failure_cuts_timestep = false,
+        tol_cnv = 1e-6,
+        linear_solver = GenericKrylov(preconditioner = AMGPreconditioner(:smoothed_aggregation), rtol = 1e-6)
+        )
+    return states[end][:Pressure]
+end
+# ## Solve the test problem with three different schemes
+# - TPFA (two-point flux approximation, inconsistent, linear)
+# - Average MPFA (consistent, linear)
+# - NTPFA (nonlinear two-point flux approximation, consistent, nonlinear)
+results = Dict()
+for m in [:tpfa, :avgmpfa, :ntpfa]
+    println("Solving $m")
+    results[m] = solve_test_problem(m);
+end
+# ## Plot the results
+# We plot the pressure solution for each of the schemes, as well as the error.
+# Note that the color axis varies between error plots. As the grid is quite
+# skewed, we observe significant errors for the TPFA scheme, with no significant
+# error for the consistent schemes.
+x = domain[:cell_centroids][1, :]
+get_ref(x) = x*1e5 + 1e5
+x_distinct = sort(unique(x))
+sol = get_ref.(x_distinct)
+
+fig = Figure(size = (1200, 600))
+for (i, (m, s)) in enumerate(results)
+    ref = get_ref.(x)
+    err = norm(s .- ref)/norm(ref)
+
+    ax = Axis(fig[1, i], title = "$m")
+    lines!(ax, x_distinct, sol, color = :red)
+    scatter!(ax, x, s, label = m, markersize = 1, alpha = 0.5, transparency = true)
+
+    ax2 = Axis(fig[2, i], title = "Error=$err")
+    Δ = s - ref
+    emin, emax = extrema(Δ)
+    largest = max(abs(emin), abs(emax))
+    crange = (-largest, largest)
+    plt = plot_cell_data!(ax2, g, Δ, colorrange = crange, colormap = :seismic)
+    Colorbar(fig[3, i], plt, vertical = false)
+    if m != :tpfa  #src
+        @test err < 1e-4 #src
+    end  #src
+end
+fig