function [ r, Times, OutputRs, Wellpara, trun ] = main1()
%% ==================================PART 1: 模型选取========================================
% 1-- classical EDFM
modelflag = 1;% 此处固定为1即可
%% ==================================PART 2: 基质网格定义=====================================
dx=[50*ones(1,20) ];nx=size(dx,2);
dy=[50*ones(1,10) ];ny=size(dy,2);
dz=[10*ones(1,1)];nz=size(dz,2);
NTG= 1*ones(nx*ny*nz,1);
r = GridProp_pre(dx,dy,dz,nx,ny,nz,NTG);
%% ==================================PART 3: 裂缝分布数据输入=====================================
% 四种输入方式: 1表示工程应用输入；2表示向量输入；3表示工程应用输入，包括基准点、倾角、方位角、抬升角；4表示从.fab文件读取
input_style = 1;
%% 输入方式1：工程应用输入
%仅能刻画具有双对称性质的矩形缝或椭圆缝
% 基准点坐标（1），方位角（2），倾角（3），抬升角（4），缝长（椭圆长轴长）（5），缝高（椭圆短轴长）（6），类型1是矩形缝、2是椭圆缝（7）
if input_style==1
input_content={
                     [255,255,5],90,90,0,200,10,1;
                     [305,255,5],90,90,0,200,10,1;
                     [355,255,5],90,90,0,200,10,1;
                     [405,255,5],90,90,0,200,10,1;
                     [455,255,5],90,90,0,200,10,1;
                     [505,255,5],90,90,0,200,10,1;
                     [555,255,5],90,90,0,200,10,1;
                     [605,255,5],90,90,0,200,10,1;
                     [655,255,5],90,90,0,200,10,1;
                     [705,255,5],90,90,0,200,10,1;
                     [755,255,5],90,90,0,200,10,1;
                   };
[f,fellip] = sort_fracture(input_content);
flowBarrierFlags = [];
m=size(f,1);n=size(fellip,1);nf=(m+n)/5;
frac_information = fractureInformation_input_engineering_vector(f,fellip,flowBarrierFlags);
end
%% 输入方式2：向量输入
%首先对矩形缝进行研究，裂缝参数包括三个向量和两个参数的取值范围，如果是其它类型缝，则是两个参数间的函数关系
%因此一个5行三列的矩阵可以确定一条裂缝
%向量分量可为非整数，以此保证参数范围的取值为整数即可，
%故可先随意确定为整数的参数范围，再根据缝长缝高去确定向量分量的取值
% 矩形缝
if input_style==2
 f=[
     0,1200,0;10,-20,0;0,0,10;0,60,0;0,3,0;
     0,1200,0;10,20,0;0,0,10;0,15,0;0,3,0;
     150,1500,0;20,-10,0;0,0,10;0,20,0;0,3,0;
     550,1300,0;10,10,0;0,0,10;0,10,0;0,3,0;
     650,1400,0;10,-20,0;0,0,10;0,55,0;0,3,0;
     605,705,0;10,0,0;0,0,10;-20,20,0;0,3,0;
 ];
% 椭圆缝
% fellip=[320,210,10;10,2,0;0,0,2;0,0,8;0,0,5;
%            280,300,10;10,2,0;0,0,2;0,0,10;0,0,4;];
%     300,205,10;10,0,0;0,0,2;0,0,10;0,0,5];
fellip=[];
m=size(f,1);n=size(fellip,1);nf=(m+n)/5;%裂缝条数
frac_information = fractureInformation_input_engineering_vector(f,fellip);
end
%% 输入方式3：二维输入
% load('fractures_generated.mat');
% fractureLines = fractures_generated;
if input_style==3
fractureLines = [
    105,150;105,350;
    205,150;205,350;
    305,150;305,350;
    405,150;405,350;
    505,150;505,350;
    605,150;605,350;
    705,150;705,350;
    805,150;805,350;
    905,150;905,350;
    ];
fellip=[]; % 椭圆缝
fractureHeights = [
    10;10;10;10;10;10;10;10;10];
% 椭圆缝
fellip=[];
flowBarrierFlags = [];% flowBarrierFlags = [1;2;3;4];
[f,frac_information] = input_fracture_2D(fractureLines,fractureHeights,flowBarrierFlags);
m=size(f,1);n=size(fellip,1);nf=(m+n)/5;%裂缝条数
end
%% 输入方式4：从.fab文件读取 Frac_PT
if input_style==4

end

%% ==================================PART 4: 基质网格及裂缝物性参数输入=====================================
% 1-- 单重介质  
% 2-- 双重介质（双孔双渗）
grid_model = 2;
if grid_model == 1
r = grid_discretization_SP_model(modelflag, r, f, frac_information, fellip, nf);
end
if grid_model == 2
r = grid_discretization_DP_model(modelflag, r, f, frac_information, fellip, nf);
end

%% 生成向量化变成所需要的算子
os = OperatorRS(r.N, r.nex, r.nc);

%% ==================================PART 5: 模型类型选取及初始状态设定=========================================
% 1-- 气水两相流  
% 2-- 油水两相流
flow_model = 3;
if flow_model == 1
[f,state0] = gas_water_flow(r);
end
if flow_model == 2
[f,state0] = oil_water_flow(r);
end
if flow_model == 3
[f,state0] = multi_component_flow(r);
end

