In [1]:
# load the package
include("../FractionalFlow/FractionalFlow.jl")
using PyPlot, SetPyPlot, Dierckx, Statistics, NLopt
import Calculus
import GR
FF = FractionalFlow
Out[1]:
In [2]:
# define the problem
# relative permeabilities
rel_perms = FF.oil_water_rel_perms(krw0=0.4, kro0=0.9,
swc=0.15, sor=0.2, nw=2.0, no = 2.0)
# FF.visualize(rel_perms)
# define the fluids
fluids = FF.oil_water_fluids(mu_water=1e-3, mu_oil=1e-3)
# define the fractional flow functions
fw, dfw = FF.fractional_flow_function(rel_perms, fluids)
# visualize the fractional flow
# FF.visualize(rel_perms, fluids, label="lowsal")
# tight_layout()
ut = 1.15e-5
phi = 0.3
L = 0.15
core_flood = FF.core_flooding(u_inj=1.15e-5, pv_inject=5,
p_back=1e5, sw_init=0.2, sw_inj=1.0, rel_perms=rel_perms)
core_props = FF.core_properties()
wf_res = FF.water_flood(core_props, fluids, rel_perms, core_flood)
fw, dfw = FF.fractional_flow_function(rel_perms, fluids)
sw_tmp = FF.linspace(0,1,100)
FF.visualize(wf_res)
Out[2]:
In [3]:
t_sec, pv, rec_fact, dp_core, x, sw_face, c_face, c_out_sal=FF.water_flood_numeric(core_props, fluids, rel_perms, core_flood, Nx=50)
Out[3]:
In [4]:
t_exp_dp = wf_res.dp_time[:,1]
dp_exp = wf_res.dp_time[:,2]
t_exp_R = wf_res.recovery_time[:,1]
R_exp = wf_res.recovery_time[:,2]
plotyy(t_exp_R, R_exp, t_exp_dp, dp_exp, fig_size = [8,5], x_label="time [s]", y1_label="R [-]", y2_label="dP [Pa]")
Out[4]:
In [5]:
plot(t_exp_dp, dp_exp, "o", t_sec, dp_core)
legend(["Analytical", "Numerical"])
Out[5]:
In [6]:
# struct
struct exp_data
t_exp_dp
dp_exp
t_exp_R
R_exp
end
exp_data1 = exp_data(t_exp_dp, dp_exp, t_exp_R, R_exp)
Out[6]:
In [7]:
"""
rel_perm_param [krw0, kro0, nw, no, swc, sor]
"""
function error_calc(rel_perm_param, exp_data, core_props, fluids, core_flood; w_p=1.0, w_R=1.0)
rel_perms = FF.oil_water_rel_perms(krw0=rel_perm_param[1], kro0=rel_perm_param[2],
swc=rel_perm_param[5], sor=rel_perm_param[6], nw=rel_perm_param[3], no = rel_perm_param[4])
wf_res = FF.water_flood(core_props, fluids, rel_perms, core_flood)
dp_calc = Spline1D(wf_res.dp_time[:,1], wf_res.dp_time[:,2], k=1, bc="nearest")
R_calc = Spline1D(wf_res.recovery_time[:,1], wf_res.recovery_time[:,2], k=1, bc="nearest")
error_dp = abs.(dp_calc(exp_data.t_exp_dp) .- exp_data.dp_exp)
# println(error_dp)
error_R = abs.(R_calc(exp_data.t_exp_R) .- exp_data.R_exp)
# println(error_R)
error_dp_norm = w_p.*error_dp./exp_data.dp_exp
error_R_norm = w_R.*error_R./(exp_data.R_exp.+eps())
return mean(error_R_norm)+mean(error_dp_norm)
end
function vis_error(rel_perm_param, exp_data, core_props, fluids, core_flood)
rel_perms = FF.oil_water_rel_perms(krw0=rel_perm_param[1], kro0=rel_perm_param[2],
swc=rel_perm_param[5], sor=rel_perm_param[6], nw=rel_perm_param[3], no = rel_perm_param[4])
wf_res = FF.water_flood(core_props, fluids, rel_perms, core_flood)
figure()
plot(wf_res.dp_time[:,1], wf_res.dp_time[:,2], exp_data.t_exp_dp, exp_data.dp_exp, "o")
xlabel("t [s]")
ylabel("dp [Pa]")
legend(["Theoretical", "Experiment"])
figure()
plot(wf_res.recovery_time[:,1], wf_res.recovery_time[:,2], exp_data.t_exp_R, exp_data.R_exp, "v")
xlabel("t [s]")
ylabel("R [-]")
legend(["Theoretical", "Experiment"])
end
# test
x_init = [0.109681, 0.201297, 3.96653, 3.0, 0.19, 0.262231]
vis_error(x_init, exp_data1, core_props, fluids, core_flood)
error_calc(x_init, exp_data1, core_props, fluids, core_flood)
Out[7]:
In [23]:
# weight factors:
w_p = ones(length(exp_data1.dp_exp))
temp_val, ind_max = findmax(exp_data1.dp_exp)
println(ind_max)
w_p[ind_max-3:ind_max+3] .= 10
w_p[end-10:end] .= 10
w_p[1]=10
w_R = ones(length(exp_data1.R_exp))
w_R[20:25] .= 10
w_R[end-3:end] .= 10
function f(x)
f_val = 0.0
try
f_val = error_calc(x, exp_data1, core_props, fluids, core_flood, w_p = w_p, w_R = w_R)
catch
f_val = 100.0
# info("Objective function did not converge!")
end
return f_val
end
function g(x)
eps1 = 1e-3
f_val = f(x)
g_val = ones(length(x))
try
# g_val = Calculus.gradient(x -> error_calc(x, exp_data1, core_props, fluids, core_flood), x)
for j in eachindex(x)
x2 = copy(x)
x2[j]+=eps1
f_val2 = f(x2)
g_val[j] = (f_val2-f_val)/eps1
end
catch
g_val = ones(length(x))
end
return g_val
end
function obj_fun(param, grad)
if length(grad)>0
grad[:] = g(param)
end
obj_fun_val = f(param)
if isnan(obj_fun_val) || isinf(obj_fun_val)
obj_fun_val = 100.0
end
return obj_fun_val
end
# test
grad_x = zeros(6)
obj_fun([1.0, 0.8, 3, 4, 0.2, 0.2], grad_x)
f([1.0, 0.8, 2, 2, 0.1, 0.2])
Out[23]:
In [9]:
grad_x
Out[9]:
In [10]:
## algorithms
# L: Local, G:global
# D: derivative-based, N: non-derivative (search-based)
# :LD_MMA
# :LN_COBYLA
# :LD_LBFGS
# :GN_DIRECT
# :GN_DIRECT_L
# GN_CRS2_LM
# G_MLSL_LDS
# GD_STOGO
# GN_ISRES
# GN_ESCH
# LN_NEWUOA_BOUND
# LN_BOBYQA
# LN_PRAXIS
# LN_NELDERMEAD
# LN_SBPLX
# LD_SLSQP
# LD_TNEWTON_PRECOND_RESTART
# LD_TNEWTON_RESTART
# LD_TNEWTON_PRECOND
In [24]:
x_init = [0.9, 0.8, 2.5, 2.5, 0.1, 0.1]
x_lb = [0.1, 0.1, 1.5, 1.5, 0.05, 0.1]
x_ub = [1.0, 1.0, 4.0, 4.0, core_flood.initial_water_saturation, 0.25]
opt_alg=:LD_SLSQP
opt1 = Opt(opt_alg, length(x_init)) # choose the algorithm
lower_bounds!(opt1, x_lb)
upper_bounds!(opt1, x_ub)
ftol_rel!(opt1, 1e-15)
ftol_abs!(opt1, 1e-15)
min_objective!(opt1, obj_fun)
(fObjOpt, paramOpt, flag) = optimize(opt1, x_init)
Out[24]:
In [25]:
x_init = paramOpt
vis_error(x_init, exp_data1, core_props, fluids, core_flood)
error_calc(x_init, exp_data1, core_props, fluids, core_flood)
Out[25]:
I'm going to try and use JuMP for history matching of core flooding data. I have used it previously for the history matching of foam flooding, but this is slightly different because there I ony had steady state data and did not need to solve a differential equation. Moreover, I could use automatic differentiation that is not possible to use here.
In [14]:
using JuMP, Ipopt
In [33]:
# weight factors:
w_p = ones(length(exp_data1.dp_exp))
temp_val, ind_max = findmax(exp_data1.dp_exp)
# println(ind_max)
w_p[ind_max-1:ind_max+1] .= 10
w_p[end-10:end] .= 1
w_p[1]=1
w_R = ones(length(exp_data1.R_exp))
w_R[20:25] .= 1
w_R[end-10:end] .= 10
# [krw0, kro0, nw, no, swc, sor]
function my_f(krw0, kro0, nw, no, swc, sor)
f_val = error_calc([krw0, kro0, nw, no, swc, sor], exp_data1, core_props,
fluids, core_flood, w_p = w_p, w_R = w_R)
return f_val
end
function ∇f(g_val, krw0, kro0, nw, no, swc, sor)
eps1 = 1e-4
x = [krw0, kro0, nw, no, swc, sor]
f_val = my_f(krw0, kro0, nw, no, swc, sor)
g_val = zeros(length(x))
for j in eachindex(x)
x2 = copy(x)
x2[j]+=eps1
f_val2 = my_f(x2...)
g_val[j] = (f_val2-f_val)/eps1
end
return g_val
end
# test
grad_x = zeros(6)
my_f(1.0, 0.8, 3, 4, 0.2, 0.2)
# ∇f(1.0, 0.8, 2, 2, 0.1, 0.2)
Out[33]:
In [34]:
model = Model(with_optimizer(Ipopt.Optimizer))
# [krw0, kro0, nw, no, swc, sor]
@variable(model, 0.1 <= kro0 <= 1.0, start = x_init[2])
@variable(model, 0.1 <= krw0 <= 1.0, start = x_init[1])
@variable(model, 0.01 <= sor <= 0.4, start = x_init[6])
@variable(model, 0.01 <= swc <= 0.4, start = x_init[5])
@variable(model, 1.0 <= no <= 4.0, start = x_init[4])
@variable(model, 1.0 <= nw <= 4.0, start = x_init[3])
register(model, :my_f, 6, my_f, ∇f)
@NLobjective(model, Min, my_f(krw0, kro0, nw, no, swc, sor))
# JuMP.register(model::Model, s::Symbol, dimension::Integer, f::Function,
# ∇f::Function, ∇²f::Function)
# my_f(x, y) = (x - 1)^2 + (y - 2)^2
# function ∇f(g, x, y)
# g[1] = 2 * (x - 1)
# g[2] = 2 * (y - 2)
# end
JuMP.optimize!(model)
In [28]:
x_init = [value(krw0), value(kro0), value(nw), value(no), value(swc), value(sor)]
vis_error(x_init, exp_data1, core_props, fluids, core_flood)
error_calc(x_init, exp_data1, core_props, fluids, core_flood)
Out[28]:
In [18]:
x_init
Out[18]: