Merge pull request #4 from slimgroup/fixerror

Restructure jutulState; add mass conservation tests

Project.toml

@@ -1,7 +1,7 @@
 name = "JutulDarcyAD"
 uuid = "41f0c4f5-9bdd-4ef1-8c3a-d454dff2d562"
 authors = ["Ziyi Yin <[email protected]>"]
-version = "0.1.1"
+version = "0.1.2"
 
 [deps]
 ChainRulesCore = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4"

src/FlowAD/Operators/jutulModeling.jl

@@ -9,8 +9,8 @@ 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};
-    state0::jutulInitState{T}=jutulInitState(S.model), visCO2::T=T(1e-4), visH2O::T=T(1e-3),
-    ρCO2::T=T(501.9), ρH2O::T=T(1053.0), info_level::Int64=-1) where {D, T, N}
+    state0::jutulState{T}=jutulState(S.model), 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)
     forces = force(S.model, f; ρCO2=ρCO2)

src/FlowAD/Operators/rrule.jl

@@ -1,14 +1,14 @@
 function rrule(S::jutulModeling{D, T}, LogTransmissibilities::AbstractVector{T}, f::jutulForce{D, N};
-    state0::jutulInitState{T}=jutulInitState(S.model), visCO2::T=T(1e-4), visH2O::T=T(1e-3),
-    ρCO2::T=T(501.9), ρH2O::T=T(1053.0), info_level::Int64=-1) where {D, T, N}
+    state0::jutulState{T}=jutulState(S.model), 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)
     forces = force(S.model, f; ρCO2=ρCO2)
     tstep = day * S.tstep
     model = model_(S.model)
     model.domain.grid.trans .= Transmissibilities
     parameters = setup_parameters(model, PhaseViscosities = [visCO2, visH2O], density = [ρCO2, ρH2O]); # 0.1 and 1 cP
-    states, rep = simulate(dict(state0), model, tstep, parameters = parameters, forces = forces, info_level = info_level, max_timestep_cuts = 1000)
+    states, _ = simulate(dict(state0), model, tstep, parameters = parameters, forces = forces, info_level = info_level, max_timestep_cuts = 1000)
     output = jutulStates(states)
     cfg = optimization_config(model, parameters, use_scaling = true, rel_min = 0.1, rel_max = 10)
     for (ki, vi) in cfg

src/FlowAD/Types/jutulForce.jl

