Merge pull request #11 from slimgroup/fixinit

set up initial state
slimgroup · Mar 8, 2023 · fe2dc98 · fe2dc98
2 parents ba2856c + 4e9b801
commit fe2dc98
Show file tree

Hide file tree

Showing 3 changed files with 20 additions and 7 deletions.
diff --git a/src/FlowAD/Operators/jutulModeling.jl b/src/FlowAD/Operators/jutulModeling.jl
@@ -22,14 +22,10 @@ function (S::jutulModeling{D, T})(LogTransmissibilities::AbstractVector{T}, f::j
     model.models.Reservoir.domain.grid.trans .= Transmissibilities
     parameters[:Reservoir][:Transmissibilities] = Transmissibilities
 
-    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, report = simulate!(sim, tstep, forces = forces, config = config, max_timestep_cuts = 1000, info_level=info_level);
     output = jutulStates(states)
     return output

diff --git a/src/FlowAD/Types/jutulState.jl b/src/FlowAD/Types/jutulState.jl
@@ -35,6 +35,10 @@ Pressure(state::jutulSimpleState) = state.state[:Pressure]
 Saturations(state::jutulSimpleOrMultiModelStates) = vcat([Saturations(state.states[i]) for i = 1:get_nt(state)]...)
 Pressure(state::jutulSimpleOrMultiModelStates) = vcat([Pressure(state.states[i]) for i = 1:get_nt(state)]...)
 
+get_Reservoir_state(state::jutulState) = state.state[:Reservoir]
+get_Reservoir_state(state::jutulSimpleState) = state.state
+get_Reservoir_state(state::jutulSimpleOrMultiModelStates) = get_Reservoir_state(state.states[end])
+
 get_nt(state::jutulSimpleOrMultiModelStates) = length(state.states)
 get_nn(state::jutulSimpleOrMultiModelState) = length(Saturations(state))
 get_nn(state::jutulSimpleOrMultiModelStates) = get_nn(state.states[1])

diff --git a/test/test_jutulModeling.jl b/test/test_jutulModeling.jl
@@ -20,7 +20,20 @@ end
     states = S(x, q1)
     for i = 1:length(states.states)
         exist_co2 = sum(Saturations(states.states[i]) .* states.states[i].state[:PhaseMassDensities][1,:] .* model.ϕ) * prod(model.d)
-        inj_co2 = JutulDarcyAD.ρCO2 * q.irate * JutulDarcyAD.day * sum(tstep[1:i])
+        inj_co2 = JutulDarcyAD.ρCO2 * q.irate * JutulDarcyAD.day * sum(S.tstep[1:i])
         @test isapprox(exist_co2, inj_co2) rtol=1e-3
     end
 end
+
+@testset "Test mass conservation for well modeling, different injection rates" begin
+    states = S(x, q)
+    pre_co2 = sum(Saturations(states.states[end]) .* states.states[end].state[:Reservoir][:PhaseMassDensities][1,:] .* model.ϕ) * prod(model.d)
+    q2 = jutulForce(q.irate * 0.5, q.loc)
+    S.tstep ./= 2
+    states_end = S(x, q2; state0=states)
+    for i = 1:length(states_end.states)
+        exist_co2 = sum(Saturations(states_end.states[i]) .* states_end.states[i].state[:Reservoir][:PhaseMassDensities][1,:] .* model.ϕ) * prod(model.d)
+        inj_co2 = JutulDarcyAD.ρCO2 * q2.irate * JutulDarcyAD.day * sum(S.tstep[1:i])
+        @test isapprox(exist_co2-pre_co2, inj_co2) rtol=1e-3
+    end
+end