function [f,state0] = oil_water_flow(r)
% =======================油相体积密度、体积系数、粘度等===============================
density_o_sc = 800; % kg/m3
% 第一种模式，直接输入参考压力、体积系数、压缩系数、粘度
oil_model =1;
if oil_model == 1
pro = 30; % reference pressure, MPa
Boi = 1.1; % gas phase volume factor
co = 1e-3 ;% gas phase compressibility, 1/MPa 
voi = 10; % gas viscosity, cp
cvo = 0; % 气相粘度随压力变化的线性系数，cp/MPa
[Ppr,BO,MUO] = cal_oil_prop(pro,Boi,co,voi,cvo);
end
if oil_model == 2
Ppr = [  0.1013
    2.0946
    4.0878
    6.0811
    8.0743
    10.0676
    12.0608
    14.0540
    16.0473
    18.0405
    20.0338
    22.0270
    24.0203
    26.0135
    28.0068
    30.0000
    ];
% BO = 0.01*ones(size(Ppr,1),1);
BO = [   1.18297
    0.0557041
    0.0278168
    0.0182594
    0.0134665
    0.0106154
    0.00874829
    0.00744907
    0.00650639
    0.0058006
    0.0052587
    0.0048336
    0.00449378
    0.00421751
    0.0039895
    0.00379871
    ];
% MUG = 0.01*ones(size(Ppr,1),1);
MUO = [0.0127683
    0.0130191
    0.0133973
    0.013872
    0.014437
    0.015088
    0.0158182
    0.0166173
    0.0174719
    0.0183671
    0.0192882
    0.0202219
    0.0211575
    0.0220865
    0.0230024
    0.0239008
    ];
end
% ========================水相体积系数、粘度===============================
density_w_sc = 1000; % kg/m3
prw = 30; % reference pressure, MPa
Bwi = 1.000; % water phase volume factor, 无因次
cw = 4.0e-4 ;% water phase compressibility, 1/MPa
vwi = 1; % water viscosity, cp
cvw = 0; % 水相粘度随压力变化的线性系数，cp/MPa
% =============================基质相渗===============================
ifpcow = 0; % 是否考虑毛管力
PRM=[
    0.0000 	0.0000 	1.0000 	3.8916 
0.2000 	0.0000 	1.0000 	3.8916 
0.3160 	0.0002 	0.6784 	0.5796 
0.4350 	0.0004 	0.6215 	0.3724 
0.5620 	0.0010 	0.5456 	0.2425 
0.6140 	0.0020 	0.3939 	0.0608 
0.7020 	0.0280 	0.1399 	0.0372 
0.8120 	0.1721 	0.0515 	0.0137 
0.8750 	0.3395 	0.0297 	0.0104 
0.9060 	0.4395 	0.0 	0.0090 
1.0 	0.4395 	0.0 	0.0090
];
SW=PRM(:,1);KRW=PRM(:,2);KRO=PRM(:,3);PCOW=PRM(:,4);

% =============================裂缝相渗===============================
PRF=[
0.0000 	0.0000 	1.0000 	3.8916 
0.2000 	0.0000 	1.0000 	3.8916 
0.3160 	0.0002 	0.6784 	0.5796 
0.4350 	0.0004 	0.6215 	0.3724 
0.5620 	0.0010 	0.5456 	0.2425 
0.6140 	0.0020 	0.3939 	0.0608 
0.7020 	0.0280 	0.1399 	0.0372 
0.8120 	0.1721 	0.0515 	0.0137 
0.8750 	0.3395 	0.0297 	0.0104 
0.9060 	0.4395 	0.0 	0.0090 
1.0 	0.4395 	0.0 	0.0090
];
SWF = PRF(:,1);KRWF = PRF(:,2);KROF = PRF(:,3);PCOWF=PRF(:,4);
%%
% ===================== 高速非达西流 Forchheimer 方程 ====================
beta_non_Darcy_flow = 0; % 1e-8
% ===================== 启动压力梯度 ====================
p_grad_threshold = 0.00; % MPa/m
%% 结合流体性质参数生成相应的数据体
f = fluidPVT_oil_water_flow(Ppr, BO, MUO, Bwi, prw, cw, vwi, cvw, SW, KRO, KRW, PCOW, SWF, KROF, KRWF, PCOWF, density_o_sc, ifpcow, r.rpt);
f.beta_non_Darcy_flow = beta_non_Darcy_flow;
f.Dwsi = density_w_sc;
f.Dosi = density_o_sc;
f.p_grad_threshold = p_grad_threshold;
%% 初值条件 Initial Condition Section
% =================压力、饱和度初值======================================= %
%压力初值要考虑重力，给出油藏下表面压力值
P = [20 * ones(r.nmc, 1); 20 * ones(r.nfc, 1)];
% plow=25; % P =plow-1e-6*800*9.8*r.z;
Sw = [0.2 * ones(r.nmc, 1); 0.2 * ones(r.nfc, 1)];
state0 = initialRS_oil_water_flow(P, Sw);
end