Skip to content

Commit

Permalink
Merge pull request #15 from slimgroup/poro
Browse files Browse the repository at this point in the history 
Add porosity into the simulation w/ rrules
  • Loading branch information
@ziyiyin97
ziyiyin97 authored Mar 15, 2023
2 parents a94e085 + f77c28d commit 24f7269
Show file tree
Hide file tree
Showing 4 changed files with 59 additions and 35 deletions.
2 changes: 1 addition & 1 deletion Project.toml
View file
Original file line number Diff line number Diff line change
@@ -1,7 +1,7 @@
name = "JutulDarcyRules"
uuid = "41f0c4f5-9bdd-4ef1-8c3a-d454dff2d562"
authors = ["Ziyi Yin <ziyi.yin@gatech.edu>"]
version = "0.2.2"
version = "0.2.3"

This comment has been minimized.

Sign in to view

Copy link
@ziyiyin97

ziyiyin97 Mar 15, 2023

Author Member

@JuliaRegistrator register

[deps]
ChainRulesCore = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4"
Expand Down
29 changes: 25 additions & 4 deletions src/FlowRules/Operators/jutulModeling.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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}

Expand All @@ -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))

Expand All @@ -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)
Expand All @@ -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
34 changes: 17 additions & 17 deletions src/FlowRules/Operators/rrule.jl
View file
Original file line number Diff line number Diff line change
@@ -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
Expand All @@ -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)

Expand All @@ -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)
Expand All @@ -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
Expand All @@ -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)
Expand All @@ -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
Expand Down
29 changes: 16 additions & 13 deletions test/test_gradient.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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

1 comment on commit 24f7269

@JuliaRegistrator
Copy link

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Registration pull request created: JuliaRegistries/General/79664

After the above pull request is merged, it is recommended that a tag is created on this repository for the registered package version.

This will be done automatically if the Julia TagBot GitHub Action is installed, or can be done manually through the github interface, or via:

git tag -a v0.2.3 -m "<description of version>" 24f7269283e1212c0f3cae96640885008e31d9de
git push origin v0.2.3

Please sign in to comment.