|
|
// netgen-2010.cpp : 定义控制台应用程序的入口点。
|
|
|
//
|
|
|
|
|
|
#include "stdafx.h"
|
|
|
#include "singlePhaseSolver.h"
|
|
|
#include "windows.h"
|
|
|
#include <iostream>
|
|
|
#include <fstream>
|
|
|
#include <vector>
|
|
|
#include <string>
|
|
|
#include <sstream>
|
|
|
#include <cstdlib>
|
|
|
#include <algorithm>
|
|
|
#include <stdexcept>
|
|
|
#include "iterator"
|
|
|
using namespace std;
|
|
|
void Ng2VTKWithData(Ng_Mesh *mesh, std::vector<double> vecField, std::string filename)
|
|
|
{
|
|
|
std::cout << "Convert Ng_Mesh to vtk legacy file version of 4.2" << std::endl;
|
|
|
int i, nseg, ne, np, matnum;
|
|
|
int nodes[3];
|
|
|
double point[2];
|
|
|
// int *EleNodeNum;
|
|
|
// double *NodeCoors;
|
|
|
|
|
|
np = Ng_GetNP_2D(mesh);
|
|
|
ne = Ng_GetNE_2D(mesh);
|
|
|
nseg = Ng_GetNSeg_2D(mesh);
|
|
|
|
|
|
std::ofstream vtkFile;
|
|
|
vtkFile.open(filename);
|
|
|
|
|
|
if (!vtkFile.is_open())
|
|
|
{
|
|
|
std::cerr << "Error opening file: " << "mesh.vtk" << std::endl;
|
|
|
return;
|
|
|
}
|
|
|
|
|
|
// Write VTK file header
|
|
|
vtkFile << "# vtk DataFile Version 4.1" << std::endl;
|
|
|
vtkFile << "Mesh data" << std::endl;
|
|
|
vtkFile << "ASCII" << std::endl;
|
|
|
vtkFile << "DATASET UNSTRUCTURED_GRID" << std::endl;
|
|
|
|
|
|
// Write points
|
|
|
vtkFile << "POINTS " << np << " float" << std::endl;
|
|
|
for (i = 1; i <= np; i++)
|
|
|
{
|
|
|
Ng_GetPoint_2D(mesh, i, point);
|
|
|
vtkFile << point[0] << " " << point[1] << " " << 0 << std::endl;
|
|
|
}
|
|
|
|
|
|
// Write cells
|
|
|
vtkFile << "CELLS " << ne << " " << ne * 4 << std::endl;
|
|
|
std::vector<int> vecMater;
|
|
|
for (int i = 1; i <= ne; i++)
|
|
|
{
|
|
|
Ng_GetElement_2D(mesh, i, nodes, &matnum);
|
|
|
vecMater.emplace_back(matnum);
|
|
|
vtkFile << 3 << " " << nodes[0] - 1 << " " << nodes[1] - 1 << " " << nodes[2] - 1 << std::endl;
|
|
|
}
|
|
|
|
|
|
// Write cell types
|
|
|
vtkFile << "CELL_TYPES " << ne << std::endl;
|
|
|
for (int i = 1; i <= ne; i++)
|
|
|
{
|
|
|
vtkFile << 5 << std::endl; // Assuming getCellType is a method that returns the VTK cell type
|
|
|
}
|
|
|
|
|
|
// Write cell data
|
|
|
vtkFile << "CELL_DATA " << ne << std::endl;
|
|
|
vtkFile << "SCALARS Material int 1" << std::endl;
|
|
|
vtkFile << "LOOKUP_TABLE default" << std::endl;
|
|
|
for (int i = 0; i < vecMater.size(); i++)
|
|
|
{
|
|
|
vtkFile << vecMater[i] << std::endl;
|
|
|
}
|
|
|
|
|
|
// Write point data
|
|
|
vtkFile << "POINT_DATA " << np << std::endl;
|
|
|
vtkFile << "SCALARS Pressure double 1" << std::endl;
|
|
|
vtkFile << "LOOKUP_TABLE default" << std::endl;
|
|
|
for (int i = 0; i < vecField.size(); i++)
|
|
|
{
|
|
|
vtkFile << vecField[i] << std::endl;
|
|
|
}
|
|
|
vtkFile.close();
|
|
|
}
|
|
|
|
|
|
int _tmain(int argc, _TCHAR *argv[])
|
|
|
{
|
|
|
/************************************************************************/
|
|
|
/* prepare data for mesh generation */
|
|
|
/************************************************************************/
|
|
|
// 多边形外边界，通过指定四条边来定义
|
|
|
CBoundShape shape;
|
|
|
shape.boundShape = POLYGON;
|
|
|
shape.nNumSegs = 4;
|
|
|
Segment seg1;
|
|
|
seg1.p1 = Point(-1500, -1500);
|
|
|
seg1.p2 = Point(1500, -1500);
|
|
|
Segment seg2;
|
|
|
seg2.p1 = Point(1500, -1500);
|
|
|
seg2.p2 = Point(1500, 1500);
|
|
|
Segment seg3;
|
|
|
seg3.p1 = Point(1500, 1500);
|
|
|
seg3.p2 = Point(-1500, 1500);
|
|
|
Segment seg4;
|
|
|
seg4.p1 = Point(-1500, 1500);
|
|
|
seg4.p2 = Point(-1500, -1500);
|
|
|
|
|
|
shape.vecSegments.push_back(seg1);
|
|
|
shape.vecSegments.push_back(seg2);
|
|
|
shape.vecSegments.push_back(seg3);
|
|
|
shape.vecSegments.push_back(seg4);
|
|
|
|
|
|
// shape.cCenter = Point(0, 0);
|
|
|
// shape.dRadius = 1000;
|
|
|
|
|
|
CBoundWell w1, w2;
|
|
|
// 第一口井为直井（默认为直井）
|
|
|
w1.cCenter = Point(0, 0); // 井中心坐标(0,0)
|
|
|
w1.dRadius = 0.09144; // 井半径0.09144
|
|
|
|
|
|
// 第二口井为垂直裂缝井
|
|
|
w2.wellType = VERTICAL_FRACTURED;
|
|
|
w2.dHf = 20; // 裂缝半长20米
|
|
|
w2.dTheata = 30; // 裂缝角度30°
|
|
|
w2.cCenter = Point(50, 0); // 井的中心位置（50，0）
|
|
|
w2.dRadius = 0.1; // 井的半径0.1
|
|
|
|
|
|
// 第三口井为多段压裂水平井
|
|
|
CBoundWell w3;
|
|
|
w3.wellType = MULTI_FRACED_HORIZONTAL;
|
|
|
w3.cCenter = Point(142, -226); // 中心点位置
|
|
|
w3.dRadius = 0.1; // 半径
|
|
|
w3.dBeta = 20; // 水平井角度
|
|
|
w3.dLength = 200; // 水平井长度
|
|
|
w3.nFracNum = 10; // 裂缝数量
|
|
|
w3.dHf = 20; // 裂缝半长
|
|
|
|
|
|
//!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
|
|
|
// 在该参数设置条件下，一次只能添加一口井进行计算求解
|
|
|
//!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
|
|
|
std::vector<CBoundWell> wells;
|
|
|
// wells.push_back(w1);
|
|
|
// wells.push_back(w2);
|
|
|
wells.push_back(w3);
|
|
|
|
|
|
CBoundPolygon limit;
|
|
|
limit.boundShape = POLYGON;
|
|
|
limit.nNumSegs = 5;
|
|
|
// 通过指定的边来定义多边形
|
|
|
limit.vecSegments.push_back(Segment(Point(-516.2335347838829, -58.37577524138442), Point(-592.72401802758, -180.9965861656542)));
|
|
|
limit.vecSegments.push_back(Segment(Point(-592.72401802758, -180.9965861656542), Point(-485.12090993338666, -280.2195626173259)));
|
|
|
limit.vecSegments.push_back(Segment(Point(-485.12090993338666, -280.2195626173259), Point(-328.9162066043042, -229.04373034669038)));
|
|
|
limit.vecSegments.push_back(Segment(Point(-328.9162066043042, -229.04373034669038), Point(-376.36137535544276, -90.33591813134171)));
|
|
|
limit.vecSegments.push_back(Segment(Point(-376.36137535544276, -90.33591813134171), Point(-516.2335347838829, -58.37577524138442)));
|
|
|
limit.bAnchor = true;
|
|
|
|
|
|
// 通过指定的点来定义多边形，逆时针排列
|
|
|
std::vector<Point> points;
|
|
|
points.emplace_back(Point(-9.23023, 31.1986));
|
|
|
points.emplace_back(Point(-27.0886, 11.26));
|
|
|
points.emplace_back(Point(-23.5281, -14.9567));
|
|
|
points.emplace_back(Point(10.0862, -27.2422));
|
|
|
points.emplace_back(Point(39.6129, -1.9348));
|
|
|
points.emplace_back(Point(32.122, 15.5807));
|
|
|
points.emplace_back(Point(22.0438, 31.371));
|
|
|
CBoundPolygon limit1(points);
|
|
|
|
|
|
std::vector<Point> points2;
|
|
|
points2.emplace_back(Point(250.074, -101.492));
|
|
|
points2.emplace_back(Point(-44, -204.034));
|
|
|
points2.emplace_back(Point(7, -336.098));
|
|
|
points2.emplace_back(Point(312.333, -234.724));
|
|
|
CBoundPolygon limit2(points2);
|
|
|
|
|
|
std::vector<CBoundPolygon> limits;
|
|
|
limits.push_back(limit);
|
|
|
limits.push_back(limit1);
|
|
|
limits.push_back(limit2);
|
|
|
|
|
|
std::vector<Frac> faults;
|
|
|
faults.push_back(Frac(Point(-267.3, 120), Point(316.7, 150)));
|
|
|
faults.push_back(Frac(Point(316.7, 150), Point(316.7, -50)));
|
|
|
|
|
|
std::vector<Frac> fracs;
|
|
|
fracs.push_back(Frac(Point(-326.159, 460.059), Point(165.83631, 459.085)));
|
|
|
fracs.push_back(Frac(Point(50, -350), Point(358.840, -289.811)));
|
|
|
|
|
|
/************************************************************************/
|
|
|
/* load netgen mesh generation library and generate mesh */
|
|
|
/************************************************************************/
|
|
|
HMODULE hMod = LoadLibrary("meshgene.dll");
|
|
|
if (NULL == hMod)
|
|
|
{
|
|
|
std::cout << "meshgene.dll加载失败！\n";
|
|
|
return 1;
|
|
|
}
|
|
|
|
|
|
typedef Ng_Mesh *(*NgMeshGen2DByIn2d)(std::vector<CBoundWell> &wells, const std::vector<CBoundPolygon> limits, const std::vector<Frac> &faults, const std::vector<Frac> &fracs, const CBoundShape &shape);
|
|
|
typedef Ng_Mesh *(*NgMeshGen2D)(std::string strFileIn2D);
|
|
|
typedef void *(*Ng2VTK)(Ng_Mesh *mesh);
|
|
|
// Get the address of the Ng_Mesh function
|
|
|
NgMeshGen2DByIn2d meshgene1 = (NgMeshGen2DByIn2d)GetProcAddress(hMod, "NgMeshGen2DByIn2d");
|
|
|
NgMeshGen2D meshgene2 = (NgMeshGen2D)GetProcAddress(hMod, "NgMeshGen2D");
|
|
|
if (NULL == meshgene1)
|
|
|
{
|
|
|
FreeLibrary(hMod);
|
|
|
std::cout << "meshgene文件加载函数地址获取失败！\n";
|
|
|
return 1;
|
|
|
}
|
|
|
Ng_Mesh *mesh = meshgene1(wells, limits, faults, fracs, shape);
|
|
|
|
|
|
Ng2VTK ng = (Ng2VTK)GetProcAddress(hMod, "Ng2VTK");
|
|
|
ng(mesh);
|
|
|
|
|
|
/************************************************************************/
|
|
|
/* solve the equation by calling the singlephasesolver */
|
|
|
/************************************************************************/
|
|
|
// Prepare the reservoir base data
|
|
|
double pBaseData[8];
|
|
|
pBaseData[0] = 30.47;
|
|
|
pBaseData[1] = 15.1454;
|
|
|
pBaseData[2] = 30.48;
|
|
|
pBaseData[3] = 1.5;
|
|
|
pBaseData[4] = 1;
|
|
|
pBaseData[5] = 0.23;
|
|
|
pBaseData[6] = 0.00043511;
|
|
|
pBaseData[7] = 1;
|
|
|
|
|
|
// prepare the well data
|
|
|
// well w1
|
|
|
// w1.dRc=0.01;
|
|
|
wells[0].dSkin = 0;
|
|
|
wells[0].dC = 0.0232056;
|
|
|
|
|
|
wells[0].bisMultFlow = true;
|
|
|
wells[0].nTimeNumQ = 5;
|
|
|
wells[0].pdTimeQ[0] = 12;
|
|
|
wells[0].pdQ[0] = 158.987;
|
|
|
wells[0].pdTimeQ[1] = 12;
|
|
|
wells[0].pdQ[1] = 190.785;
|
|
|
wells[0].pdTimeQ[2] = 12;
|
|
|
wells[0].pdQ[2] = 222.582;
|
|
|
wells[0].pdTimeQ[3] = 48;
|
|
|
wells[0].pdQ[3] = 238.481;
|
|
|
wells[0].pdTimeQ[4] = 72;
|
|
|
wells[0].pdQ[4] = 0;
|
|
|
|
|
|
// 裂缝井必须指定裂缝的导流能力，通过wFrac变量设置
|
|
|
wells[0].wFrac[0].dFc = 1000;
|
|
|
wells[0].wFrac[1].dFc = 1000;
|
|
|
|
|
|
// 多段压裂水平井必须指定裂缝的导流能力1000
|
|
|
// 导流能力的设置需要在网格划分后完成
|
|
|
for (int i = 0; i < wells[0].vecFracs.size(); ++i)
|
|
|
{
|
|
|
wells[0].vecFracs[i].dFc = 1000; //
|
|
|
}
|
|
|
|
|
|
// prepare the data for limits
|
|
|
limits[0].dComKr = 1;
|
|
|
limits[0].dComW = 1;
|
|
|
|
|
|
HMODULE hMod_solver = LoadLibrary("singlePhaseSolverDll.dll");
|
|
|
if (NULL == hMod_solver)
|
|
|
{
|
|
|
std::cout << "singlePhaseSolverDll.dll加载失败！\n";
|
|
|
return 1;
|
|
|
}
|
|
|
TimeDis timeDis;
|
|
|
timeDis.dTB = 1.0e-5; //
|
|
|
|
|
|
typedef int (*Solver)(const TimeDis time, double *m_pBaseData, nglib::Ng_Mesh *mesh, std::vector<CBoundWell> &m_wWell, const std::vector<CBoundPolygon> &comPolygon, const std::vector<Frac> &faultPolygon, std::vector<Frac> &fracs, CBoundShape &outBound, std::vector<std::pair<double, std::vector<double>>> &vecNodePre);
|
|
|
Solver solverFun = (Solver)GetProcAddress(hMod_solver, "singlePhaseSolverNgmesh");
|
|
|
|
|
|
if (NULL == solverFun)
|
|
|
{
|
|
|
FreeLibrary(hMod_solver);
|
|
|
std::cout << "singlePhaseSolverNgmesh failed!\n";
|
|
|
return 1;
|
|
|
}
|
|
|
// std::vector<std::vector<Point>> vecBP;
|
|
|
std::vector<std::pair<double, std::vector<double>>> vecFieldPre;
|
|
|
solverFun(timeDis, pBaseData, mesh, wells, limits, faults, fracs, shape, vecFieldPre);
|
|
|
|
|
|
// export bottom hole pressure history
|
|
|
std::cout << "Export log-log plot and history plot data to files\n";
|
|
|
std::ofstream historyFile;
|
|
|
historyFile.precision(8);
|
|
|
historyFile.open("output/history.txt");
|
|
|
if (!historyFile.is_open())
|
|
|
{
|
|
|
std::cerr << "Error opening file: " << "history.txt" << std::endl;
|
|
|
return 1;
|
|
|
}
|
|
|
int nWell = 0;
|
|
|
for (int i = 0; i < wells[nWell].m_PressureData.size(); i++)
|
|
|
{
|
|
|
historyFile << wells[nWell].m_PressureData[i].x << "\t" << wells[nWell].m_PressureData[i].y << std::endl;
|
|
|
}
|
|
|
historyFile.close();
|
|
|
|
|
|
// calculate the pressure derivative
|
|
|
typedef bool (*PreLog)(const std::vector<Point> &pressure, const int &nSection, double *TimeQ, double *dQ, int nTimeQ, std::vector<Point> &logPre);
|
|
|
PreLog preLogFun = (PreLog)GetProcAddress(hMod_solver, "logLogPre");
|
|
|
|
|
|
if (NULL == preLogFun)
|
|
|
{
|
|
|
FreeLibrary(hMod_solver);
|
|
|
std::cout << "preLogFun failed!\n";
|
|
|
return 1;
|
|
|
}
|
|
|
|
|
|
std::vector<Point> logPre;
|
|
|
preLogFun(wells[0].m_PressureData, 5, wells[0].pdTimeQ, wells[0].pdQ, wells[0].nTimeNumQ, logPre);
|
|
|
|
|
|
// export the pressure derivative
|
|
|
std::ofstream logFile;
|
|
|
logFile.precision(8);
|
|
|
logFile.open("output/logPre.txt");
|
|
|
if (!logFile.is_open())
|
|
|
{
|
|
|
std::cerr << "Error opening file: " << "logPre.txt" << std::endl;
|
|
|
return 1;
|
|
|
}
|
|
|
for (int i = 0; i < logPre.size(); i++)
|
|
|
{
|
|
|
logFile << logPre[i].x << "\t" << logPre[i].y << "\t" << logPre[i].z << "\t" << std::endl;
|
|
|
}
|
|
|
logFile.close();
|
|
|
|
|
|
// export with pressure field with vtk format (version 4.2)
|
|
|
for (int i = 0; i < vecFieldPre.size(); ++i)
|
|
|
{
|
|
|
if (i % 10 == 0)
|
|
|
{
|
|
|
Ng2VTKWithData(mesh, vecFieldPre[i].second, "output/" + std::to_string(long double(i)) + ".vtk");
|
|
|
}
|
|
|
}
|
|
|
|
|
|
FreeLibrary(hMod);
|
|
|
FreeLibrary(hMod_solver);
|
|
|
}
|