%% ==================================PART 6: 井制度schedule设定================================================
%将直井、水平井与多段压裂水平井分开处理
%% 直井、常规水平井定义
% Wellcpara: wellname(1) nperf，射孔点数量(2) index(3) rw(4) skin(5) welltype(6) 
% welltype=1 直井，welltype=2 水平井沿x方向，welltype=3 水平井沿y方向 welltype=4 多段压裂水平井 
% 直井、水平井，该种井型的处理是将射孔段安排在基质网格
% 对于直井、斜井、水平井，采取与eclipse一致的方式，给出其在基质网格中位置
% 每个射孔点由下px,py,pz确定，分别表示该射孔点在x方向，y方向的网格编号，及层数即z反方向上的网格编号（最初是以垂直向上建立的基质网格）
well1={
           % 'w1', 2, [20 10 1;20 10 2;], 0.178/2, 0, 1;
           % 'w2', 2, [80 10 1;80 10 2;], 0.178/2, 0, 1;
           % 'w3', 1, [20 40 1;], 0.178/2, 0, 1;
           % 'w4', 1, [80 40 1;], 0.178/2, 0, 1;
%            'w2', 3, [67 20 1;67 20 2;67 20 3], 0.178/2, 0, 1;
           };
well1 = handle_well1(well1,r);
%% 压裂水平井定义
% 多段压裂水平井：射孔段设置在裂缝单元上，可以给出该射孔点所在的坐标，然后去寻找该点所在的裂缝单元编号
num_fracture_wells = 1;
welloc = cell(num_fracture_wells,1);
perfnum = cell(num_fracture_wells,1);
welloc{1,1}=[
    255,255,5;
    305,255,5;
    355,255,5;
    405,255,5;
    455,255,5;
    505,255,5;
    555,255,5;
    605,255,5;
    655,255,5;
    705,255,5;
    755,255,5;
             ];
for i = 1:num_fracture_wells
perfnum{i,1} = findWelloc(r, welloc{i,1});%在用以Wellc0之中
end
% perfnum2 = findWelloc(r, welloc2);
well2={
    'w1', length(perfnum{1,1}), perfnum{1,1}, 0.178/2, 0,4;
           };
Wellc = handle_well1_well2(well1,well2,welloc,r); %加井和射孔位置
%% 井制度设置（开关井、定压/定产、便于吞吐等）
number_phases = 3;
well_schedules = cell(number_phases,1);
time = zeros(number_phases,1);
% 不同方式对应的收敛难度、需要不一样，需分阶段调控，
% 例如在注入阶段需要小一些，如果需要做早期试井，dtmax、dtmin都需要小一些
dtmax = zeros(number_phases,1);
dtmin = zeros(number_phases,1);

% (1)well state: 'open'表示该井是开的；'close'表明该井一直是关的，如果只是阶段性关井，就用'open'定流量0生产
% (2)well type: 'pro'表示是生产井；'inj'表示是注入井;
% (3)protype：'const_q'是定流量；'const_pwf'则是定压
% (4)value：定的流量值或者井底压力值
% (5)value：此代码该值与（3）值相同
% (8)value：注入表活剂的浓度 kg/m3
% (10)value：注入水/压裂液的盐浓度（矿化度），一般是 1kg/m3
time(1) = 10; dtmax(1) = 0.1; dtmin(1) = 0.001;
well_schedules{1,1} = {
                       % 'w1','open','inj','const_pwf',30,30;
                       % 'w2','open','inj','const_pwf',30,30;
                       % 'w3','open','inj','const_pwf',30,30;
                       % 'w4','open','inj','const_pwf',30,30;
                       'w1','open','inj','const_pwf',40,40,'Cs_inj',0.5,'Cb_inj',1};
time(2) = 30; dtmax(2) = 0.5; dtmin(2) = 0.001;
well_schedules{2,1} = {
                       'w1','open','pro','const_q',0,0;
                       % 'w1','open','pro','const_pwf',10,10;
                       % 'w2','open','pro','const_q',0,0;
                       };
time(3) = 200; dtmax(3) = 5; dtmin(3) = 0.001;
well_schedules{3,1} = {
                       'w1','open','pro','const_pwf',10,10;
                       % 'w2','open','pro','const_pwf',10,10;
                       };

% time(1) = 200;
% well_schedules{1,1} = {
%                        'w1','open','pro','const_pwf',10,10;
%                        % 'w2','open','pro','const_pwf',10,10;
%                        };

%% ======================================PART 7: Solver设置===========================================
yitap = 5; % 50-500psi
yitas = 0.04; % 0.05-0.5
omega = 0.5; % 0-1
Nmax = 50;
epsave = 1e-6;
epsmax = 1e-6;

%% ======================================PART 8: 正式计算===========================================
[Times, OutputRs, Wellpara, trun] = solver_NR(r, flow_model, f, os, state0, yitap,yitas,omega,Nmax,epsave,epsmax,dtmin,dtmax,Wellc,time,well_schedules);
run_time=toc;
trun.run_time=run_time;
end