diff --git a/Project.toml b/Project.toml index a0e4f87..a19771e 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "JutulDarcyAD" uuid = "41f0c4f5-9bdd-4ef1-8c3a-d454dff2d562" authors = ["Ziyi Yin "] -version = "0.2.0" +version = "0.2.1" [deps] ChainRulesCore = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4" diff --git a/examples/Vwell_2D.jl b/examples/Vwell_2D.jl new file mode 100644 index 0000000..b10baef --- /dev/null +++ b/examples/Vwell_2D.jl @@ -0,0 +1,57 @@ +## A simple 2D example for fluid-flow simulation + +using DrWatson +@quickactivate "JutulDarcyAD-example" + +using JutulDarcyAD +using LinearAlgebra +using PyPlot + +## grid size +n = (30, 1, 15) +d = (30.0, 30.0, 30.0) + +## permeability +K0 = 200 * md * ones(n) +K = deepcopy(K0) +K[:,:,1:2:end] .*= 100 + +ϕ = 0.25 +model = jutulModel(n, d, ϕ, K1to3(K)) + +## simulation time steppings +tstep = 50 * ones(10) +tot_time = sum(tstep) + +## injection & production +inj_loc = (15, 1, 10) .* d +irate = 5e-3 +q = jutulVWell(irate, (450., 30.); startz = 270., endz = 330.) + +## set up modeling operator +S = jutulModeling(model, tstep) + +## simulation +Trans = KtoTrans(CartesianMesh(model), K1to3(K)) +Trans0 = KtoTrans(CartesianMesh(model), K1to3(K0)) +@time states = S(log.(Trans), q) +@time states0 = S(log.(Trans0), q) + +## plotting +fig=figure(figsize=(20,12)); +subplot(1,2,1); +imshow(reshape(Saturations(states.states[end]), n[1], n[end])'); colorbar(); title("saturation") +subplot(1,2,2); +imshow(reshape(Pressure(states.states[end]), n[1], n[end])'); colorbar(); title("pressure") + +## plotting +fig=figure(figsize=(20,12)); +subplot(1,2,1); +imshow(reshape(Saturations(states0.states[end]), n[1], n[end])'); colorbar(); title("saturation") +subplot(1,2,2); +imshow(reshape(Pressure(states0.states[end]), n[1], n[end])'); colorbar(); title("pressure") + +exist_co2 = sum(Saturations(states.states[end]) .* states.states[end].state[:Reservoir][:PhaseMassDensities][1,:] .* model.ϕ) * prod(model.d) +inj_co2 = JutulDarcyAD.ρCO2 * q.irate * JutulDarcyAD.day * sum(tstep) + +norm(exist_co2-inj_co2)/norm(exist_co2+inj_co2) diff --git a/examples/compass.jl b/examples/compass.jl index e44bd8c..e385bc9 100644 --- a/examples/compass.jl +++ b/examples/compass.jl @@ -166,4 +166,4 @@ for j=1:niterations end -JLD2.@save "compass-inv.jld2" logK logK0 state +JLD2.@save "compass-inv.jld2" logK logK0 state \ No newline at end of file diff --git a/src/FlowAD/FlowAD.jl b/src/FlowAD/FlowAD.jl index 082dff6..7224ce1 100644 --- a/src/FlowAD/FlowAD.jl +++ b/src/FlowAD/FlowAD.jl @@ -1,6 +1,7 @@ ### types include("Types/jutulModel.jl") include("Types/jutulForce.jl") +include("Types/jutulVWell.jl") include("Types/jutulSource.jl") include("Types/jutulState.jl") include("Types/type_utils.jl") diff --git a/src/FlowAD/Operators/jutulModeling.jl b/src/FlowAD/Operators/jutulModeling.jl index 14d253b..2aff594 100644 --- a/src/FlowAD/Operators/jutulModeling.jl +++ b/src/FlowAD/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::jutulForce{D, N}; +function (S::jutulModeling{D, T})(LogTransmissibilities::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} diff --git a/src/FlowAD/Operators/rrule.jl b/src/FlowAD/Operators/rrule.jl index 6fa07c5..2e43a48 100644 --- a/src/FlowAD/Operators/rrule.jl +++ b/src/FlowAD/Operators/rrule.jl @@ -1,4 +1,4 @@ -function rrule(S::jutulModeling{D, T}, LogTransmissibilities::AbstractVector{T}, f::jutulForce{D, N}; +function rrule(S::jutulModeling{D, T}, LogTransmissibilities::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} diff --git a/src/FlowAD/Types/jutulVWell.jl b/src/FlowAD/Types/jutulVWell.jl new file mode 100644 index 0000000..1bd341b --- /dev/null +++ b/src/FlowAD/Types/jutulVWell.jl @@ -0,0 +1,17 @@ +export jutulVWell + +struct jutulVWell{D, T} + irate::T + name::Vector{Symbol} + loc::Vector + startz + endz +end + +jutulVWell(irate::T, loc::NTuple{D, T}; startz=nothing, endz=nothing) where {D, T} = jutulVWell(irate, [loc]; startz=startz, endz=endz) +jutulVWell(irate::T, loc::Vector{NTuple{D, T}}; startz=nothing, endz=nothing) where {D, T}= jutulVWell{D+1, T}(irate, vcat(:Injector, [:Producer for i = 1:length(loc)]), loc, startz, endz) + +display(f::jutulVWell{D, T}) where {D, T} = + println("$(D)D jutulVWell structure with $(length(f.loc)) injection/production wells and rate $(f.irate) m^3/s") + +==(A::jutulVWell{D, T}, B::jutulVWell{D, T}) where {D,T} = (A.irate == B.irate && A.name == B.name && A.loc == B.loc) \ No newline at end of file diff --git a/src/FlowAD/Types/type_utils.jl b/src/FlowAD/Types/type_utils.jl index 768d886..4c4ee95 100644 --- a/src/FlowAD/Types/type_utils.jl +++ b/src/FlowAD/Types/type_utils.jl @@ -2,9 +2,9 @@ function force(M::jutulModel{D, T}, w::jutulForce{D, T}, tstep::Vector{T}; ρCO2::T=T(ρCO2), ρH2O::T=T(ρH2O), g::T=T(10.0)) where {D, T} ## set up well information - cell_loc = [Int.(round.(w.loc[i] ./ M.d)) for i = 1:length(w.loc)] + cell_loc = [Int.(round.(w.loc[i] ./ M.d[1:length(w.loc[1])])) for i = 1:length(w.loc)] Is = [setup_well(CartesianMesh(M), M.K, [cell_loc[i]], name = w.name[i]) for i = 1:length(w.loc)] - ctrls = [w.name[i]==:Injector ? InjectorControl(TotalRateTarget(w.irate), [1.0, 0.0], density = ρCO2) : ProducerControl(BottomHolePressureTarget(2.0 * ρH2O * g * w.loc[i][3])) for i = 1:length(w.loc)] + ctrls = [w.name[i]==:Injector ? InjectorControl(TotalRateTarget(w.irate), [1.0, 0.0], density = ρCO2) : ProducerControl(BottomHolePressureTarget(50*bar)) for i = 1:length(w.loc)] controls = Dict() for i = 1:length(w.loc) controls[w.name[i]] = ctrls[i] @@ -12,7 +12,23 @@ function force(M::jutulModel{D, T}, w::jutulForce{D, T}, tstep::Vector{T}; return Is, controls end -function setup_well_model(M::jutulModel{D, T}, f::jutulForce{D, T}, tstep::Vector{T}; +function force(M::jutulModel{D, T}, w::jutulVWell{D, T}, tstep::Vector{T}; + ρCO2::T=T(ρCO2), ρH2O::T=T(ρH2O), g::T=T(10.0)) where {D, T} + + ## set up well information + cell_loc = [Int.(round.(w.loc[i] ./ M.d[1:length(w.loc[1])])) for i = 1:length(w.loc)] + heel = [isnothing(w.startz) ? 1 : Int(div(w.startz[i], M.d[end])) for i = 1:length(w.loc)] + toe = [isnothing(w.endz) ? M.n[end] : Int(div(w.endz[i], M.d[end])) for i = 1:length(w.loc)] + Is = [setup_vertical_well(CartesianMesh(M), M.K, cell_loc[i]..., name = w.name[i]; heel = heel[i], toe = toe[i]) for i = 1:length(w.loc)] + ctrls = [w.name[i]==:Injector ? InjectorControl(TotalRateTarget(w.irate), [1.0, 0.0], density = ρCO2) : ProducerControl(BottomHolePressureTarget(50*bar)) for i = 1:length(w.loc)] + controls = Dict() + for i = 1:length(w.loc) + controls[w.name[i]] = ctrls[i] + end + return Is, controls +end + +function setup_well_model(M::jutulModel{D, T}, f::Union{jutulForce{D, T}, jutulVWell{D, T}}, tstep::Vector{T}; visCO2::T=T(visCO2), visH2O::T=T(visH2O), ρCO2::T=T(ρCO2), ρH2O::T=T(ρH2O), g::T=T(10.0)) where {D, T} ### set up well controls diff --git a/test/test_gradient.jl b/test/test_gradient.jl index cc0f0ee..7e36093 100644 --- a/test/test_gradient.jl +++ b/test/test_gradient.jl @@ -1,4 +1,4 @@ -model, model0, q, q1, init_state, init_state1, tstep = test_config(); +model, model0, q, q1, q2, state0, state1, tstep = test_config(); ## set up modeling operator S = jutulModeling(model0, tstep) @@ -27,3 +27,12 @@ g1 = gradient(()->misfit1(x0), Flux.params(x0)) @testset "Taylor-series gradient test of simple jutulModeling" begin grad_test(misfit1, x0, dx/1.5, g1[x0]) end + +states2 = S(x, q2) + +misfit2(x0) = 0.5 * norm(S(x0, q2) - states2).^2 +g2 = gradient(()->misfit2(x0), Flux.params(x0)) + +@testset "Taylor-series gradient test of jutulModeling with vertical wells" begin + grad_test(misfit2, x0, dx/1.5, g2[x0]) +end diff --git a/test/test_jutulModeling.jl b/test/test_jutulModeling.jl index 6511bc1..b59fbe9 100644 --- a/test/test_jutulModeling.jl +++ b/test/test_jutulModeling.jl @@ -1,4 +1,4 @@ -model, model0, q, q1, init_state, init_state1, tstep = test_config(); +model, model0, q, q1, q2, init_state, init_state1, tstep = test_config(); ## set up modeling operator S = jutulModeling(model, tstep) @@ -37,3 +37,12 @@ end @test isapprox(exist_co2-pre_co2, inj_co2) rtol=1e-3 end end + +@testset "Test mass conservation for vertical well modeling" begin + states = S(x, q2) + for i = 1:length(states.states) + exist_co2 = sum(Saturations(states.states[i]) .* states.states[i].state[:Reservoir][:PhaseMassDensities][1,:] .* model.ϕ) * prod(model.d) + inj_co2 = JutulDarcyAD.ρCO2 * q.irate * JutulDarcyAD.day * sum(tstep[1:i]) + @test isapprox(exist_co2, inj_co2) rtol=1e-3 + end +end diff --git a/test/test_jutulState.jl b/test/test_jutulState.jl index d975d4f..ca1a97f 100644 --- a/test/test_jutulState.jl +++ b/test/test_jutulState.jl @@ -1,4 +1,4 @@ -model, model0, q, q1, init_state, init_state_simple, tstep = test_config(); +model, model0, q, q1, q2, init_state, init_state_simple, tstep = test_config(); @testset "Test jutulState" begin @info "length, size" diff --git a/test/test_jutulVWell.jl b/test/test_jutulVWell.jl new file mode 100644 index 0000000..ba13645 --- /dev/null +++ b/test/test_jutulVWell.jl @@ -0,0 +1,15 @@ +@testset "Test jutulVWell" begin + + d = (1.0, 20.0, 3.0) + inj_loc = (3, 1) .* d[1:2] + prod_loc = (28, 1) .* d[1:2] + irate = rand() + q = jutulVWell(irate, [inj_loc, prod_loc]) + q1 = jutulVWell(irate, [inj_loc]) + + @test q != q1 + @test q.irate == q1.irate + @test q.loc[1] == q1.loc[1] + @test q.name[1] == q1.name[1] + @test q.name[2] == :Producer +end \ No newline at end of file diff --git a/test/test_utils.jl b/test/test_utils.jl index 14bcd6f..8a508df 100644 --- a/test/test_utils.jl +++ b/test/test_utils.jl @@ -18,11 +18,12 @@ function test_config() inj_loc = (15, 1, 10) .* d prod_loc = (30, 1, 10) .* d irate = 5e-3 - q = jutulForce(irate, [inj_loc]) + q = jutulForce(irate, inj_loc) q1 = jutulSource(irate, [inj_loc, prod_loc]) + q2 = jutulVWell(irate, inj_loc[1:2]; startz = 9 * d[3], endz = 11 * d[3]) state0 = jutulState(JutulDarcyAD.setup_well_model(model, q, tstep)[3]) state1 = jutulSimpleState(model) - return model, model0, q, q1, state0, state1, tstep + return model, model0, q, q1, q2, state0, state1, tstep end mean(x) = sum(x)/length(x)