diff --git a/Project.toml b/Project.toml index 5f5bc461..dc6e11a0 100644 --- a/Project.toml +++ b/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" diff --git a/docs/make.jl b/docs/make.jl index 49f3e800..b1c44e95 100644 --- a/docs/make.jl +++ b/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", diff --git a/docs/src/refs.bib b/docs/src/refs.bib index 561f5f01..630ee16f 100644 --- a/docs/src/refs.bib +++ b/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} +} diff --git a/examples/co2_props.jl b/examples/co2_props.jl index f999a8ed..5ad13989 100644 --- a/examples/co2_props.jl +++ b/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. diff --git a/examples/consistent_avgmpfa.jl b/examples/consistent_avgmpfa.jl new file mode 100644 index 00000000..45cfeaea --- /dev/null +++ b/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 diff --git a/examples/mpfa_weno_discretizations.jl b/examples/mpfa_weno_discretizations.jl new file mode 100644 index 00000000..db415948 --- /dev/null +++ b/examples/mpfa_weno_discretizations.jl @@ -0,0 +1,232 @@ +# # Alternative discretizations: WENO, NTPFA and AvgMPFA +# The following example demonstrates how to set up a reservoir model with +# different discretizations for flow and transport for a skewed grid. The model +# setup is taken from 6.1.2 in the MRST book [mrst-book-i](@cite). +# +# The example is a two-dimensional model with a single layer. The domain itself +# is completely symmetric, with producer in either lower corner and a single +# injector in the middle of the domain. From a flow perspective, the model is +# completely symmetric, and exactly the same results should be obtained for both +# producer wells. The grid is intentionally distorted to be skewed towards one +# producer, which will lead to consistency issues for the standard two-point +# flux approximation used by the solvers. +# +# The fluid model is also quite simple, with two immiscible phases that have +# linear relative permeability functions. A consequence of this is that the +# fluid fronts are not self-sharpening and are expected to be smeared out due to +# the numerical diffusion of the standard first-point upwind scheme. +using Jutul, JutulDarcy, GLMakie +Darcy, kg, meter, day, bar = si_units(:darcy, :kg, :meter, :day, :bar) +ny = 20 +nx = 2*ny + 1 +nz = 1 + +gcart = CartesianMesh((nx, ny, nz), (2.0, 1.0, 1.0)) +g = UnstructuredMesh(gcart) + +for (i, pt) in enumerate(g.node_points) + x, y, z = pt + pt += [0.4*(1.0-(x-1.0)^2)*(1.0-y), 0, 0] + g.node_points[i] = pt .* [1000.0, 1000.0, 1.0] +end + +c_i = cell_index(g, (nx÷2+1, ny, 1)) +c_p1 = cell_index(g, (1, 1, 1)) +c_p2 = cell_index(g, (nx, 1, 1)) + +fig = Figure() +ax = Axis(fig[1, 1]) +Jutul.plot_mesh_edges!(ax, g) +plot_mesh!(ax, g, cells = c_i, color = :red) +plot_mesh!(ax, g, cells = [c_p1, c_p2], color = :blue) +fig +# ## Define wells, fluid system and controls +reservoir = reservoir_domain(g, permeability = 0.1Darcy, porosity = 0.1) + +c_i1 = cell_index(g, (nx÷2+1, ny, 1)) + +I1 = setup_vertical_well(reservoir, nx÷2+1, ny, name = :I1, simple_well = true) +P1 = setup_vertical_well(reservoir, 1, 1, name = :P1, simple_well = true) +P2 = setup_vertical_well(reservoir, nx, 1, name = :P2, simple_well = true) + +phases = (AqueousPhase(), LiquidPhase()) +rhoWS = rhoLS = 1000.0 +rhoS = [rhoWS, rhoLS] .* kg/meter^3 +sys = ImmiscibleSystem(phases, reference_densities = rhoS) + +dt = repeat([30.0]*day, 100) +pv = pore_volume(reservoir) +inj_rate = sum(pv)/sum(dt) + +rate_target = TotalRateTarget(inj_rate) +I_ctrl = InjectorControl(rate_target, [1.0, 0.0], density = rhoWS) +bhp_target = BottomHolePressureTarget(50*bar) +P_ctrl = ProducerControl(bhp_target) +controls = Dict() +controls[:I1] = I_ctrl +controls[:P1] = P_ctrl +controls[:P2] = P_ctrl + +# ## Define functions to perform the simulations +# We are going to perform a number of simulations with different discretizations +# and it is therefore convenient to setup functions that can be called with the +# desired discretizations as arguments. The results are stored in a dictionary +# for easy retrieval after the simulations are done. +all_results = Dict() + +# ### Function to perform simulation +# The function `simulate_with_discretizations` sets up the reservoir model with +# the requested discretizations. The default behavior mirrors the default +# behavior of JutulDarcy by using the industry standard single-point upwind, +# two-point flux approximation. +function simulate_with_discretizations(upwind = :spu, kgrad = :tpfa) + model, parameters = setup_reservoir_model(reservoir, sys, + kgrad = kgrad, + upwind = upwind, + wells = [P1, P2, I1] + ) + kr = BrooksCoreyRelativePermeabilities(sys, 1.0, 0.0, 1.0) + replace_variables!(model, RelativePermeabilities = kr) + forces = setup_reservoir_forces(model, control = controls) + state0 = setup_reservoir_state(model, Pressure = 100*bar, Saturations = [0.0, 1.0]) + result = simulate_reservoir(state0, model, dt, info_level = 0, parameters = parameters, forces = forces); +end +# ### Function to plot the results +# We create a function `plot_discretization_result` that takes the result from a +# simultation and plots the saturation field after 75 steps. +function plot_discretization_result(result, name) + ws, states = result + sg = states[75][:Saturations][1, :] + fig = Figure() + ax = Axis(fig[1, 1], title = name) + plot_cell_data!(ax, g, sg, colormap = :seaborn_icefire_gradient) + Jutul.plot_mesh_edges!(ax, g) + return fig +end +# ### Function to simulate and plot discretizations +# The function `simulate_and_plot_discretizations` is a convenience function +# that calls the `simulate_with_discretizations` function and then plots the result. + +function simulate_and_plot_discretizations(upwind, kgrad) + result = simulate_with_discretizations(upwind, kgrad) + ustr = uppercase("$upwind") + if kgrad == :avgmpfa + kgstr = "AverageMPFA" + else + kgstr = uppercase("$kgrad") + end + descr = "$ustr with $kgstr" + all_results[descr] = result + plot_discretization_result(result, descr) +end + +# ### Simulate SPU-TPFA +# We start by simulating the model with the standard single-point upwind and the +# two-point flux approximation. This is the default discretization used by +# JutulDarcy. The schemes are robust (easy to implement and obtain nonlinear +# convergence) but suffer for consistency issues on skewed and non-K-orthogonal +# grids. We can observe both significant smearing (the fluid front is not sharp) +# and a significant bias towards the producer in the lower right corner. +simulate_and_plot_discretizations(:spu, :tpfa) +# ### Simulate SPU-AvgMPFA +# The next simulation uses a type of multi-point flux approximation (AverageMPFA +# or AvgMPFA), which is one class of consistent discretizations. As the scheme +# is consistent, the effect of the skewed grid is significantly reduced, with +# the remaining error being due to the coarse mesh chosen for illustrative +# purposes. Using AverageMPFA does not alleviate the smearing of the fluid +# front, however. +simulate_and_plot_discretizations(:spu, :avgmpfa) +# ### Simulate SPU-TPFA +# The next simulation uses a nonlinear two-point flux approximation (NTPFA). +# This scheme is consistent and monotone, and can be derived from the +# AverageMPFA scheme by allowing the ratio between the two half-face fluxes at +# each interface to vary based on pressure. +simulate_and_plot_discretizations(:spu, :ntpfa) +# ### Simulate WENO-TPFA +# The penultimate simulation uses a high-resolution scheme (WENO) for the flow, +# but reverts back to the inconsistent TPFA scheme for the pressure gradient. We +# now see significantly sharper fronts, but the bias towards the producer in the +# lower right corner is now back. +simulate_and_plot_discretizations(:weno, :tpfa) +# ### Simulate WENO-AvgMPFA +# The final simulation combines the high-resolution WENO scheme with the +# consistent AverageMPFA scheme. The result is a sharp fluid front with no bias +# towards either producer. This is the most accurate simulation of the four for +# this intentionally badly gridded problem. +simulate_and_plot_discretizations(:weno, :avgmpfa) +# ## Compare the results +# We can now compare the results of the different discretizations. We start by +# plotting the well curves in the two producers. Recall that the model is, from +# a flow perspective, completely symmetric, and the water cut in the two +# producers should be identical if the problem is fully resolved. We observe +# this for the consistent schemes, and delayed breakthrough for the SPU solves. +fig = Figure() +ax1 = Axis(fig[1, 1], title = "P1 water cut") +ax2 = Axis(fig[2, 1], title = "P2 water cut") +colors = Makie.wong_colors() +for (i, pr) in enumerate(all_results) + name, result = pr + ws, states = result + wcut1 = ws[:P1][:wcut] + wcut2 = ws[:P2][:wcut] + lines!(ax1, wcut1, label = "$name", color = colors[i]) + lines!(ax2, wcut2, label = "$name", color = colors[i]) +end +axislegend(ax1, position = :lt) +fig +# ## Compare the saturation fields as contours +# We can also compare the saturation fields for the different discretizations to +# see the spatial differences. To make life a bit easier, we make a function to +# do the comparison plots for us. +function compare_contours(name1, name2, title) + cmp = :seaborn_icefire_gradient + r1 = all_results[name1] + r2 = all_results[name2] + fig = Figure(size = (1200, 800)) + ax = Axis(fig[1, 1], title = title) + s1 = reshape(r1.states[75][:Saturations][1, :], nx, ny) + s2 = reshape(r2.states[75][:Saturations][1, :], nx, ny) + + contourf!(ax, s2, colormap = cmp, label = name2) + contour!(ax, s1, color = :white, linewidth = 5) + contour!(ax, s1, colormap = cmp, linewidth = 3, label = name1) + axislegend(position = :ct) + + axd = Axis(fig[1, 2], title = "Difference") + plt = contourf!(axd, s1-s2, + colormap = :seismic, + label = name2, + levels = range(-1, 1, length = 10), + colorscale = (-1.0, 1.0) + ) + Colorbar(fig[1, 3], limits = (-1, 1), colormap = :seismic) + Colorbar(fig[2, 1:3], limits = (-1, 1), colormap = cmp, vertical = false) + + ax1 = Axis(fig[3, 1], title = name1) + contourf!(ax1, s1, colormap = cmp) + + ax2 = Axis(fig[3, 2], title = name2) + contourf!(ax2, s2, colormap = cmp) + fig +end +# ### Compare SPU-TPFA and SPU-AvgMPFA +compare_contours("SPU with AverageMPFA", "SPU with TPFA", "Consistent (AvgMPFA) vs inconsistent (TPFA)") +# ### Compare WENO-TPFA and SPU-TPFA +compare_contours("WENO with TPFA", "SPU with TPFA", "High resolution (WENO) vs first order (SPU)") +# ### Compare WENO-AvgMPFA and SPU-TPFA +compare_contours("WENO with AverageMPFA", "SPU with TPFA", "High resolution + consistent vs first order + inconsistent") +# ## Conclusion +# The example demonstrates how different discretizations can affect the results +# of a simulation. The example is intentionally set up to show the differences +# in the discretizations. The results show that the consistent discretizations +# (AvgMPFA and NTPFA) can reduce the bias from grid orientation, and, while not +# present here, permeability tensor effects. The high-resolution WENO scheme can +# be employed to reduce smearing of the fluid fronts, which is important for +# problems that are not self-sharpening (e.g. compositional flow and miscible +# displacements). +# +# These schemes are not a pancea, however, as they require more computational +# effort to linearize and can have robustness issues where nonlinear convergence +# rates deteriorate and shorter timesteps may be required. They are accessible +# in JutulDarcy through the high-level interface and can be applied to any mesh +# and physics combination supported by the rest of the solvers. diff --git a/src/cpr.jl b/src/cpr.jl index 2483bbc6..19db4f54 100644 --- a/src/cpr.jl +++ b/src/cpr.jl @@ -563,8 +563,13 @@ function update_weights!(cpr, cpr_storage::CPRStorage, model, res_storage, stora wno = nc+i w_i = view(w, :, wno) well = wr_map.wells[i] - acc_i = storage[well].state.TotalMasses - true_impes!(w_i, acc_i, r, 1, ncomp, ps, scaling) + wstate = storage[well].state + if strategy == :analytical + cpr_weights_no_partials!(w_i, model, wstate, nothing, 1, ncomp, scaling) + else + acc_i = wstate.TotalMasses + true_impes!(w_i, acc_i, r, 1, ncomp, ps, scaling) + end end end return w diff --git a/src/ext/ext_amgx.jl b/src/ext/ext_amgx.jl index d59884be..1f602aa4 100644 --- a/src/ext/ext_amgx.jl +++ b/src/ext/ext_amgx.jl @@ -59,7 +59,8 @@ end function JutulDarcy.gpu_update_preconditioner!(cpr::AMGXCPR, lsys, model, storage, recorder, executor, krylov, J_bsr, r_cu, op) T = eltype(J_bsr) - @tic "CPU cpr work" Jutul.update_preconditioner!(cpr, lsys, model, storage, recorder, executor, update_system_precond = false, T = T) + ctx = model.context + @tic "CPU cpr work" Jutul.update_preconditioner!(cpr, lsys, ctx, model, storage, recorder, executor, update_system_precond = false, T = T) # Transfer pressure system to GPU @tic "update system precond" JutulDarcy.gpu_update_preconditioner!(cpr.system_precond, lsys, model, storage, recorder, executor, krylov, J_bsr, r_cu, op) @tic "update pressure system" JutulDarcy.update_amgx_pressure_system!(cpr.pressure_precond, cpr.storage.A_p, eltype(J_bsr), cpr, recorder) diff --git a/src/ext/ext_cuda.jl b/src/ext/ext_cuda.jl index f6c8ead2..96c2157a 100644 --- a/src/ext/ext_cuda.jl +++ b/src/ext/ext_cuda.jl @@ -18,6 +18,7 @@ end function Jutul.linear_solve!(lsys::Jutul.LSystem, krylov::CUDAReservoirKrylov{Tv, Ti}, + context, model, storage = nothing, dt = nothing, diff --git a/src/facility/wells/wells.jl b/src/facility/wells/wells.jl index 1e944b2c..3c05f89a 100644 --- a/src/facility/wells/wells.jl +++ b/src/facility/wells/wells.jl @@ -108,7 +108,11 @@ function setup_well(g, K, reservoir_cells::AbstractVector; centers = cell_centers[:, reservoir_cells] if isnothing(reference_depth) - reference_depth = centers[3, 1] + if size(centers, 1) == 2 + reference_depth = 0.0 + else + reference_depth = centers[3, 1] + end end volumes = zeros(n) WI_computed = zeros(n) diff --git a/src/flux.jl b/src/flux.jl index dfe12b01..a987de27 100644 --- a/src/flux.jl +++ b/src/flux.jl @@ -127,7 +127,6 @@ pressure_gradient(state, disc) = gradient(state.Pressure, disc) return F(up) end - @inline function upwind(upw::SPU, X::AbstractArray, q) flag = q < zero(q) if flag @@ -143,6 +142,10 @@ end return upwind(upw, F, q) end +@inline function upwind(upw::Jutul.WENO.WENOFaceDiscretization, F, q) + return Jutul.WENO.weno_upwind(upw, F, q) +end + @inline capillary_gradient(::Nothing, c_l, c_r, ph, ph_ref) = 0.0 @inline function capillary_gradient(pc, c_l, c_r, ph, ph_ref) if ph == ph_ref diff --git a/src/forces/bc.jl b/src/forces/bc.jl index f93e1f83..3c8f0cc8 100644 --- a/src/forces/bc.jl +++ b/src/forces/bc.jl @@ -49,6 +49,13 @@ function FlowBoundaryCondition( dir = findfirst(isequal(dir), (:x, :y, :z)) end dist = cell_dims(G, cell)[dir] + # Approximate area since we don't know the face. + A = 1.0 + for i in 1:D + if i != dir + A *= cell_dims(G, cell)[i] + end + end K = domain[:permeability] cond = domain[:rock_thermal_conductivity] @@ -61,7 +68,7 @@ function FlowBoundaryCondition( # Take the diagonal ki = Jutul.expand_perm(ki, D)[dir, dir] # Distance to boundary is half the cell width - return ki*(dist/2.0) + return A*ki/(dist/2.0) end T_flow = local_trans(K) diff --git a/src/multiphase.jl b/src/multiphase.jl index a173b896..d4c825f1 100644 --- a/src/multiphase.jl +++ b/src/multiphase.jl @@ -492,7 +492,7 @@ function cpr_weights_no_partials!(w, model::SimulationModel{R, S}, state, r, n, tb = minbatch(model.context, nc) M = global_map(model.domain) density = Jutul.active_view(ρ, M, for_variables = false) - @batch minbatch = tb for i in axes(w, 2) + @batch minbatch = tb for i in 1:n for ph in axes(w, 1) @inbounds w[ph, i] = 1/value(density[ph, i]) end diff --git a/src/porousmedia.jl b/src/porousmedia.jl index 6f59bfcc..ebc36c29 100644 --- a/src/porousmedia.jl +++ b/src/porousmedia.jl @@ -156,7 +156,10 @@ end function discretized_domain_tpfv_flow(domain::Jutul.DataDomain; general_ad = false, kgrad = nothing, - upwind = nothing + upwind = nothing, + weno_threshold = 0.0, + weno_do_clamp = false, + weno_epsilon = 1e-10 ) N = domain[:neighbors] g = physical_representation(domain) @@ -189,6 +192,19 @@ function discretized_domain_tpfv_flow(domain::Jutul.DataDomain; else @assert isnothing(kgrad) || kgrad isa AbstractVector end + if upwind isa Symbol + if upwind == :spu + upwind = nothing + elseif upwind == :weno + upwind = Jutul.WENO.weno_discretize(domain, + threshold = weno_threshold, + do_clamp = weno_do_clamp, + epsilon = weno_epsilon + ) + else + error("Unknown upwind scheme $upwind, must be :spu or :weno") + end + end if general_ad ad_flag = :generic else diff --git a/src/utils.jl b/src/utils.jl index 662b7fe4..bcfe47a6 100644 --- a/src/utils.jl +++ b/src/utils.jl @@ -191,6 +191,9 @@ reservoir and that facility. # Keyword arguments + +## Basic model setup + - `wells=[]`: Vector of wells (e.g. from [`setup_well`](@ref)) that are to be used in the model. Each well must have a unique name. - `extra_out=true`: Return both the model and the parameters instead of just the @@ -206,6 +209,9 @@ reservoir and that facility. meshes with anisotropic perm or non-orthogonal cells than `:tpfa`. - `:ntpfa` gives a consistent nonlinear MPFA scheme (nonlinear version of `:avgmpfa` that preserves monotonicity) +- `upwind=nothing`: Type of upwinding to use. Can be `:spu` or `nothing` for + standard upwinding or `:weno` for a second-order weighted essentially + non-oscillatory scheme. - `extra_outputs=Symbol[]`: Extra output variables for reservoir model. Defaults to "typical" values seen in reservoir simulation. Valid values: Vector of symbols to be output, `true` for all variables and `false` for the minimal set @@ -238,6 +244,8 @@ impact simulation speed. - `general_ad=false`: Use more general form of AD. Will result in slower execution speed than if set to true, but can be useful when working with custom discretizations. +- `discretization_arg=NamedTuple()`: Additional keyword arguments passed onto + `discretized_domain_tpfv_flow` when setting up discretizations. ## Increment and variable options These options govern the range of values and the maximum allowable change of @@ -313,9 +321,11 @@ function setup_reservoir_model(reservoir::DataDomain, system::JutulSystem; nthreads = Threads.nthreads(), minbatch = 1000, kgrad = nothing, + upwind = nothing, immutable_model = false, wells_systems = missing, - wells_as_cells = false + wells_as_cells = false, + discretization_arg = NamedTuple() ) # Deal with wells, make sure that multisegment wells come last. if !(wells isa AbstractArray) @@ -355,7 +365,8 @@ function setup_reservoir_model(reservoir::DataDomain, system::JutulSystem; system; context = reservoir_context, general_ad = general_ad, - kgrad = kgrad + kgrad = kgrad, + upwind = upwind ) if thermal rmodel = add_thermal_to_model!(rmodel) @@ -603,6 +614,8 @@ end - `initial_dt=si_unit(:day)`: initial timestep in seconds (one day by default) - `target_ds=Inf`: target saturation change over a timestep used by timestepper. +- `target_dz=Inf`: target mole fraction change over a timestep used by + timestepper (compositional only). - `target_its=8`: target number of nonlinear iterations per time step - `offset_its=1`: dampening parameter for time step selector where larger values lead to more pessimistic estimates. @@ -663,6 +676,7 @@ function setup_reservoir_simulator(case::JutulCase; rtol = nothing, initial_dt = si_unit(:day), target_ds = Inf, + target_dz = Inf, target_its = 8, offset_its = 1, tol_cnv = 1e-3, @@ -740,6 +754,13 @@ function setup_reservoir_simulator(case::JutulCase; ) push!(sel, t_sat) end + if isfinite(target_dz) + @assert mode == :default "target_dz is only supported in serial." + t_sat = VariableChangeTimestepSelector( + :OverallMoleFractions, target_dz, relative = false, reduction = :max, model = :Reservoir + ) + push!(sel, t_sat) + end else @assert isnothing(timesteps) || timesteps == :none end diff --git a/test/discretizations.jl b/test/discretizations.jl new file mode 100644 index 00000000..ce3ea327 --- /dev/null +++ b/test/discretizations.jl @@ -0,0 +1,20 @@ +using Test, JutulDarcy +@testset "Discretizations" begin + pth = JutulDarcy.GeoEnergyIO.test_input_file_path("SPE1", "SPE1.DATA") + subs = 1:2 + for kgrad in [:tpfa, :tpfa_test, :avgmpfa, :ntpfa] + for upw in [:spu, :weno] + @testset "$kgrad-$upw" begin + case = setup_case_from_data_file(pth, + kgrad = kgrad, + upwind = upw, + ) + _, states = simulate_reservoir(case[subs], + info_level = -1, + failure_cuts_timestep = false + ) + @test length(states) == length(subs) + end + end + end +end diff --git a/test/gpu.jl b/test/gpu.jl index 2e490fe9..877eb8de 100644 --- a/test/gpu.jl +++ b/test/gpu.jl @@ -70,16 +70,17 @@ if CUDA.functional() @testset "High-level CUDA" begin spe1_pth = JutulDarcy.GeoEnergyIO.test_input_file_path("SPE1", "SPE1.DATA") case = setup_case_from_data_file(spe1_pth) + sim_kwarg = (info_level = -1, failure_cuts_timestep = false) @testset "CuSPARSE-ILU(0)" begin - res_ilu = simulate_reservoir(case, info_level = -1, precond = :ilu0) - res_cuilu = simulate_reservoir(case, linear_solver_backend = :cuda, precond = :ilu0, info_level = -1); + res_ilu = simulate_reservoir(case; precond = :ilu0, sim_kwarg...) + res_cuilu = simulate_reservoir(case; linear_solver_backend = :cuda, precond = :ilu0, sim_kwarg...); spe1_gpu_compare(res_ilu, res_cuilu) end if Sys.islinux() @testset "AMGX-CPR" begin using AMGX - res_cpr = simulate_reservoir(case, info_level = -1, precond = :cpr) - res_cucpr = simulate_reservoir(case, linear_solver_backend = :cuda, precond = :cpr, info_level = -1); + res_cpr = simulate_reservoir(case; precond = :cpr, sim_kwarg...) + res_cucpr = simulate_reservoir(case; linear_solver_backend = :cuda, precond = :cpr, sim_kwarg...); spe1_gpu_compare(res_cpr, res_cucpr) end end diff --git a/test/runtests.jl b/test/runtests.jl index 6fe06e39..1e412d63 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -60,6 +60,13 @@ end include("co2props.jl") end +@testitem "Discretizations" begin + include("discretizations.jl") +end + +@testitem "GPU" begin + include("gpu.jl") +end + @run_package_tests nothing -# include("gpu.jl")