notebook.community

Edit and run 





In [11]:



    


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



    


















    






WARNING: replacing module FractionalFlow.
┌ Info: Precompiling LsqFit [2fda8390-95c7-5789-9bda-21331edee243]
└ @ Base loading.jl:1242










    
Out[11]:








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 0x7fe5daa88898>

















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










    
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 0x7fe5daa05908>, PyObject <matplotlib.axes._subplots.AxesSubplot object at 0x7fe5dcfff198>)

















In [5]:



    


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



    


















    

















    
Out[5]:








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

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

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 R_calc(exp_data.t_exp_R), dp_calc(exp_data.t_exp_dp)
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, dp = value_calc(rel_perm_param, exp_data1, core_props, fluids, core_flood)
    return [R; dp]
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]:








([0.0, 0.008646351148232072, 0.0064260485727136735, 0.02178008221754791, 0.008889753417706652, 0.00328637382985146, 0.006381408862702492, 0.005629587593630362, 0.0019165291156212428, 0.0017686791221853662  …  0.160452544751863, 0.1603675150354671, 0.16028290415948443, 0.1601987362875763, 0.16011493393053486, 0.16003149240752199, 0.15994847275163437, 0.1598658502313511, 0.15978360284862386, 0.15970173766610962], [3.0385398682327907, 3.429805043534102, 3.763808734383482, 4.131554264525912, 4.3657008237836905, 4.69335773717772, 4.892060417814959, 5.075108747358546, 5.3013263386155085, 5.473758002272222  …  3.881257491179571, 3.8784675962961668, 3.8756974594095777, 3.8729473372196592, 3.87021555041089, 3.8675018376741734, 3.8648072828195636, 3.8621312441243703, 3.8594730239586883, 3.856832638721718])

















In [17]:



    


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]

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] .= 1
w_p[1]=1
w_R = ones(length(exp_data1.R_exp))
w_R[20:25] .= 1
w_R[end-3:end] .= 10
w = [w_R; w_p]

fit = curve_fit(obj_lsq, [exp_data1.t_exp_R; exp_data1.t_exp_dp], 
    [exp_data1.R_exp; exp_data1.dp_exp], w, x_init, lower = x_lb, upper = x_ub)



    


















    






23










    
Out[17]:








LsqFit.LsqFitResult{Array{Float64,1},Array{Float64,1},Array{Float64,2},Array{Float64,1}}([0.1832917386613496, 0.8288809397320573, 4.0, 1.5, 0.1962529923175783, 0.1], [0.0, 0.0008408267192457075, 0.001054058486717989, 0.0012653660096467734, -0.0015542168313464888, -0.0013469781493674615, 0.0031779181088801567, 0.00035926902442270925, 0.0005664203895845987, 0.0007742877931425363  …  4980.409172658182, 4980.969990442932, 4981.525506327905, 4982.075876366665, 4982.6211112031815, 4983.161320050926, 4983.696558864536, 4984.226925534889, 4984.7524305117995, 4985.273142167621], [0.0 0.0 … 0.0 0.0; -0.002242574509618174 0.0004959040123152895 … 0.0009225772409935119 -0.046187882739342614; … ; -51316.37344061732 -13.008856872708051 … -15.193652452057554 -14.90556042254398; -51316.66156448802 -12.881656574784818 … -15.046038797156099 -14.752018899276612], false, [1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0  …  1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0])

















In [18]:



    


x_init = fit.param
vis_error(x_init, exp_data1, core_props, fluids, core_flood)
error_calc(x_init, exp_data1, core_props, fluids, core_flood)



    


















    

















    

















    
Out[18]:








([0.0, 0.02635505791636104, 0.016925139240102544, 0.013657049958201432, 0.012329283679967877, 0.008610293984659855, 0.01701093480403033, 0.0016313973393230993, 0.002260246358751108, 0.002755669688417556  …  0.184852334849844, 0.18476151679162003, 0.18467160498276597, 0.18458258212829948, 0.18449443749374575, 0.1844071572244876, 0.1843207286529585, 0.18423513840754438, 0.18415037586425612, 0.1840664295226036], [0.06450326434286709, 0.00012409419912773388, 0.05336204590195966, 0.10134481429243704, 0.1374940862378707, 0.17713911284835612, 0.22640082722211577, 0.2521086781457729, 0.28237503008359116, 0.3103542308729185  …  1.122550150176599, 1.1228335328109924, 1.1231142339610527, 1.1233923123240246, 1.1236677898385097, 1.123940710655756, 1.1242111081386539, 1.1244790237965736, 1.1247444796292978, 1.1250075086564255])

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