notebook.community

Edit and run 





In [1]:



    


# load the package
include("../FractionalFlow/FractionalFlow.jl")
using PyPlot, SetPyPlot, Dierckx, Statistics, NLopt
import Calculus
import GR
FF = FractionalFlow



    


















    






┌ Info: Recompiling stale cache file /home/ali/.julia/compiled/v1.2/SetPyPlot/Rkvls.ji for SetPyPlot [d6c70c59-9b85-50b1-926c-19fb5cf24e7d]
└ @ Base loading.jl:1240
┌ Warning: Package SetPyPlot does not have PyPlot in its dependencies:
│ - If you have SetPyPlot checked out for development and have
│   added PyPlot as a dependency but haven't updated your primary
│   environment's manifest file, try `Pkg.resolve()`.
│ - Otherwise you may need to report an issue with SetPyPlot
└ Loading PyPlot into SetPyPlot from project dependency, future warnings for SetPyPlot are suppressed.










    
Out[1]:








Main.FractionalFlow

Water-flooding





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]:








PyObject <matplotlib.legend.Legend object at 0x7f69c7852588>

















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)



    


















    






Progress: 100%|█████████████████████████████████████████| Time: 0:00:08










    
Out[3]:








([0.0, 26.086956521739133, 26.607186199524936, 27.135507511278284, 27.67240634320934, 28.218400516364337, 28.774042266160855, 29.33992094890637, 29.91666599949902, 30.50495016746274  …  5937.3177179224285, 6564.964864736294, 7310.971927867449, 8207.808975006947, 9300.41617919656, 10652.48751343839, 12356.809746482852, 14552.934447722162, 17458.74357803524, 19565.217391304348], [0.0, 0.006666666666666668, 0.006799614250989706, 0.006934629697326673, 0.007071837176597942, 0.0072113690208486635, 0.007353366356907773, 0.0074979797980538505, 0.0076453701998719726, 0.0077957094872404785  …  1.5173145279135094, 1.6777132432103863, 1.868359492677237, 2.0975511825017756, 2.3767730235724542, 2.7223023645453663, 3.1578513796567287, 3.7190832477512195, 4.461678914386783, 5.0], [0.0, 0.008307692307685011, 0.008473365451226038, 0.008641615468969048, 0.008812597096983947, 0.008986475241357903, 0.009163425767831633, 0.009343636363721174, 0.009527307479832943, 0.009714653361015308  …  0.7125107850115477, 0.7155383707111747, 0.7185246684648909, 0.7214619180623228, 0.7243419625168532, 0.7271559553398653, 0.7298940269376648, 0.7325448834514464, 0.7350953087125995, 0.7365900794227], [2242.5, 2305.734685888674, 2307.355536004805, 2308.986310678418, 2310.62636906479, 2312.2751296576316, 2313.9320886554488, 2315.596840519138, 2317.2691006082314, 2318.9487302833586  …  4681.681946972254, 4649.764171730654, 4618.228540283337, 4587.175690185526, 4556.708717095651, 4526.936048760283, 4497.974758699333, 4469.954585332351, 4443.022951688239, 4427.254353064796], [0.0, 0.003, 0.006, 0.009000000000000001, 0.012, 0.015, 0.018000000000000002, 0.021, 0.024, 0.027  …  0.123, 0.126, 0.129, 0.132, 0.135, 0.138, 0.14100000000000001, 0.14400000000000002, 0.147, 0.15], [1.0, 0.7989258599876767, 0.7984393358425793, 0.7979800290450736, 0.797534952877075, 0.7970991749018759, 0.7966702153052446, 0.7962466292295635, 0.7958274918248032, 0.795412170534303  …  0.7831889142012902, 0.7828304209027822, 0.7824730746993335, 0.782116860602225, 0.7817617640725233, 0.7814077709896464, 0.7810548676244633, 0.7807030406155645, 0.780352276944878, 0.7801771593382116], [1.0, 1.0000000000001026, 0.9999999999999544, 0.999999999999509, 0.9999999999992694, 0.9999999999992175, 0.9999999999990471, 0.9999999999988729, 0.9999999999988551, 0.9999999999990851  …  0.9999999999994613, 0.9999999999992298, 0.999999999999132, 0.9999999999992364, 0.9999999999993301, 0.9999999999994139, 0.9999999999993046, 0.9999999999991385, 0.9999999999992722, 0.9999999999993806], [0.0, -1.711739514045568e-108, -1.746569850414833e-108, -1.782673111283491e-108, -1.820133500718822e-108, -1.8590428523785183e-108, -1.8995014747665467e-108, -1.9416191078534836e-108, -1.9855160081927466e-108, -2.0313241826985282e-108  …  0.9999985007861333, 0.9999998363369871, 0.9999999844679713, 0.9999999987771271, 0.9999999999139335, 0.999999999995234, 1.0000000000008638, 1.0000000000007336, 1.000000000000335, 0.9999999999993805])

