Merge pull request #12 from slimgroup/verticalwell

vertical wells

Project.toml

@@ -1,7 +1,7 @@
 name = "JutulDarcyAD"
 uuid = "41f0c4f5-9bdd-4ef1-8c3a-d454dff2d562"
 authors = ["Ziyi Yin <ziyi.yin@gatech.edu>"]
-version = "0.2.0"
+version = "0.2.1"
 
 [deps]
 ChainRulesCore = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4"

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)

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

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")

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}
 

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}
 

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)

src/FlowAD/Types/type_utils.jl

@@ -2,17 +2,33 @@ 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]
     end
     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

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

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

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"

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

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)