diff --git a/Project.toml b/Project.toml index e0fb615..ca94bd7 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.2" +version = "0.2.3" [deps] ChainRulesCore = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4" diff --git a/src/FlowRules/Operators/jutulModeling.jl b/src/FlowRules/Operators/jutulModeling.jl index 2aff594..3630fe4 100644 --- a/src/FlowRules/Operators/jutulModeling.jl +++ b/src/FlowRules/Operators/jutulModeling.jl @@ -8,7 +8,7 @@ end display(M::jutulModeling{D, T}) where {D, T} = println("$(D)D jutulModeling structure with $(sum(M.tstep)) days in $(length(M.tstep)) time steps") -function (S::jutulModeling{D, T})(LogTransmissibilities::AbstractVector{T}, f::Union{jutulForce{D, N}, jutulVWell{D, N}}; +function (S::jutulModeling{D, T})(LogTransmissibilities::AbstractVector{T}, ϕ::AbstractVector{T}, f::Union{jutulForce{D, N}, jutulVWell{D, N}}; state0=nothing, visCO2::T=T(visCO2), visH2O::T=T(visH2O), ρCO2::T=T(ρCO2), ρH2O::T=T(ρH2O), info_level::Int64=-1) where {D, T, N} @@ -20,7 +20,9 @@ function (S::jutulModeling{D, T})(LogTransmissibilities::AbstractVector{T}, f::U ### 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) .* ϕ parameters[:Reservoir][:Transmissibilities] = Transmissibilities + parameters[:Reservoir][:FluidVolume] .= prod(S.model.d) .* ϕ isnothing(state0) || (state0_[:Reservoir] = get_Reservoir_state(state0)) @@ -31,8 +33,8 @@ function (S::jutulModeling{D, T})(LogTransmissibilities::AbstractVector{T}, f::U return output end -function (S::jutulModeling{D, T})(LogTransmissibilities::AbstractVector{T}, f::jutulSource{D, N}; - state0::jutulSimpleState{T}=jutulSimpleState(S.model), visCO2::T=T(visCO2), visH2O::T=T(visH2O), +function (S::jutulModeling{D, T})(LogTransmissibilities::AbstractVector{T}, ϕ::AbstractVector{T}, f::jutulSource{D, N}; + state0=nothing, visCO2::T=T(visCO2), visH2O::T=T(visH2O), ρCO2::T=T(ρCO2), ρH2O::T=T(ρH2O), info_level::Int64=-1) where {D, T, N} Transmissibilities = exp.(LogTransmissibilities) @@ -43,7 +45,26 @@ function (S::jutulModeling{D, T})(LogTransmissibilities::AbstractVector{T}, f::j 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) .* ϕ parameters = setup_parameters(model, PhaseViscosities = [visCO2, visH2O]); - states, _ = simulate(dict(state0), model, tstep, parameters = parameters, forces = forces, info_level = info_level, max_timestep_cuts = 1000) + 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) return jutulSimpleStates(states) +end + +function (S::jutulModeling{D, T})(f::Union{jutulForce{D, N}, jutulVWell{D, N}, jutulSource{D, N}}; + LogTransmissibilities::AbstractVector{T}=KtoTrans(CartesianMesh(S.model), S.model.K), ϕ::AbstractVector{T}=S.model.ϕ, + state0=nothing, visCO2::T=T(visCO2), visH2O::T=T(visH2O), + ρ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 + +function (S::jutulModeling{D, T})(LogTransmissibilities::AbstractVector{T}, f::Union{jutulForce{D, N}, jutulVWell{D, N}, jutulSource{D, N}}; + ϕ::AbstractVector{T}=S.model.ϕ, + state0=nothing, visCO2::T=T(visCO2), visH2O::T=T(visH2O), + ρ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 diff --git a/src/FlowRules/Operators/rrule.jl b/src/FlowRules/Operators/rrule.jl index 2e43a48..88e8e59 100644 --- a/src/FlowRules/Operators/rrule.jl +++ b/src/FlowRules/Operators/rrule.jl @@ -1,7 +1,7 @@ -function rrule(S::jutulModeling{D, T}, LogTransmissibilities::AbstractVector{T}, f::Union{jutulForce{D, N}, jutulVWell{D, N}}; +function rrule(S::jutulModeling{D, T}, LogTransmissibilities::AbstractVector{T}, ϕ::AbstractVector{T}, f::Union{jutulForce{D, N}, jutulVWell{D, N}}; state0=nothing, visCO2::T=T(visCO2), visH2O::T=T(visH2O), ρCO2::T=T(ρCO2), ρH2O::T=T(ρH2O), info_level::Int64=-1) where {D, T, N} - + Transmissibilities = exp.(LogTransmissibilities) ### set up simulation time @@ -10,16 +10,14 @@ 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) .* ϕ parameters[:Reservoir][:Transmissibilities] = Transmissibilities + parameters[:Reservoir][:FluidVolume] .= prod(S.model.d) .* ϕ - if isnothing(state0) - state0 = state0_ - else - state0 = dict(state0) - end + isnothing(state0) || (state0_[:Reservoir] = get_Reservoir_state(state0)) ### simulation - sim, config = setup_reservoir_simulator(model, state0, parameters); + sim, config = setup_reservoir_simulator(model, state0_, parameters); states, reports = simulate!(sim, tstep, forces = forces, config = config, max_timestep_cuts = 1000, info_level=info_level); output = jutulStates(states) @@ -33,15 +31,15 @@ function rrule(S::jutulModeling{D, T}, LogTransmissibilities::AbstractVector{T}, states_ref = dict(states_ref_) mass_mismatch = (m, state, dt, step_no, forces) -> loss_per_step(m, state, dt, step_no, forces, states_ref) 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); + 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)], NoTangent() + return NoTangent(), g[1:length(LogTransmissibilities)], g[length(LogTransmissibilities)+1:length(LogTransmissibilities)+prod(S.model.n)] * prod(S.model.d), NoTangent() end return output, pullback end -function rrule(S::jutulModeling{D, T}, LogTransmissibilities::AbstractVector{T}, f::jutulSource{D, N}; - state0::jutulSimpleState{T}=jutulSimpleState(S.model), visCO2::T=T(visCO2), visH2O::T=T(visH2O), +function rrule(S::jutulModeling{D, T}, LogTransmissibilities::AbstractVector{T}, ϕ::AbstractVector{T}, f::jutulSource{D, N}; + state0=nothing, visCO2::T=T(visCO2), visH2O::T=T(visH2O), ρCO2::T=T(ρCO2), ρH2O::T=T(ρH2O), info_level::Int64=-1) where {D, T, N} Transmissibilities = exp.(LogTransmissibilities) @@ -52,10 +50,13 @@ function rrule(S::jutulModeling{D, T}, LogTransmissibilities::AbstractVector{T}, 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) .* ϕ parameters = setup_parameters(model, PhaseViscosities = [visCO2, visH2O]); - states, reports = simulate(dict(state0), model, tstep, parameters = parameters, forces = forces, info_level = info_level, max_timestep_cuts = 1000) + 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) output = jutulSimpleStates(states) - cfg = optimization_config(model, parameters, use_scaling = true, rel_min = 0.1, rel_max = 10) + cfg = optimization_config(model, parameters, use_scaling = false, rel_min = 0., rel_max = Inf) for (ki, vi) in cfg if ki in [:TwoPointGravityDifference, :PhaseViscosities] vi[:active] = false @@ -64,7 +65,6 @@ function rrule(S::jutulModeling{D, T}, LogTransmissibilities::AbstractVector{T}, vi[:scaler] = :log end end - cfg[:Transmissibilities][:use_scaling] = false function pullback(dy) states_dy = output(dy) @@ -82,9 +82,9 @@ function rrule(S::jutulModeling{D, T}, LogTransmissibilities::AbstractVector{T}, mass_mismatch = (m, state, dt, step_no, forces) -> loss_per_step_simple(m, state, dt, step_no, forces, states_ref) Jutul.evaluate_objective(mass_mismatch, model, states_ref, tstep, forces) 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); + 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)], NoTangent() + return NoTangent(), g[1:length(LogTransmissibilities)], g[length(LogTransmissibilities)+1:length(LogTransmissibilities)+prod(S.model.n)] * prod(S.model.d), NoTangent() end return output, pullback end diff --git a/test/test_gradient.jl b/test/test_gradient.jl index 7e36093..f29a15c 100644 --- a/test/test_gradient.jl +++ b/test/test_gradient.jl @@ -6,33 +6,36 @@ S = jutulModeling(model0, tstep) ## simulation x = log.(KtoTrans(CartesianMesh(model), model.K)) x0 = log.(KtoTrans(CartesianMesh(model0), model0.K)) +ϕ = S.model.ϕ -states = S(x, q) +states = S(x, ϕ, q) -misfit(x0) = 0.5 * norm(S(x0, q) - states).^2 -g = gradient(()->misfit(x0), Flux.params(x0)) +misfit(x0, ϕ, q, states) = 0.5 * norm(S(x0, ϕ, q) - states).^2 +g = gradient(()->misfit(x0, ϕ, q, states), Flux.params(x0, ϕ)) dx = randn(MersenneTwister(2023), length(x0)) dx = dx/norm(dx) * norm(x0)/5.0 +dϕ = randn(MersenneTwister(2023), length(ϕ)) +dϕ = dϕ/norm(dϕ) * norm(ϕ)/2.75e11 + @testset "Taylor-series gradient test of jutulModeling with wells" begin - grad_test(misfit, x0, dx, g[x0]) + grad_test(x0->misfit(x0, ϕ, q, states), x0, dx, g[x0]) + grad_test(ϕ->misfit(x0, ϕ, q, states), ϕ, dϕ, g[ϕ]) end -states1 = S(x, q1) - -misfit1(x0) = 0.5 * norm(S(x0, q1) - states1).^2 -g1 = gradient(()->misfit1(x0), Flux.params(x0)) +states1 = S(x, ϕ, q1) +g1 = gradient(()->misfit(x0, ϕ, q1, states1), Flux.params(x0, ϕ)) @testset "Taylor-series gradient test of simple jutulModeling" begin - grad_test(misfit1, x0, dx/1.5, g1[x0]) + grad_test(x0->misfit(x0, ϕ, q1, states1), x0, dx/1.5, g1[x0]) + grad_test(ϕ->misfit(x0, ϕ, q1, states1), ϕ, dϕ/1e2, g1[ϕ]) end states2 = S(x, q2) - -misfit2(x0) = 0.5 * norm(S(x0, q2) - states2).^2 -g2 = gradient(()->misfit2(x0), Flux.params(x0)) +g2 = gradient(()->misfit(x0, ϕ, q2, states2), Flux.params(x0, ϕ)) @testset "Taylor-series gradient test of jutulModeling with vertical wells" begin - grad_test(misfit2, x0, dx/1.5, g2[x0]) + grad_test(x0->misfit(x0, ϕ, q2, states2), x0, dx/1.5, g2[x0]) + grad_test(ϕ->misfit(x0, ϕ, q2, states2), ϕ, dϕ/4551., g2[ϕ]) end