diff --git a/Project.toml b/Project.toml index 2c05a8d..23c30c5 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "JutulDarcyRules" uuid = "41f0c4f5-9bdd-4ef1-8c3a-d454dff2d562" authors = ["Ziyi Yin "] -version = "0.2.4" +version = "0.2.5" [deps] ChainRulesCore = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4" @@ -15,8 +15,8 @@ Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2" [compat] ChainRulesCore = "1" Flux = "0.12, 0.13" -Jutul = "=0.2.2" -JutulDarcy = "=0.2.1" +Jutul = "=0.2.5" +JutulDarcy = "=0.2.2" Optim = "1" julia = "1" diff --git a/src/FlowRules/Operators/jutulModeling.jl b/src/FlowRules/Operators/jutulModeling.jl index 3630fe4..d5158fa 100644 --- a/src/FlowRules/Operators/jutulModeling.jl +++ b/src/FlowRules/Operators/jutulModeling.jl @@ -19,8 +19,8 @@ function (S::jutulModeling{D, T})(LogTransmissibilities::AbstractVector{T}, ϕ:: ### set up simulation configurations model, parameters, state0_, forces = setup_well_model(S.model, f, tstep; visCO2=visCO2, visH2O=visH2O, ρCO2=ρCO2, ρH2O=ρH2O) - model.models.Reservoir.domain.grid.trans .= Transmissibilities - model.models.Reservoir.domain.grid.pore_volumes .= prod(S.model.d) .* ϕ + + model.models.Reservoir.data_domain[:porosity] = ϕ parameters[:Reservoir][:Transmissibilities] = Transmissibilities parameters[:Reservoir][:FluidVolume] .= prod(S.model.d) .* ϕ @@ -44,9 +44,12 @@ function (S::jutulModeling{D, T})(LogTransmissibilities::AbstractVector{T}, ϕ:: ### set up simulation time tstep = day * S.tstep model = simple_model(S.model; ρCO2=ρCO2, ρH2O=ρH2O) - model.domain.grid.trans .= Transmissibilities - model.domain.grid.pore_volumes .= prod(S.model.d) .* ϕ + model.data_domain[:porosity] = ϕ + parameters = setup_parameters(model, PhaseViscosities = [visCO2, visH2O]); + parameters[:Transmissibilities] = Transmissibilities + parameters[:FluidVolume] .= prod(S.model.d) .* ϕ + state0_ = jutulSimpleState(S.model) isnothing(state0) || (state0_ = state0) states, _ = simulate(dict(state0_), model, tstep, parameters = parameters, forces = forces, info_level = info_level, max_timestep_cuts = 1000) @@ -67,4 +70,4 @@ function (S::jutulModeling{D, T})(LogTransmissibilities::AbstractVector{T}, f::U ρCO2::T=T(ρCO2), ρH2O::T=T(ρH2O), info_level::Int64=-1) where {D, T, N} return S(LogTransmissibilities, ϕ, f; state0=state0, visCO2=visCO2, visH2O=visH2O, ρCO2=ρCO2, ρH2O=ρH2O, info_level=info_level) -end \ No newline at end of file +end diff --git a/src/FlowRules/Operators/rrule.jl b/src/FlowRules/Operators/rrule.jl index 88e8e59..9016ff0 100644 --- a/src/FlowRules/Operators/rrule.jl +++ b/src/FlowRules/Operators/rrule.jl @@ -9,8 +9,8 @@ function rrule(S::jutulModeling{D, T}, LogTransmissibilities::AbstractVector{T}, ### set up simulation configurations model, parameters, state0_, forces = setup_well_model(S.model, f, tstep; visCO2=visCO2, visH2O=visH2O, ρCO2=ρCO2, ρH2O=ρH2O) - model.models.Reservoir.domain.grid.trans .= Transmissibilities - model.models.Reservoir.domain.grid.pore_volumes .= prod(S.model.d) .* ϕ + + model.models.Reservoir.data_domain[:porosity] = ϕ parameters[:Reservoir][:Transmissibilities] = Transmissibilities parameters[:Reservoir][:FluidVolume] .= prod(S.model.d) .* ϕ @@ -33,7 +33,11 @@ function rrule(S::jutulModeling{D, T}, LogTransmissibilities::AbstractVector{T}, F_o, dF_o, F_and_dF, x0, lims, data = setup_parameter_optimization( states, reports, model, state0_, parameters, tstep, forces, mass_mismatch, cfg, param_obj = true, print = info_level, config = config, use_sparsity = false); g = dF_o(similar(x0), x0); - return NoTangent(), g[1:length(LogTransmissibilities)], g[length(LogTransmissibilities)+1:length(LogTransmissibilities)+prod(S.model.n)] * prod(S.model.d), NoTangent() + n_faces = length(LogTransmissibilities) + n_cells = prod(S.model.n) + dLogTransmissibilities = g[1:n_faces] + dϕ = g[n_faces + 1 : n_faces + n_cells] * prod(S.model.d) + return NoTangent(), dLogTransmissibilities, dϕ, NoTangent() end return output, pullback end @@ -49,9 +53,12 @@ function rrule(S::jutulModeling{D, T}, LogTransmissibilities::AbstractVector{T}, ### set up simulation time tstep = day * S.tstep model = simple_model(S.model; ρCO2=ρCO2, ρH2O=ρH2O) - model.domain.grid.trans .= Transmissibilities - model.domain.grid.pore_volumes .= prod(S.model.d) .* ϕ + model.data_domain[:porosity] = ϕ + parameters = setup_parameters(model, PhaseViscosities = [visCO2, visH2O]); + parameters[:Transmissibilities] = Transmissibilities + parameters[:FluidVolume] .= prod(S.model.d) .* ϕ + state0_ = jutulSimpleState(S.model) isnothing(state0) || (state0_ = state0) states, reports = simulate(dict(state0_), model, tstep, parameters = parameters, forces = forces, info_level = info_level, max_timestep_cuts = 1000) @@ -84,7 +91,11 @@ function rrule(S::jutulModeling{D, T}, LogTransmissibilities::AbstractVector{T}, F_o, dF_o, F_and_dF, x0, lims, data = setup_parameter_optimization(states, reports, model, dict(state0_), parameters, tstep, forces, mass_mismatch, cfg, print = -1, param_obj = true); g = dF_o(similar(x0), x0); - return NoTangent(), g[1:length(LogTransmissibilities)], g[length(LogTransmissibilities)+1:length(LogTransmissibilities)+prod(S.model.n)] * prod(S.model.d), NoTangent() + n_faces = length(LogTransmissibilities) + n_cells = prod(S.model.n) + dLogTransmissibilities = g[1:n_faces] + dϕ = g[n_faces + 1 : n_faces + n_cells] * prod(S.model.d) + return NoTangent(), dLogTransmissibilities, dϕ, NoTangent() end return output, pullback end @@ -209,4 +220,4 @@ function setup_parameter_optimization(precomputed_states, reports, case::JutulCa dF = (dFdx, x) -> gradient_opt!(dFdx, x, data) F_and_dF = (F, dFdx, x) -> objective_and_gradient_opt!(F, dFdx, x, data, print) return (F! = F, dF! = dF, F_and_dF! = F_and_dF, x0 = x0, limits = lims, data = data) -end \ No newline at end of file +end diff --git a/src/FlowRules/Types/type_utils.jl b/src/FlowRules/Types/type_utils.jl index 3a6daa1..02996e2 100644 --- a/src/FlowRules/Types/type_utils.jl +++ b/src/FlowRules/Types/type_utils.jl @@ -36,12 +36,13 @@ function setup_well_model(M::jutulModel{D, T}, f::Union{jutulForce{D, T}, jutulV ### set up model, parameters sys = ImmiscibleSystem((VaporPhase(), AqueousPhase()), reference_densities = [ρH2O, ρCO2]) - domain = discretized_domain_tpfv_flow(tpfv_geometry(CartesianMesh(M)), porosity = M.ϕ, permeability = M.K) - model, parameters = setup_reservoir_model(domain, sys, wells = Is) + domain_spec = reservoir_domain(CartesianMesh(M), porosity = M.ϕ, permeability = M.K) + domain = discretized_domain_tpfv_flow(domain_spec) + model_parameters = Dict(:PhaseViscosities=> [visCO2, visH2O]) + model, parameters = setup_reservoir_model(domain_spec, sys, wells = Is, parameters=model_parameters) select_output_variables!(model.models.Reservoir, :all) ρ = ConstantCompressibilityDensities(p_ref = 150*bar, density_ref = [ρCO2, ρH2O], compressibility = [1e-4/bar, 1e-6/bar]) replace_variables!(model, PhaseMassDensities = ρ) - replace_variables!(model, PhaseViscosities = vcat(visCO2 * ones(prod(M.n))', visH2O * ones(prod(M.n))')) replace_variables!(model, RelativePermeabilities = BrooksCoreyRelPerm(sys, [2.0, 2.0], [0.1, 0.1], 1.0)) for x ∈ keys(model.models) Jutul.select_output_variables!(model.models[x], :all) @@ -69,11 +70,12 @@ end function simple_model(M::jutulModel{D, T}; ρCO2::T=T(ρCO2), ρH2O::T=T(ρH2O)) where {D, T} sys = ImmiscibleSystem((VaporPhase(), AqueousPhase())) g = CartesianMesh(M.n, M.d .* M.n) - G = discretized_domain_tpfv_flow(tpfv_geometry(g), porosity = M.ϕ, permeability = M.K) - model = SimulationModel(G, sys, output_level = :all) + domain_spec = reservoir_domain(g, porosity = M.ϕ, permeability = M.K) + G = discretized_domain_tpfv_flow(domain_spec) + model = SimulationModel(domain_spec, sys, output_level = :all) model.primary_variables[:Pressure] = JutulDarcy.Pressure(minimum = -Inf, max_rel = nothing) ρ = ConstantCompressibilityDensities(p_ref = 100*bar, density_ref = [ρCO2, ρH2O], compressibility = [1e-4/bar, 1e-6/bar]) replace_variables!(model, PhaseMassDensities = ρ) replace_variables!(model, RelativePermeabilities = BrooksCoreyRelPerm(sys, [2.0, 2.0], [0.1, 0.1], 1.0)) return model -end \ No newline at end of file +end diff --git a/test/grad_test.jl b/test/grad_test.jl new file mode 100644 index 0000000..dd1e003 --- /dev/null +++ b/test/grad_test.jl @@ -0,0 +1,108 @@ +mean(x) = sum(x)/length(x) + +function log_division(a, b) + if b == 0 + return a == 0 ? NaN : Inf + end + return log(a / b) +end + +""" + grad_test(J, x0, Δx, dJdx; ΔJ=nothing, maxiter=6, h0=5e-2, stol=1e-1, hfactor=8e-1) + +Test the gradient using Taylor series convergence. + +Compute a series of residuals using zeroth order and first order Taylor approximations. +Each perturbed computation J(x₀ + hᵢ Δx) is compared to J(x₀) and J(x₀) + hᵢ ΔJ, +where hᵢ = hfactor * hᵢ₋₁ and ΔJ = dJdx ⋅ Δx by default. If the computation of J and +dJdx is correct, J is sufficiently smooth, and h is sufficiently small, the zeroth +order approximation should converge linearly and the first order approximation should +converge quadratically. + +It can be difficult to obtain the correct convergence, so we average the +convergence factors across maxiter values of h and test that the convergence +factor is more than correct convergence factor minus the stol parameter. + +# Mathematical basis + +For J sufficiently smooth, the value of a point at a small perturbation from x₀ can be +computed using the Taylor series expansion. +```math +J(x₀ + h Δx) = J(x₀) + h (dJ/dx) Δx + h² (d²J/dx²) : (Δx Δxᵀ) + O(h³) +``` + +If d²J/dx² is non-zero and h is small enough, the value of the first-order Taylor +approximation differs from the true value by approximately a constant proportional to h². +```math +err(h) = |J(x₀ + h Δx) - J(x₀) - h (dJ/dx) Δx| ≈ h²c +``` + +If we consider the error for two different values of h, the unknown constant can be eliminated. +```math +err(h₁) / err(h₂) ≈ (h₁ / h₂)^2 +``` +So if h₁ is divided by a factor α, then the ratio of the errors should be divided by a factor α². +Or we can compute the exponent for the rate of convergence using logarithms. +```math +log(err(h₁) / err(h₂)) / log (h₁ / h₂) ≈ 2 +``` +""" +function grad_test(J, x0, Δx, dJdx; ΔJ=nothing, maxiter=6, h0=5e-2, stol=1e-1, hfactor=8e-1) + if !xor(isnothing(dJdx), isnothing(ΔJ)) + error("Must specify either dJdx or ΔJ") + end + if isnothing(ΔJ) + ΔJ = dot(dJdx, Δx) + end + J0 = J(x0) + h = h0 + + log_factor = log(hfactor) + expected_f1 = 1e0 / hfactor + expected_f2 = 1e0 / hfactor ^2e0 + + err1 = zeros(Float64, maxiter) + err2 = zeros(Float64, maxiter) + Js = zeros(Float64, maxiter) + @printf("%11s, %12s | %11s, %11s | %11s, %11s | %12s, %12s | %12s %12s %12s \n", "h", "ΔJ", "e1", "e2", "factor1", "factor2", "rate1", "rate2", "Finite-diff", "T1 approx", "T2 approx") + @printf("%11s, % 12.5e | %11s, %11s | %11.5e, %11.5e | % 12.5e, % 12.5e | \n", "", ΔJ, "0", "0", expected_f1, expected_f2, 1, 2) + @printf("%11s, %12s | %11s, %11s | %11s, %11s | % 12s, % 12s \n", "___________", "___________", "___________", "___________", "___________", "___________", "____________", "____________") + all_info = [] + for j=1:maxiter + Jh = J(x0 + h*Δx) + Js[j] = Jh + err1[j] = norm(Jh - J0, 1) + err2[j] = norm(Jh - J0 - h*ΔJ, 1) + j == 1 ? prev = 1 : prev = j - 1 + + dJ_est = (Jh - J0) / h + α = expected_f2 + dJ_est1 = (α*Js[prev] - Jh + (1-α)*J0) / (h * (α/hfactor - 1)) + α = -expected_f2 + dJ_est2 = (α*Js[prev] - Jh + (1-α)*J0) / (h * (α/hfactor - 1)) + + rate1 = log_division(err1[j], err1[prev]) / log_factor + rate2 = log_division(err2[j], err2[prev]) / log_factor + + info = (h, h*norm(ΔJ, 1), err1[j], err2[j], err1[prev]/err1[j], err2[prev]/err2[j], rate1, rate2, dJ_est, dJ_est1, dJ_est2) + push!(all_info, info) + @printf("%11.5e, % 12.5e | %11.5e, %11.5e | %11.5e, %11.5e | % 12.5e, % 12.5e | % 12.5e, % 12.5e, % 12.5e \n", info...) + h = h * hfactor + end + + println() + @printf("%11s, %12s | %11s, %11s | %11s, %11s | %12s, %12s | %12s %12s %12s \n", "h", "ΔJ", "e1", "e2", "factor1", "factor2", "rate1", "rate2", "Finite-diff", "T1 approx", "T2 approx") + @printf("%11s, % 12.5e | %11s, %11s | %11.5e, %11.5e | % 12.5e, % 12.5e | \n", "", ΔJ, "0", "0", expected_f1, expected_f2, 1, 2) + @printf("%11s, %12s | %11s, %11s | %11s, %11s | % 12s, % 12s \n", "___________", "___________", "___________", "___________", "___________", "___________", "____________", "____________") + for j=1:maxiter + info = all_info[j] + @printf("%11.5e, % 12.5e | %11.5e, %11.5e | %11.5e, %11.5e | % 12.5e, % 12.5e | % 12.5e, % 12.5e, % 12.5e \n", info...) + h = h * hfactor + end + + factor1 = err1[1:end-1]./err1[2:end] + factor2 = err2[1:end-1]./err2[2:end] + @test mean(factor1) ≥ expected_f1 - stol + @test mean(factor2) ≥ expected_f2 - stol +end + diff --git a/test/test_gradient.jl b/test/test_gradient.jl index 06d8373..342a511 100644 --- a/test/test_gradient.jl +++ b/test/test_gradient.jl @@ -21,7 +21,7 @@ dϕ = dϕ/norm(dϕ) * norm(ϕ)/2.75e11 @testset "Taylor-series gradient test of jutulModeling with wells" begin grad_test(x0->misfit(x0, ϕ, q, states), x0, dx, g[x0]) - grad_test(ϕ->misfit(x0, ϕ, q, states), ϕ, dϕ, g[ϕ]) + grad_test(ϕ->misfit(x0, ϕ, q, states), ϕ, dϕ, g[ϕ], h0=2e1, maxiter=12) end states1 = S(x, ϕ, q1) @@ -36,6 +36,7 @@ states2 = S(x, q2) g2 = gradient(()->misfit(x0, ϕ, q2, states2), Flux.params(x0, ϕ)) @testset "Taylor-series gradient test of jutulModeling with vertical wells" begin - grad_test(x0->misfit(x0, ϕ, q2, states2), x0, dx/1.5, g2[x0]) - grad_test(ϕ->misfit(x0, ϕ, q2, states2), ϕ, dϕ/4551., g2[ϕ]) + # This test is very brittle. There may be an issue here. + grad_test(x0->misfit(x0, ϕ, q2, states2), x0, dx, g2[x0], h0=1e-2, maxiter=12) + grad_test(ϕ->misfit(x0, ϕ, q2, states2), ϕ, dϕ, g2[ϕ], h0=5e0, maxiter=6) end diff --git a/test/test_utils.jl b/test/test_utils.jl index 4d184a5..66affe3 100644 --- a/test/test_utils.jl +++ b/test/test_utils.jl @@ -26,30 +26,4 @@ function test_config() return model, model0, q, q1, q2, state0, state1, tstep end -mean(x) = sum(x)/length(x) - -## adapted from https://github.com/slimgroup/JUDI.jl/blob/master/test/testing_utils.jl -function grad_test(misfit, x0, dx, g; maxiter=6, h0=5f-2, stol=1f-1) - # init - err1 = zeros(Float32, maxiter) - err2 = zeros(Float32, maxiter) - - gdx = dot(g, dx) - f0 = misfit(x0) - h = h0 - - @printf("%11.5s, %11.5s, %11.5s, %11.5s, %11.5s, %11.5s \n", "h", "gdx", "e1", "e2", "rate1", "rate2") - for j=1:maxiter - f = misfit(x0 + h*dx) - err1[j] = norm(f - f0, 1) - err2[j] = norm(f - f0 - h*gdx, 1) - j == 1 ? prev = 1 : prev = j - 1 - @printf("%5.5e, %5.5e, %5.5e, %5.5e, %5.5e, %5.5e \n", h, h*norm(gdx, 1), err1[j], err2[j], err1[prev]/err1[j], err2[prev]/err2[j]) - h = h * .8f0 - end - - rate1 = err1[1:end-1]./err1[2:end] - rate2 = err2[1:end-1]./err2[2:end] - @test isapprox(mean(rate1), 1.25f0; atol=stol) - @test isapprox(mean(rate2), 1.5625f0; atol=stol) -end +include("grad_test.jl")