#include"pch.h"
void Write2DVectorToCSV(const std::vector<std::vector<double>>& data, const std::string& filename) {
	std::ofstream file(filename.c_str());  // VS2010需使用.c_str()
	if (file.is_open()) {
		for (size_t row = 0; row < data.size(); ++row) {
			for (size_t col = 0; col < data[row].size(); ++col) {
				// 设置固定小数格式和精度
				file << std::fixed << std::setprecision(15) << data[row][col];
				// 非最后一列时添加逗号
				if (col != data[row].size() - 1) {
					file << ",";
				}
			}
			file << "\n";  // 换行符
		}
		file.close();
	}
}
void Write1DVectorToCSV(const std::vector<double>& data, const std::string& filename, int precision = 15) {
	std::ofstream file(filename.c_str());
	if (file.is_open()) {
		file << std::fixed << std::setprecision(precision);
		for (size_t i = 0; i < data.size(); ++i) {
			file << data[i] << "\n";  // 每个元素单独一行
		}
		file.close();
	}
}

int main()
{

	HX_NWTM_GRID_INPUT p0;
	HX_NWTM_GRID_OUTPUT1 p1;
	HX_NWTM_GRID_OUTPUT2 p2;
	

	//不同井型算例
	int welltype;//1为一口直井，2为一口压裂直井，3为一口多段压裂水平井，4为五口井(含直井，压裂直井，多段压裂水平井，断层)，5为50口直井
	welltype = 1;
	//模型算例
	int flowtype;//1为油单相常数pvt（一口井），2为油单相常数pvt（五口井），3为油单相常数pvt（五十口井），4为油单相变化pvt（一口井），5为水单相常数pvt（一口井），6为水单相变化pvt（一口井），7为气单相变化pvt（一口井），8为气单相拟压力（一口井），9为油水两相（一口井）
	flowtype = 1;
	//非均质性
	int feijunzhi = 0;//0为不考虑储层非均质，1为考虑储层非均质

	//不同井型设置
	dVec1 a(3), b(4), c(5);	
	if (welltype == 1) {
		//一口直井
		p0.Boundary.resize(4);
		b[0] = -1500.0; b[1] = -1500.0; b[2] = -1500.0; b[3] = 1500.0;  p0.Boundary[0] = b;
		b[0] = -1500.0; b[1] = 1500.0;  b[2] = 1500.0;  b[3] = 1500.0;  p0.Boundary[1] = b;
		b[0] = 1500.0;  b[1] = 1500.0;  b[2] = 1500.0;  b[3] = -1500.0; p0.Boundary[2] = b;
		b[0] = 1500.0;  b[1] = -1500.0; b[2] = -1500.0; b[3] = -1500.0; p0.Boundary[3] = b;
		p0.VerticalWell.resize(1);
		a[0] = 0;    a[1] = 0;    a[2] = 0.1;  p0.VerticalWell[0] = a;
		p0.HorizontalWell.resize(0);
		p0.FractureVerticalWell.resize(0);
		p0.MultistageFracturedHorizontalWell.resize(0);
		p0.InclinedWell.resize(0);
		p0.Fault.resize(0);
	}
	else if (welltype == 2) {
		//一口压裂直井
		p0.Boundary.resize(4);
		b[0] = -1500.0; b[1] = -1500.0; b[2] = -1500.0; b[3] = 1500.0;  p0.Boundary[0] = b;
		b[0] = -1500.0; b[1] = 1500.0;  b[2] = 1500.0;  b[3] = 1500.0;  p0.Boundary[1] = b;
		b[0] = 1500.0;  b[1] = 1500.0;  b[2] = 1500.0;  b[3] = -1500.0; p0.Boundary[2] = b;
		b[0] = 1500.0;  b[1] = -1500.0; b[2] = -1500.0; b[3] = -1500.0; p0.Boundary[3] = b;
		p0.VerticalWell.resize(0);
		p0.HorizontalWell.resize(0);
		p0.FractureVerticalWell.resize(1);
		c[0] = -200; c[1] = 0; c[2] = 200; c[3] = 0; c[4] = 0.05; p0.FractureVerticalWell[0] = c;
		p0.MultistageFracturedHorizontalWell.resize(0);
		p0.InclinedWell.resize(0);
		p0.Fault.resize(0);
	}
	else if (welltype == 3) {
		//一口多段压裂水平井
		p0.Boundary.resize(4);
		b[0] = -1500.0; b[1] = -1500.0; b[2] = -1500.0; b[3] = 1500.0;  p0.Boundary[0] = b;
		b[0] = -1500.0; b[1] = 1500.0;  b[2] = 1500.0;  b[3] = 1500.0;  p0.Boundary[1] = b;
		b[0] = 1500.0;  b[1] = 1500.0;  b[2] = 1500.0;  b[3] = -1500.0; p0.Boundary[2] = b;
		b[0] = 1500.0;  b[1] = -1500.0; b[2] = -1500.0; b[3] = -1500.0; p0.Boundary[3] = b;
		p0.VerticalWell.resize(0);
		p0.HorizontalWell.resize(0);
		p0.FractureVerticalWell.resize(0);
		p0.MultistageFracturedHorizontalWell.resize(1);
		p0.MultistageFracturedHorizontalWell[0].resize(7, dVec1(5));
		c[0] = -600; c[1] = -200; c[2] = -600; c[3] = 200; c[4] = 0.05; p0.MultistageFracturedHorizontalWell[0][0] = c;
		c[0] = -400; c[1] = -200; c[2] = -400; c[3] = 200; c[4] = 0.05; p0.MultistageFracturedHorizontalWell[0][1] = c;
		c[0] = -200; c[1] = -200; c[2] = -200; c[3] = 200; c[4] = 0.05; p0.MultistageFracturedHorizontalWell[0][2] = c;
		c[0] = 0; c[1] = -200; c[2] =0; c[3] = 200; c[4] = 0.05; p0.MultistageFracturedHorizontalWell[0][3] = c;
		c[0] = 200; c[1] = -200; c[2] = 200; c[3] = 200; c[4] = 0.05; p0.MultistageFracturedHorizontalWell[0][4] = c;
		c[0] = 400; c[1] = -200; c[2] = 400; c[3] = 200; c[4] = 0.05; p0.MultistageFracturedHorizontalWell[0][5] = c;
		c[0] = 600; c[1] = -200; c[2] = 600; c[3] = 200; c[4] = 0.05; p0.MultistageFracturedHorizontalWell[0][6] = c;
		p0.InclinedWell.resize(0);
		p0.Fault.resize(0);
	}
	else if (welltype == 4) {
		//五口井(含直井，压裂直井，多段压裂水平井，断层)
		p0.Boundary.resize(4);
		b[0] = -1500.0; b[1] = -1500.0; b[2] = -1500.0; b[3] = 1500.0;  p0.Boundary[0] = b;
		b[0] = -1500.0; b[1] = 1500.0;  b[2] = 1500.0;  b[3] = 1500.0;  p0.Boundary[1] = b;
		b[0] = 1500.0;  b[1] = 1500.0;  b[2] = 1500.0;  b[3] = -1500.0; p0.Boundary[2] = b;
		b[0] = 1500.0;  b[1] = -1500.0; b[2] = -1500.0; b[3] = -1500.0; p0.Boundary[3] = b;
		p0.VerticalWell.resize(3);
		a[0] = 0;    a[1] = 0;    a[2] = 0.1;  p0.VerticalWell[0] = a;
		a[0] = 1000; a[1] = 1000; a[2] = 0.1;  p0.VerticalWell[1] = a;
		a[0] = -1000; a[1] = -1000; a[2] = 0.1;  p0.VerticalWell[2] = a;
		p0.HorizontalWell.resize(0);
		p0.FractureVerticalWell.resize(1);
		c[0] = -200; c[1] = -200; c[2] = 200; c[3] = -200; c[4] = 0.05; p0.FractureVerticalWell[0] = c;
		p0.MultistageFracturedHorizontalWell.resize(1);
		p0.MultistageFracturedHorizontalWell[0].resize(3, dVec1(5));
		c[0] = -600; c[1] = 600; c[2] = -400; c[3] = 600; c[4] = 0.1; p0.MultistageFracturedHorizontalWell[0][0] = c;
		c[0] = -600; c[1] = 400; c[2] = -400; c[3] = 400; c[4] = 0.1; p0.MultistageFracturedHorizontalWell[0][1] = c;
		c[0] = -600; c[1] = 200; c[2] = -400; c[3] = 200; c[4] = 0.1; p0.MultistageFracturedHorizontalWell[0][2] = c;
		p0.InclinedWell.resize(0);
		p0.Fault.resize(1);
		c[0] = -500; c[1] = 1000; c[2] = 500; c[3] = 500; p0.Fault[0] = c;
	}
	else if (welltype == 5) {
		//50口直井
		p0.Boundary.resize(4);
		b[0] = -1000.0; b[1] = -1000.0; b[2] = -1000.0; b[3] = 1000.0;  p0.Boundary[0] = b;
		b[0] = -1000.0; b[1] = 1000.0;  b[2] = 1000.0;  b[3] = 1000.0;  p0.Boundary[1] = b;
		b[0] = 1000.0;  b[1] = 1000.0;  b[2] = 1000.0;  b[3] = -1000.0; p0.Boundary[2] = b;
		b[0] = 1000.0;  b[1] = -1000.0; b[2] = -1000.0; b[3] = -1000.0; p0.Boundary[3] = b;
		p0.VerticalWell.resize(50);
		a[0] = 766.32207856101854;    a[1] = -408.4253100821785;    a[2] = 0.108;  p0.VerticalWell[0] = a;
		a[0] = -359.11086885711029;   a[1] = -242.17185988435722;   a[2] = 0.108;  p0.VerticalWell[1] = a;
		a[0] = -396.28113134064483;   a[1] = 618.18974488936828;    a[2] = 0.108;  p0.VerticalWell[2] = a;
		a[0] = -844.38929572547772;   a[1] = 859.25724767620932;    a[2] = 0.108;  p0.VerticalWell[3] = a;
		a[0] = -631.69279373636357;   a[1] = 620.26791301684102;    a[2] = 0.108;  p0.VerticalWell[4] = a;
		a[0] = -858.84439780240791;   a[1] = 333.48071142559911;    a[2] = 0.108;  p0.VerticalWell[5] = a;
		a[0] = -844.38929572547772;   a[1] = 609.87707237947734;    a[2] = 0.108;  p0.VerticalWell[6] = a;
		a[0] = -639.95285206603785;   a[1] = 867.56992018610026;    a[2] = 0.108;  p0.VerticalWell[7] = a;
		a[0] = 764.25706397860017;    a[1] = -190.2176566975379;    a[2] = 0.108;  p0.VerticalWell[8] = a;
		a[0] = -410.73623341757479;   a[1] = 867.56992018610026;    a[2] = 0.108;  p0.VerticalWell[9] = a;
		a[0] = -121.63419187897284;   a[1] = 98.647713021176514;    a[2] = 0.108;  p0.VerticalWell[10] = a;
		a[0] = -152.6094106152517;    a[1] = 364.65323333769061;    a[2] = 0.108;  p0.VerticalWell[11] = a;
		a[0] = -152.6094106152517;    a[1] = 620.26791301684102;    a[2] = 0.108;  p0.VerticalWell[12] = a;
		a[0] = -361.17588343952877;   a[1] = 102.80404927612199;    a[2] = 0.108;  p0.VerticalWell[13] = a;
		a[0] = -844.38929572547772;   a[1] = -231.78101924699342;   a[2] = 0.108;  p0.VerticalWell[14] = a;
		a[0] = -156.73943978008879;   a[1] = 875.88259269599121;    a[2] = 0.108;  p0.VerticalWell[15] = a;
		a[0] = 138.55764550576896;    a[1] = 890.42976958830059;    a[2] = 0.108;  p0.VerticalWell[16] = a;
		a[0] = 446.24481828613784;    a[1] = -439.59783199426988;   a[2] = 0.108;  p0.VerticalWell[17] = a;
		a[0] = -846.45431030789632;   a[1] = 61.240686726666581;    a[2] = 0.108;  p0.VerticalWell[18] = a;
		a[0] = -142.28433770315883;   a[1] = -223.46834673710214;   a[2] = 0.108;  p0.VerticalWell[19] = a;
		a[0] = -375.63098551645885;   a[1] = -472.84852203383423;   a[2] = 0.108;  p0.VerticalWell[20] = a;
		a[0] = -602.78258958250342;   a[1] = 354.2623927003267;     a[2] = 0.108;  p0.VerticalWell[21] = a;
		a[0] = 431.78971620920788;    a[1] = 624.42424927178672;    a[2] = 0.108;  p0.VerticalWell[22] = a;
		a[0] = -834.06422281338484;   a[1] = -474.92669016130685;   a[2] = 0.108;  p0.VerticalWell[23] = a;
		a[0] = -836.12923739580344;   a[1] = -747.16671486023938;   a[2] = 0.108;  p0.VerticalWell[24] = a;
		a[0] = 813.81741395664631;    a[1] = -693.13434354594744;   a[2] = 0.108;  p0.VerticalWell[25] = a;
		a[0] = 458.63490578064966;    a[1] = -718.07236107562062;   a[2] = 0.108;  p0.VerticalWell[26] = a;
		a[0] = 409.07455580260353;    a[1] = 873.80442456851847;    a[2] = 0.108;  p0.VerticalWell[27] = a;
		a[0] = 124.10254342883854;    a[1] = -454.14500888657915;   a[2] = 0.108;  p0.VerticalWell[28] = a;
		a[0] = 138.55764550576896;    a[1] = 96.569544893703778;    a[2] = 0.108;  p0.VerticalWell[29] = a;
		a[0] = 101.38738302223419;    a[1] = 611.95524050694985;    a[2] = 0.108;  p0.VerticalWell[30] = a;
		a[0] = -373.56597093404037;   a[1] = 360.49689708274514;    a[2] = 0.108;  p0.VerticalWell[31] = a;
		a[0] = -604.8476041649219;    a[1] = 79.944199873921661;    a[2] = 0.108;  p0.VerticalWell[32] = a;
		a[0] = -611.0426479121777;    a[1] = -219.31201048215667;   a[2] = 0.108;  p0.VerticalWell[33] = a;
		a[0] = 448.30983286855667;    a[1] = -202.68666546237455;   a[2] = 0.108;  p0.VerticalWell[34] = a;
		a[0] = 735.34685982474002;    a[1] = 880.03892895093713;    a[2] = 0.108;  p0.VerticalWell[35] = a;
		a[0] = 107.58242676948976;    a[1] = -724.30686545803894;   a[2] = 0.108;  p0.VerticalWell[36] = a;
		a[0] = 150.94773300028032;    a[1] = -196.45216107995623;   a[2] = 0.108;  p0.VerticalWell[37] = a;
		a[0] = 762.19204939618135;    a[1] = 113.19488991348589;    a[2] = 0.108;  p0.VerticalWell[38] = a;
		a[0] = -125.76422104381015;   a[1] = -470.77035390636138;   a[2] = 0.108;  p0.VerticalWell[39] = a;
		a[0] = 452.43986203339387;    a[1] = 366.73140146516357;    a[2] = 0.108;  p0.VerticalWell[40] = a;
		a[0] = -350.85081052743578;   a[1] = -751.32305111518485;   a[2] = 0.108;  p0.VerticalWell[41] = a;
		a[0] = 122.03752884641995;    a[1] = 354.2623927003267;     a[2] = 0.108;  p0.VerticalWell[42] = a;
		a[0] = -608.9776333297591;    a[1] = -447.91050450416094;   a[2] = 0.108;  p0.VerticalWell[43] = a;
		a[0] = 407.0095412201847;     a[1] = 113.19488991348589;    a[2] = 0.108;  p0.VerticalWell[44] = a;
		a[0] = 731.21683065990283;    a[1] = 636.89325803662314;    a[2] = 0.108;  p0.VerticalWell[45] = a;
		a[0] = -121.63419187897284;   a[1] = -724.30686545803894;   a[2] = 0.108;  p0.VerticalWell[46] = a;
		a[0] = 289.3037100223255;     a[1] = -925.88917382289742;   a[2] = 0.108;  p0.VerticalWell[47] = a;
		a[0] = 737.41187440715862;    a[1] = 366.73140146516357;    a[2] = 0.108;  p0.VerticalWell[48] = a;
		a[0] = -611.0426479121777;    a[1] = -730.54136984045726;   a[2] = 0.108;  p0.VerticalWell[49] = a;
		p0.HorizontalWell.resize(0);
		p0.FractureVerticalWell.resize(0);
		p0.MultistageFracturedHorizontalWell.resize(0);
		p0.InclinedWell.resize(0);
		p0.Fault.resize(0);
	}
	p0.GridControl = 150.0;
	p0.D = 2;

	//网格计算
	HX_NWTM_GRID(p1, p2, p0, "HX_license.dat");
	HX_NWTM_MODEL_INPUT p3(p2);
	HX_NWTM_MODEL_OUTPUT p4;
	
	//模型设置
	if (flowtype == 1) {
		//油单相常数pvt（一口井）	
		p3.T = 1;
		p3.Rate.t.resize(1);
		p3.Rate.qo.resize(1);
		p3.Rate.t[0].resize(2);  p3.Rate.t[0][0] = 2000;    p3.Rate.t[0][1] = 500;
		p3.Rate.qo[0].resize(2); p3.Rate.qo[0][0] = 10;      p3.Rate.qo[0][1] = 0;
		p3.CS.C.resize(1);
		p3.CS.C[0] = 0.1;
		p3.CS.S.resize(1);
		p3.CS.S[0] = 0.1;
		p3.PVT.p = dVec1(200, 0); for (int i = 0; i < 200; ++i) { p3.PVT.p[i] = (i + 1.0); }
		p3.PVT.Bo = dVec1(200, 1.2);//所有数为一个值
		p3.PVT.miuo = dVec1(200, 0.5);//所有数为一个值
		p3.Base.Pi = 40.0;
		p3.Base.Cti = 1e-3;
		p3.Base.k = dVec1(p2.Trinodexy.size(), 0.001);
		p3.Base.phi = dVec1(p2.Trinodexy.size(), 0.1);
		p3.Base.h = dVec1(p2.Trinodexy.size(), 10);
		p3.Base.d = 1.05;
		p3.Base.dt_Min = 0.0025;
		p3.Base.dt_Max = 12.5;
	}
	else if (flowtype == 2) {
		//油单相常数pvt（五口井）	
		p3.T = 1;
		p3.Rate.t.resize(5);
		p3.Rate.qo.resize(5);
		p3.Rate.qw.resize(5);
		p3.Rate.qg.resize(5);
		p3.Rate.t[0].resize(2);  p3.Rate.t[0][0] = 2000;    p3.Rate.t[0][1] = 500;
		p3.Rate.qo[0].resize(2); p3.Rate.qo[0][0] = 10;      p3.Rate.qo[0][1] = 0;
		p3.Rate.t[1].resize(0);
		p3.Rate.qo[1].resize(0);
		p3.Rate.t[2].resize(3);  p3.Rate.t[2][0] = 1000;    p3.Rate.t[2][1] = 1000;   p3.Rate.t[2][2] = 500;
		p3.Rate.qo[2].resize(3); p3.Rate.qo[2][0] = 30;      p3.Rate.qo[2][1] = 40;    p3.Rate.qo[2][2] = 20;
		p3.Rate.t[3].resize(2);  p3.Rate.t[3][0] = 1500;     p3.Rate.t[3][1] = 1000;
		p3.Rate.qo[3].resize(2); p3.Rate.qo[3][0] = 30;      p3.Rate.qo[3][1] = 20;
		p3.Rate.t[4].resize(2);  p3.Rate.t[4][0] = 1000;    p3.Rate.t[4][1] = 1500;
		p3.Rate.qo[4].resize(2); p3.Rate.qo[4][0] = -50;     p3.Rate.qo[4][1] = -60;
		p3.Pressure.t.resize(0);
		p3.Pressure.p.resize(0);
		p3.CS.C.resize(5);
		p3.CS.C[0] = 0.1; p3.CS.C[1] = 0.1; p3.CS.C[2] = 0.1; p3.CS.C[3] = 0.1; p3.CS.C[4] = 0.1;
		p3.CS.S.resize(5);
		p3.CS.S[0] = 0.1; p3.CS.S[1] = 0.1; p3.CS.S[2] = 0.1; p3.CS.S[3] = 0.1; p3.CS.S[4] = 0.1;
		p3.PVT.p = dVec1(200, 0); for (int i = 0; i < 200; ++i) { p3.PVT.p[i] = (i + 1.0); }
		p3.PVT.Bo = dVec1(200, 1.2);//所有数为一个值
		p3.PVT.miuo = dVec1(200, 0.5);//所有数为一个值
		p3.Base.Pi = 40.0;
		p3.Base.Cti = 1e-3;
		p3.Base.k = dVec1(p2.Trinodexy.size(), 0.001);
		p3.Base.phi = dVec1(p2.Trinodexy.size(), 0.1);
		p3.Base.h = dVec1(p2.Trinodexy.size(), 10);
		p3.Base.d = 1.05;
		p3.Base.dt_Min = 0.0025;
		p3.Base.dt_Max = 12.5;
	}
	else if (flowtype == 3) {
		//油单相常数pvt（五十口井）	
		p3.T = 1;
		dVec1 t;
		t.push_back(12);
		t.push_back(12);
		t.push_back(12);
		t.push_back(48);
		t.push_back(72);
		p3.Rate.t.assign(50, t);
		dVec1 q;
		q.push_back(158.98699999999999);
		q.push_back(190.785);
		q.push_back(222.58199999999999);
		q.push_back(238.48099999999999);
		q.push_back(0.0);
		p3.Rate.qo.assign(50, q);
		p3.CS.C = dVec1(50, 0.01);
		p3.CS.S = dVec1(50, 0);


		p3.PVT.p = dVec1(200, 0); for (int i = 0; i < 200; ++i) { p3.PVT.p[i] = (i + 1.0); }
		p3.PVT.Bo = dVec1(200, 1.07);//所有数为一个值
		p3.PVT.miuo = dVec1(200, 0.79);//所有数为一个值
		p3.Base.Pi = 40;
		p3.Base.Cti = 0.43e-3;
		p3.Base.k = dVec1(p2.Trinodexy.size(), 0.025);
		p3.Base.phi = dVec1(p2.Trinodexy.size(), 0.1);
		p3.Base.h = dVec1(p2.Trinodexy.size(), 10);
		p3.Base.d = 1.05;
		p3.Base.dt_Min = 0.0025;
		p3.Base.dt_Max = 12.5;
	}
	else if (flowtype == 4) {
		//油单相变化pvt（一口井）	
		p3.T = 2;
		p3.Rate.t.resize(1);
		p3.Rate.qo.resize(1);
		p3.Rate.t[0].resize(2);  p3.Rate.t[0][0] = 2000;    p3.Rate.t[0][1] = 500;
		p3.Rate.qo[0].resize(2); p3.Rate.qo[0][0] = 10;      p3.Rate.qo[0][1] = 0;
		p3.CS.C.resize(1);
		p3.CS.C[0] = 0.1;
		p3.CS.S.resize(1);
		p3.CS.S[0] = 0.1;
		p3.PVT.p = dVec1(200, 0); for (int i = 0; i < 200; ++i) { p3.PVT.p[i] = (i + 1.0); }
		p3.PVT.Bo = dVec1(200, 1.2);//数值随压力变化
		p3.PVT.miuo = dVec1(200, 0.5);//数值随压力变化
		p3.PVT.Co = dVec1(200, 0.001);//数值随压力变化
		p3.Base.Pi = 40.0;
		p3.Base.Cf = 1e-3;
		p3.Base.k = dVec1(p2.Trinodexy.size(), 0.001);
		p3.Base.phi = dVec1(p2.Trinodexy.size(), 0.1);
		p3.Base.h = dVec1(p2.Trinodexy.size(), 10);
		p3.Base.d = 1.05;
		p3.Base.dt_Min = 0.0025;
		p3.Base.dt_Max = 12.5;
	}
	else if (flowtype == 5) {
		//水单相常数pvt（一口井）	
		p3.T = 3;
		p3.Rate.t.resize(1);
		p3.Rate.qw.resize(1);
		p3.Rate.t[0].resize(2);  p3.Rate.t[0][0] = 2000;    p3.Rate.t[0][1] = 500;
		p3.Rate.qw[0].resize(2); p3.Rate.qw[0][0] = -10;      p3.Rate.qw[0][1] = 0;
		p3.CS.C.resize(1);
		p3.CS.C[0] = 0.1;
		p3.CS.S.resize(1);
		p3.CS.S[0] = 0.1;
		p3.PVT.p = dVec1(200, 0); for (int i = 0; i < 200; ++i) { p3.PVT.p[i] = (i + 1.0); }
		p3.PVT.Bw = dVec1(200, 1.05);//所有数为一个值
		p3.PVT.miuw = dVec1(200, 0.8);//所有数为一个值
		p3.Base.Pi = 40.0;
		p3.Base.Cti = 1e-3;
		p3.Base.k = dVec1(p2.Trinodexy.size(), 0.001);
		p3.Base.phi = dVec1(p2.Trinodexy.size(), 0.1);
		p3.Base.h = dVec1(p2.Trinodexy.size(), 10);
		p3.Base.d = 1.05;
		p3.Base.dt_Min = 0.0025;
		p3.Base.dt_Max = 12.5;
	}
	else if (flowtype == 6) {
		//水单相变化pvt（一口井）	
		p3.T = 4;
		p3.Rate.t.resize(1);
		p3.Rate.qw.resize(1);
		p3.Rate.t[0].resize(2);  p3.Rate.t[0][0] = 2000;    p3.Rate.t[0][1] = 500;
		p3.Rate.qw[0].resize(2); p3.Rate.qw[0][0] = -10;      p3.Rate.qw[0][1] = 0;
		p3.CS.C.resize(1);
		p3.CS.C[0] = 0.1;
		p3.CS.S.resize(1);
		p3.CS.S[0] = 0.1;
		p3.PVT.p = dVec1(200, 0); for (int i = 0; i < 200; ++i) { p3.PVT.p[i] = (i + 1.0); }
		p3.PVT.Bw = dVec1(200, 1.05);//数值随压力变化
		p3.PVT.miuw = dVec1(200, 0.8);//数值随压力变化
		p3.PVT.Cw = dVec1(200, 0.0001);//数值随压力变化
		p3.Base.Pi = 40.0;
		p3.Base.Cf = 1e-3;
		p3.Base.k = dVec1(p2.Trinodexy.size(), 0.001);
		p3.Base.phi = dVec1(p2.Trinodexy.size(), 0.1);
		p3.Base.h = dVec1(p2.Trinodexy.size(), 10);
		p3.Base.d = 1.05;
		p3.Base.dt_Min = 0.0025;
		p3.Base.dt_Max = 12.5;
	}
	else if (flowtype == 7) {
		//气单相变化pvt（一口井）	
		p3.T = 5;
		p3.Rate.t.resize(1);
		p3.Rate.qg.resize(1);
		p3.Rate.t[0].resize(2);  p3.Rate.t[0][0] = 2000;    p3.Rate.t[0][1] = 500;
		p3.Rate.qg[0].resize(2); p3.Rate.qg[0][0] = 50000;      p3.Rate.qg[0][1] = 0;
		p3.CS.C.resize(1);
		p3.CS.C[0] = 0.1;
		p3.CS.S.resize(1);
		p3.CS.S[0] = 0.1;
		p3.PVT.p = dVec1(200, 0); for (int i = 0; i < 200; ++i) { p3.PVT.p[i] = (i + 1.0); }
		p3.PVT.Bg = dVec1(200, 5e-3);//数值随压力变化
		p3.PVT.miug = dVec1(200, 2e-2);//数值随压力变化
		p3.PVT.Cg = dVec1(200, 2e-2);//数值随压力变化
		p3.Base.Pi = 40.0;
		p3.Base.Cf = 1e-3;
		p3.Base.k = dVec1(p2.Trinodexy.size(), 0.001);
		p3.Base.phi = dVec1(p2.Trinodexy.size(), 0.1);
		p3.Base.h = dVec1(p2.Trinodexy.size(), 10);
		p3.Base.d = 1.05;
		p3.Base.dt_Min = 0.0025;
		p3.Base.dt_Max = 12.5;
	}
	else if (flowtype == 8) {
		//气单相拟压力（一口井）	
		p3.T = 6;
		p3.Rate.t.resize(1);
		p3.Rate.qg.resize(1);
		p3.Rate.t[0].resize(2);  p3.Rate.t[0][0] = 2000;    p3.Rate.t[0][1] = 500;
		p3.Rate.qg[0].resize(2); p3.Rate.qg[0][0] = 50000;      p3.Rate.qg[0][1] = 0;
		p3.CS.C.resize(1);
		p3.CS.C[0] = 0.1;
		p3.CS.S.resize(1);
		p3.CS.S[0] = 0.1;
		p3.PVT.p = dVec1(200, 0); for (int i = 0; i < 200; ++i) { p3.PVT.p[i] = (i + 1.0); }
		p3.PVT.Bg = dVec1(200, 5e-3);//数值随压力变化
		p3.PVT.miug = dVec1(200, 2e-2);//数值随压力变化
		p3.PVT.Cg = dVec1(200, 2e-2);//数值随压力变化
		p3.Base.Pi = 40.0;
		p3.Base.Cf = 1e-3;
		p3.Base.k = dVec1(p2.Trinodexy.size(), 0.001);
		p3.Base.phi = dVec1(p2.Trinodexy.size(), 0.1);
		p3.Base.h = dVec1(p2.Trinodexy.size(), 10);
		p3.Base.d = 1.05;
		p3.Base.dt_Min = 0.0025;
		p3.Base.dt_Max = 12.5;
	}
	else if (flowtype == 9) {
		//油水两相（一口井）	
		p3.T = 8;
		p3.Rate.t.resize(1);
		p3.Rate.qg.resize(1);
		p3.Rate.t[0].resize(2);  p3.Rate.t[0][0] = 2000;    p3.Rate.t[0][1] = 500;
		//定产油量
		p3.Rate.qo[0].resize(2); p3.Rate.qo[0][0] = 10;      p3.Rate.qo[0][1] = 0;
		p3.Rate.qw[0].resize(2); p3.Rate.qw[0][0] = 0;      p3.Rate.qw[0][1] = 0;
		////定产液量
		//p3.Rate.qo[0].resize(2); p3.Rate.qo[0][0] = 10;      p3.Rate.qo[0][1] = 0;
		//p3.Rate.qw[0].resize(2); p3.Rate.qw[0][0] = 10;      p3.Rate.qw[0][1] = 0;
		////定注水量
		//p3.Rate.qo[0].resize(2); p3.Rate.qo[0][0] = 0;      p3.Rate.qo[0][1] = 0;
		//p3.Rate.qw[0].resize(2); p3.Rate.qw[0][0] = -10;      p3.Rate.qw[0][1] = 0;

		p3.CS.C.resize(1);
		p3.CS.C[0] = 0.1;
		p3.CS.S.resize(1);
		p3.CS.S[0] = 0.1;
		p3.PVT.p = dVec1(200, 0); for (int i = 0; i < 200; ++i) { p3.PVT.p[i] = (i + 1.0); }
		p3.PVT.Bo = dVec1(200, 1.2);//数值随压力变化
		p3.PVT.miuo = dVec1(200, 0.5);//数值随压力变化
		p3.PVT.Bw = dVec1(200, 1.05);//数值随压力变化
		p3.PVT.miuw = dVec1(200, 0.8);//数值随压力变化
		p3.PVT.So = dVec1(81, 0); for (int i = 0; i < 81; ++i) { p3.PVT.So[i] = (0.1 + i * 0.01); }
		p3.PVT.Kro = dVec1(81, 0); for (int i = 0; i < 81; ++i) { p3.PVT.Kro[i] = (i*0.0125); }
		p3.PVT.Krw = dVec1(81, 0); for (int i = 0; i < 81; ++i) { p3.PVT.Krw[i] = (1 - i * 0.0125); }
		p3.Base.Pi = 40.0;
		p3.Base.Cf = 1e-3;
		p3.Base.Swi = 0.2;
		p3.Base.k = dVec1(p2.Trinodexy.size(), 0.001);
		p3.Base.phi = dVec1(p2.Trinodexy.size(), 0.1);
		p3.Base.h = dVec1(p2.Trinodexy.size(), 10);
		p3.Base.d = 1.05;
		p3.Base.dt_Min = 0.0025;
		p3.Base.dt_Max = 12.5;
	}

	//非均质性设置
	if (feijunzhi == 1) {
		dVec2 k;
		a[0] = -1000; a[1] = -800; a[2] = 0.001; k.push_back(a);
		a[0] = -800; a[1] = 1000; a[2] = 0.01; k.push_back(a);
		a[0] = 0; a[1] = 100; a[2] = 0.005; k.push_back(a);
		a[0] = 500; a[1] = -1000; a[2] = 0.008; k.push_back(a);
		a[0] = 1000; a[1] = 800; a[2] = 0.02; k.push_back(a);
		HX_KRING_INPUT k0(0.01, 100, 1000, 0, p2.Trinodexy, k);
		HX_KRING_OUTPUT k1;
		HX_NWTM_KRINGING(k1, k0, "HX_license.dat");
		p3.Base.k = k1.v;

		dVec2 phi;
		a[0] = -600; a[1] = 0; a[2] = 0.1; phi.push_back(a);
		a[0] = 0; a[1] = 100; a[2] = 0.2; phi.push_back(a);
		a[0] = 1000; a[1] = -200; a[2] = 0.05; phi.push_back(a);
		HX_KRING_INPUT phi0(0.01, 100, 1000, 0, p2.Trinodexy, phi);
		HX_KRING_OUTPUT phi1;
		HX_NWTM_KRINGING(phi1, phi0, "HX_license.dat");
		p3.Base.phi = phi1.v;

		dVec2 h;
		a[0] = -1000; a[1] = 100; a[2] = 10; h.push_back(a);
		a[0] = 0; a[1] = 0; a[2] = 5; h.push_back(a);
		a[0] = 1000; a[1] = 300; a[2] = 12; h.push_back(a);
		HX_KRING_INPUT h0(0.01, 100, 1000, 0, p2.Trinodexy, h);
		HX_KRING_OUTPUT h1;
		HX_NWTM_KRINGING(h1, h0, "HX_license.dat");
		p3.Base.h = h1.v;

	/*	std::ofstream file("k_phi_h.csv");
		file << std::fixed << std::setprecision(15);
		for (int i = 0; i < p2.Trinodexy.size(); ++i)
		{
			file << p2.Trinodexy[i][0] << "," << p2.Trinodexy[i][1] << "," << p3.Base.k[i] << "," << p3.Base.phi[i] << "," << p3.Base.h[i] <<"\n";
		}

		file.close();*/
	}
	
	//模型计算
	HX_NWTM_MODEL(p4, p3, "HX_license.dat");


	//数据导出
	Write1DVectorToCSV(p4.t, "t.csv");
	Write2DVectorToCSV(p4.pw, "pw.csv");
	Write2DVectorToCSV(p4.p, "C2.csv");
	if (p3.T == 8) {
		Write2DVectorToCSV(p4.So, "So.csv");
	}
		
	return 0;
}