@@ -10,7 +10,7 @@ jutulForce(irate::T, loc::Vector{NTuple{D, T}}) where {D, T} = jutulForce([(-T(1
 display(f::jutulForce{D, T}) where {D, T} =
     println("$(D)D jutulForce structure with $(length(f.loc)) injection/production wells and rate $(f.irate) m^3/s")
 
-function force(M::jutulModel{D, T}, f::jutulForce{D, T}; ρCO2::T=T(501.9)) where {D, T}
+function force(M::jutulModel{D, T}, f::jutulForce{D, T}; ρCO2::T=T(ρCO2)) where {D, T}
     model = model_(M)
     cell_loc = [Int.(round.(f.loc[i] ./ M.d)) for i = 1:length(f.loc)]
     cell = [sum([(cell_loc[i][d]-1) * prod(M.n[1:d-1]) for d = length(cell_loc[i]):-1:1]) + 1 for i = 1:length(cell_loc)]

src/FlowAD/Types/jutulModel.jl

@@ -15,7 +15,7 @@ CartesianMesh(M::jutulModel{D, T}) where {D, T} = CartesianMesh(M.n, M.d .* M.n)
 function model_(M::jutulModel{D, T}) where {D, T}
     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)
+    model = SimulationModel(G, sys, output_level = :all)
     replace_variables!(model, RelativePermeabilities = kr)
     return model
 end

src/FlowAD/Types/jutulState.jl

@@ -1,8 +1,11 @@
-export jutulState, jutulInitState, jutulStates, dict
+export jutulState, jutulStates, dict
 
 abstract type jutulAllState{T} <: DenseVector{T} end
 
 struct jutulState{T} <: jutulAllState{T}
+    PhaseMassMobilities::Matrix{T}
+    PhaseMassDensities::Matrix{T}
+    RelativePermeabilities::Matrix{T}
     Saturations::Vector{T}
     Pressure::Vector{T}
     TotalMasses::Matrix{T}
@@ -12,40 +15,29 @@ struct jutulStates{T} <: jutulAllState{T}
     states::Vector{jutulState{T}}
 end
 
-struct jutulInitState{T} <: jutulAllState{T}
-    PhaseMassMobilities::Matrix{T}
-    PhaseMassDensities::Matrix{T}
-    RelativePermeabilities::Matrix{T}
-    Saturations::Vector{T}
-    Pressure::Vector{T}
-    TotalMasses::Matrix{T}
-end
-
 display(state::jutulAllState{T}) where T = println("$(typeof(state))")
 
-jutulState(state::Dict{Symbol, T}) where T = jutulState(state[:Saturations][1,:], state[:Pressure], state[:TotalMasses])
+jutulState(state::Dict{Symbol, T}) where T =
+    jutulState(state[:PhaseMassMobilities], state[:PhaseMassDensities], state[:RelativePermeabilities], 
+    state[:Saturations][1,:], state[:Pressure], state[:TotalMasses])
 
 jutulStates(states::Vector{Dict{Symbol, T}}) where T = jutulStates([jutulState(states[i]) for i = 1:length(states)])
 
-jutulInitState(state::Dict{Symbol, T}) where T = jutulInitState(state[:PhaseMassMobilities],
-    state[:PhaseMassDensities], state[:RelativePermeabilities], state[:Saturations][1,:],
-    state[:Pressure], state[:TotalMasses])
-
-function jutulInitState(M::jutulModel{D, T}; ρH2O::T=T(1053.0), g::T=T(10.0)) where {D, T}
+function jutulState(M::jutulModel{D, T}; ρH2O::T=T(ρH2O), g::T=T(10.0)) where {D, T}
     ## default state at time 0 with all water
     Z = repeat((1:M.n[end])*M.d[end], inner = prod(M.n[1:2]))
     p0 = ρH2O * g * Z # rho * g * h
     state0 = setup_state(model_(M), Pressure = p0, Saturations = [0.0, 1.0])
-    return jutulInitState(state0)
+    return jutulState(state0)
 end
 
 get_nt(state::jutulStates) = length(state.states)
-get_nn(state::jutulAllState) = length(state.Saturations)
+get_nn(state::jutulState) = length(state.Saturations)
 get_nn(state::jutulStates) = get_nn(state.states[1])
 
 ###### turn jutulStates to state dictionary
 
-function dict(state::jutulAllState{T}) where T
+function dict(state::jutulState{T}) where T
     dict_ = Dict(fieldnames(typeof(state)) .=> getfield.(Ref(state), fieldnames(typeof(state))))
     dict_[:Saturations] = hcat(dict_[:Saturations], T(1) .- dict_[:Saturations])'
     return dict_
@@ -55,15 +47,15 @@ dict(state::jutulStates{T}) where T = [dict(state.states[i]) for i = 1:get_nt(st
 
 ###### AbstractVector
 
-length(state::jutulAllState) = length(state.Saturations) + length(state.Pressure)
+length(state::jutulState) = length(state.Saturations) + length(state.Pressure)
 length(state::jutulStates) = sum([length(state.states[i]) for i = 1:get_nt(state)])
 size(state::jutulAllState) = (length(state),)
 
 vec(state::jutulStates) = vcat(vcat([state.states[i].Saturations for i = 1:get_nt(state)]...), vcat([state.states[i].Pressure for i = 1:get_nt(state)]...))
-vec(state::jutulAllState) = vcat(state.Saturations, state.Pressure)
+vec(state::jutulState) = vcat(state.Saturations, state.Pressure)
 
-IndexStyle(state::jutulAllState) = IndexLinear()
-function getindex(state::jutulAllState, i::Int)
+IndexStyle(::jutulAllState) = IndexLinear()
+function getindex(state::jutulState, i::Int)
     if i <= get_nn(state)
         ## getindex for saturation
         return state.Saturations[i]
@@ -85,7 +77,7 @@ function getindex(state::jutulStates, i::Int)
     end
 end
 
-function setindex!(state::jutulAllState, v, i::Int)
+function setindex!(state::jutulState, v, i::Int)
     if i <= get_nn(state)
         ## setindex for saturation
         state.Saturations[i] = v

src/JutulDarcyAD.jl

@@ -19,6 +19,11 @@ module JutulDarcyAD
 
     JutulDarcyADPATH = dirname(pathof(JutulDarcyAD))
 
+    visCO2 = 1e-4
+    visH2O = 1e-3
+    ρCO2 = 501.9
+    ρH2O = 1053.0
+
     const Darcy = 9.869232667160130e-13
     const md = Darcy * 1e-3
 

test/runtests.jl

@@ -12,8 +12,10 @@ Random.seed!(2023)
 include("test_utils.jl")
 
 include("test_conversion.jl")
+
 include("test_jutulState.jl")
 include("test_jutulForce.jl")
 include("test_jutulModel.jl")
+include("test_jutulModeling.jl")
 
 include("test_gradient.jl")

test/test_gradient.jl

@@ -2,7 +2,6 @@ model, model0, q, tstep = test_config()
 
 ## set up modeling operator
 S = jutulModeling(model0, tstep)
-S = jutulModeling(model0, tstep)
 
 ## simulation
 x = log.(KtoTrans(CartesianMesh(model), model.K))

test/test_jutulModeling.jl

@@ -0,0 +1,18 @@
+model, model0, q, tstep = test_config()
+
+## set up modeling operator
+S = jutulModeling(model, tstep)
+
+## simulation
+x = log.(KtoTrans(CartesianMesh(model), model.K))
+x0 = log.(KtoTrans(CartesianMesh(model0), model0.K))
+
+states = S(x, q)
+
+@testset "Test mass conservation" begin
+    for i = 1:length(states.states)
+        exist_co2 = sum(states.states[i].Saturations .* states.states[i].PhaseMassDensities[1,:]) * prod(model.d) * model.ϕ
+        inj_co2 = JutulDarcyAD.ρCO2 * q.irate[1] * JutulDarcyAD.day * sum(tstep[1:i])
+        @test isapprox(exist_co2, inj_co2) rtol=1e-3
+    end
+end

test/test_utils.jl

@@ -8,7 +8,7 @@ function setup_model_state()
     K = 200 * rand() * md * ones(n)
     ϕ = rand()
     model = jutulModel(n, d, ϕ, K1to3(K))
-    init_state = jutulInitState(model)
+    init_state = jutulState(model)
     return model, init_state
 
 end