synthetic experimental data





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]:








(Figure(PyObject <Figure size 800x500 with 2 Axes>), PyObject <matplotlib.axes._subplots.AxesSubplot object at 0x7f69c76eeba8>, PyObject <matplotlib.axes._subplots.AxesSubplot object at 0x7f69c9dbdf60>)

















In [5]:



    


plot(t_exp_dp, dp_exp, "o", t_sec, dp_core)
legend(["Analytical", "Numerical"])



    


















    

















    
Out[5]:








PyObject <matplotlib.legend.Legend object at 0x7f69c36a37b8>

define the objective function





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]:








exp_data(Real[0.0, 98.31767533318768, 196.63535066637536, 294.953025999563, 393.2707013327507, 491.58837666593837, 589.906051999126, 688.2237273323137, 786.5414026655014, 884.8590779986891  …  18680.35831330566, 18778.675988638846, 18876.993663972033, 18975.311339305223, 19073.62901463841, 19171.946689971595, 19270.264365304785, 19368.58204063797, 19466.89971597116, 19565.217391304348], Real[2242.5, 2382.128775673707, 2512.9953521661364, 2643.8423673312104, 2792.5201557768337, 2923.3335362920134, 3054.1548119246854, 3202.8353234859637, 3333.644906448923, 3464.460588605819  …  4436.691912494659, 4436.071639197636, 4435.457548034828, 4434.849537166554, 4434.247520718975, 4433.651413110157, 4433.061124183338, 4432.476569199722, 4431.89766279601, 4431.32432788954], Real[0.0, 98.31767533318768, 196.63535066637536, 294.953025999563, 393.2707013327507, 491.58837666593837, 589.906051999126, 688.2237273323137, 786.5414026655014, 884.8590779986891  …  18680.35831330566, 18778.675988638846, 18876.993663972033, 18975.311339305223, 19073.62901463841, 19171.946689971595, 19270.264365304785, 19368.58204063797, 19466.89971597116, 19565.217391304348], Real[0.0, 0.031903808442163284, 0.062277684795673145, 0.09265295312820346, 0.12605897241797712, 0.15643811370055924, 0.18681619472947614, 0.2202216564676862, 0.2506011733595147, 0.2809798998758703  …  0.7383707061041689, 0.7384285104024199, 0.7384857406425488, 0.7385424062711412, 0.7385985152353743, 0.7386540754616802, 0.7387090953066585, 0.7387635826197505, 0.7388175453119816, 0.7388709905210578])

Calculating error values

The optimizer needs to provide two different time series to the objective function: the first time series is for the recovery factor and the second one is for the pressure drops. There is no parameter that indicates the end of the first time series, but it can be easily inferred from the jump to zero. One other trick is to just supply the data inside the objective function and accept the time series as a dummy input variable.





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 error_R_norm, error_dp_norm
end

function value_calc(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)
    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")
    return dp_calc(exp_data.t_exp_dp), R_calc(exp_data.t_exp_R)
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

function obj_lsq(t, rel_perm_param)
    R_er, dp_er = error_calc(rel_perm_param, exp_data1, core_props, fluids, core_flood)
    return [R_er; dp_er]
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]:








4.576103795950522

















In [ ]:



    


fit = curve_fit(obj_lsq, xdata, ydata, p0)

define the objective function and gradients and weight factors





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



    


















    






23










    
Out[23]:








1.0713534659832842

















In [9]:



    


grad_x



    


















    
Out[9]:








6-element Array{Float64,1}:
  0.24796568682372566
  0.41012431561071594
 -1.6544442529164982 
 -1.665444599678878  
 -1.8706704817250674 
  1.8432281272295903 

















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]:








(0.01409400817359404, [0.3978714543738633, 0.9538120333900926, 2.8772028814937194, 1.8677706182694838, 0.12606550988980245, 0.2034362090288226], :FTOL_REACHED)

















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]:








0.007544830494364934

Content source: simulkade/peteng

Similar notebooks:

notebook.community | gallery | about