📄 boastview.cpp
字号:
// BOASTView.cpp : implementation of the CBOASTView class
//
#include "stdafx.h"
#include "BOAST.h"
#include "BOASTDoc.h"
#include "BOASTView.h"
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#define lll 100
#define mmm 100
#define nnn 11
#define maxwell 180
#define maxpvt 12
#define maxword 80
#define NN 20
#define MM 20
#define PI 3.1415926
/* bubble data */
float pbo,vslope,bslope,rslope,pmaxt,rhosco;
float rhoscg,rhoscw,pbot[lll][mmm][nnn];
int msat,mpot,mpwt,mpgt,ireprs;
/* coef data */
float aw[lll][mmm][nnn],ae[lll][mmm][nnn],an[lll][mmm][nnn],as[lll][mmm][nnn];
float ab[lll][mmm][nnn],at[lll][mmm][nnn],e[lll][mmm][nnn],b[lll][mmm][nnn];
/* sarray */
float pn[lll][mmm][nnn],iocode,son[lll][mmm][nnn],swn[lll][mmm][nnn];
/* sarrcy */
float sgn[lll][mmm][nnn],so1[lll][mmm][nnn],sw1[lll][mmm][nnn],sg1[lll][mmm][nnn];
/* sarrby */
float a1[lll][mmm][nnn],a2[lll][mmm][nnn],a3[lll][mmm][nnn],sum[lll][mmm][nnn],gam[lll][mmm][nnn],qs[lll][mmm][nnn];
/* sparm1 */
float kz[lll][mmm][nnn],el[lll][mmm][nnn],tx[lll+1][mmm][nnn],ty[lll][mmm+1][nnn],tz[lll][mmm][nnn+1];
/* sprtps */
/* spvt */
/* srate */
float pid[maxwell][nnn],pwf[maxwell][nnn],pwfc[maxwell][nnn];
float gmo[maxwell][nnn],gmw[maxwell][nnn],gmg[maxwell][nnn],qvo[maxwell];
float qvw[maxwell],qvg[maxwell],qvt[maxwell],cumo[maxwell][nnn],cumw[maxwell][nnn],cumg[maxwell][nnn];
int kip[maxwell],layer[maxwell];
/* ssoln */
float bo[lll][mmm][nnn],bw[lll][mmm][nnn],bg[lll][mmm][nnn];
float qo[lll][mmm][nnn],qw[lll][mmm][nnn],qg[lll][mmm][nnn];
float gowt[lll][mmm][nnn],gwwt[lll][mmm][nnn],ggwt[lll][mmm][nnn];
float ow[lll+1][mmm][nnn],oe[lll+1][mmm][nnn],ww[lll+1][mmm][nnn],we[lll+1][mmm][nnn];
float os[lll][mmm+1][nnn],on[lll][mmm+1][nnn],ws[lll][mmm+1][nnn],wn[lll][mmm+1][nnn];
float ot[lll][mmm][nnn+1],ob[lll][mmm][nnn+1],wt[lll][mmm][nnn+1],wb[lll][mmm][nnn+1];
float qowg[lll][mmm][nnn],vp[lll][mmm][nnn],ct[lll][mmm][nnn];
/* vector */
int iqn1[maxwell],iqn2[maxwell],iqn3[maxwell];
char wellid[maxwell][8];
int pjnum;
float ppj[2000],fwj[2000],timej[2000];
/* codes sub prog is here .solution method, debug print, and time-step control parameters */
static int nn=8000,nvqn;
static float fact1=1.25;
static float fact2=0.75;
static float tmax=9500.;
static float wormax=95.;
static float gormax=10000.;
static float pamin=(float)0.01406/0.0068948;
static float pamax=(float)21.093/0.0068948;
static float dsmax=0.05;
static float dpmax=(float)0.2031/0.0068948;
static FILE *fr,*fw;
static FILE *oko, *okw, *okp;
/******************** Global Parameter *************************/
int Now,ii, jj, kk, nmax;
float Tarry[10000], Parry[10000], Rarry[10000], kx[lll][mmm][nnn],ky[lll][mmm][nnn],varel[lll][mmm];
float dx[lll][mmm][nnn],dy[lll][mmm][nnn],dz[lll][mmm][nnn],phi[lll][mmm][nnn];
float p[lll][mmm][nnn],so[lll][mmm][nnn],sw[lll][mmm][nnn],sg[lll][mmm][nnn];
/************************************************************/
float pavg, pavg0, cwor, cgor, gor, wor;
float opr,wpr,gpr,wir,gir,cop=0.,cwp=0.,cgp=0.;
float cwi=0.,cgi=0.,mbeo,mbew,mbeg;
#ifdef _DEBUG
#define new DEBUG_NEW
#undef THIS_FILE
static char THIS_FILE[] = __FILE__;
#endif
/////////////////////////////////////////////////////////////////////////////
// CBOASTView
IMPLEMENT_DYNCREATE(CBOASTView, CView)
BEGIN_MESSAGE_MAP(CBOASTView, CView)
//{{AFX_MSG_MAP(CBOASTView)
ON_COMMAND(ID_BOAST_CALCULATION, OnBoastCalculation)
//}}AFX_MSG_MAP
// Standard printing commands
ON_COMMAND(ID_FILE_PRINT, CView::OnFilePrint)
ON_COMMAND(ID_FILE_PRINT_DIRECT, CView::OnFilePrint)
ON_COMMAND(ID_FILE_PRINT_PREVIEW, CView::OnFilePrintPreview)
END_MESSAGE_MAP()
/////////////////////////////////////////////////////////////////////////////
// CBOASTView construction/destruction
CBOASTView::CBOASTView()
{
// TODO: add construction code here
}
CBOASTView::~CBOASTView()
{
}
BOOL CBOASTView::PreCreateWindow(CREATESTRUCT& cs)
{
// TODO: Modify the Window class or styles here by modifying
// the CREATESTRUCT cs
return CView::PreCreateWindow(cs);
}
/////////////////////////////////////////////////////////////////////////////
// CBOASTView drawing
void CBOASTView::OnDraw(CDC* pDC)
{
CBOASTDoc* pDoc = GetDocument();
ASSERT_VALID(pDoc);
// TODO: add draw code for native data here
}
/////////////////////////////////////////////////////////////////////////////
// CBOASTView printing
BOOL CBOASTView::OnPreparePrinting(CPrintInfo* pInfo)
{
// default preparation
return DoPreparePrinting(pInfo);
}
void CBOASTView::OnBeginPrinting(CDC* /*pDC*/, CPrintInfo* /*pInfo*/)
{
// TODO: add extra initialization before printing
}
void CBOASTView::OnEndPrinting(CDC* /*pDC*/, CPrintInfo* /*pInfo*/)
{
// TODO: add cleanup after printing
}
/////////////////////////////////////////////////////////////////////////////
// CBOASTView diagnostics
#ifdef _DEBUG
void CBOASTView::AssertValid() const
{
CView::AssertValid();
}
void CBOASTView::Dump(CDumpContext& dc) const
{
CView::Dump(dc);
}
CBOASTDoc* CBOASTView::GetDocument() // non-debug version is inline
{
ASSERT(m_pDocument->IsKindOf(RUNTIME_CLASS(CBOASTDoc)));
return (CBOASTDoc*)m_pDocument;
}
#endif //_DEBUG
/////////////////////////////////////////////////////////////////////////////
// CBOASTView message handlers
void CBOASTView::cal()
{
char words[maxword];
int i,j,k,n;
char ssn[50];
float sat[maxpvt],krot[maxpvt],krwt[maxpvt],krgt[maxpvt],pcowt[maxpvt],pcgot[maxpvt],pot[maxpvt];
float muot[maxpvt],bot[maxpvt],bopt[maxpvt],rsot[maxpvt],rsopt[maxpvt];
float pwt[maxpvt],muwt[maxpvt],bwt[maxpvt],bwpt[maxpvt],rswt[maxpvt],rswpt[maxpvt];
float pgt[maxpvt],mugt[maxpvt],bgt[maxpvt],bgpt[maxpvt],crt[maxpvt];
float mcfgi=0.0;
float mcfg,mcfg1,mcfgt=0.0;
float ooip[nnn],owip[nnn],odgip[nnn],ofgip[nnn];
int im=60,jm=30,km=8;
float eti=0.,ft=0.,ftmax=0.;
float d5615=1.0/5.615;
float d288=1.0/288.;
int iwlcng,ichang,iwlrep;
int isumry,ipmap,isomap,iswmap,isgmap,ipbmap;
float day,dtmin,dtmax;
/* --- add the vary -----*/
float delt,delt0,dsmc=0.0001,dpmc=0.001,div1,resvol;
float scfo,scfw,scfg,scfg1,pp,bpt;
float rso=1.,rsw=1.,ff1,ff2,boder=1.,rsoder=1.,bwder=1.,rswder=1.,bgder=1.;
float co,cw,cg,cr=1.,tooip,towip,todgip,tofgip,stbo,stbw;
float stboi,stbwi;
int itflag;
/* in the cal so -------*/
float ppn,bbo=0.,bbw=0.,bbg=0.,vpp=0.,dp1,dp2,dp3,dp4,dp5,dp6,daodp=0.,dawdp=0.;
float ppm,som,swm,sgm,dpo,dso,dsw,dsg;
int ip,jp,kp;
float rsonew,pbonew=1.;
float flag,tor,tgr,twr,toc,tgc,twc,qoo,qww,qgg;
int ij,iq1,iq2,iq3,lay,nlp;
int numft=10;
pjnum=(int)0;
CDC* pDC=GetDC();
pDC->SelectStockObject(WHITE_BRUSH);
pDC->Rectangle(0,0,640,380);
if((fr=fopen("bos.dat","r"))==NULL)
{AfxMessageBox("error on bostar.dat\n");
return;
}
if((fw=fopen("bos.out","w"))==NULL)
{ printf("error on bostar.dat\n");
exit(1);
}
if((oko=fopen("figoil.out","w"))==NULL)
{ printf("error on figoil.out\n");
exit(1);
}
if((okw=fopen("figwat.out","w"))==NULL)
{ printf("error on figwater.out\n");
exit(1);
}
if((okp=fopen("figprs.out","w"))==NULL)
{ printf("error on figpres.out\n");
exit(1);
}
initvary();
grid1(); /*--------------- dx dy dz ------------ ---*/
parm1(); /*--------------- phi kx,ky,kz---------- --*/
tran1(); /*-caculate interblock transmissibilities--*/
/*-----------read pvt data ----------------*/
table(sat,krot,krwt,krgt,pcowt,pcgot,pot,muot,bot,bopt,rsot,
rsopt,pwt,muwt,bwt,bwpt,rswt,rswpt,pgt,mugt,bgt,bgpt,crt);
uinit1(); /**** establish initial conditions*********/
/* ----- calculate the inittial reservoir volume--------*/
resvol=0.0;
scfo=0.0;
scfw=0.0;
scfg=0.0;
scfg1=0.0;
for(k=1;k<=kk;k++)
{
ooip[k]=0.0;
owip[k]=0.0;
odgip[k]=0.0;
ofgip[k]=0.0;
for(j=1;j<=jj;j++)
for(i=1;i<=ii;i++)
{
pp=p[i][j][k];
bpt=pbot[i][j][k];
vp[i][j][k]=vp[i][j][k]*dx[i][j][k]*dy[i][j][k]*dz[i][j][k];
resvol=resvol+vp[i][j][k];
/**** note: we are assuming initial phi is at initial reservoir pressure */
rso=intpvt(bpt,rslope,pot,rsot,mpot,pp,rso);
rsw=interp(pwt,rswt,mpwt,pp,rsw);
bo[i][j][k]=intpvt(bpt,bslope,pot,bot,mpot,pp,bo[i][j][k]);
bw[i][j][k]=interp(pwt,bwt,mpwt,pp,bw[i][j][k]);
bg[i][j][k]=interp(pgt,bgt,mpgt,pp,bg[i][j][k]);
ff1=so[i][j][k]/bo[i][j][k];
ff2=sw[i][j][k]/bw[i][j][k];
scfo=scfo+vp[i][j][k]*ff1;
scfw=scfw+vp[i][j][k]*ff2;
scfg=scfg+vp[i][j][k]*sg[i][j][k]/bg[i][j][k];
scfg1=scfg1+vp[i][j][k]*(rso*ff1+rsw*ff2);
boder=interp(pot,bopt,mpot,pp,boder);
rsoder=interp(pot,rsopt,mpot,pp,rsoder);
bwder=interp(pwt,bwpt,mpwt,pp,bwder);
rswder=interp(pwt,rswpt,mpwt,pp,rswder);
bgder=interp(pgt,bgpt,mpgt,pp,bgder);
if(pp>pbot[i][j][k]) {boder=bslope;
rsoder=rslope; }
co=-(boder-bg[i][j][k]*rsoder)/bo[i][j][k];
cw=-(bwder-bg[i][j][k]*rswder)/bw[i][j][k];
cg=-bgder/bg[i][j][k];
cr=interp(pgt,crt,mpgt,pp,cr);
ct[i][j][k]=cr+co*so[i][j][k]+cw*sw[i][j][k]+cg*sg[i][j][k];
ooip[k]=ooip[k]+d5615*0.000001*son[i][j][k]*vp[i][j][k]/bo[i][j][k];
owip[k]=owip[k]+d5615*0.000001*swn[i][j][k]*vp[i][j][k]/bw[i][j][k];
odgip[k]=odgip[k]+0.001*0.000001*(rso*son[i][j][k]*vp[i][j][k]/
bo[i][j][k]+rsw*swn[i][j][k]*vp[i][j][k]/bw[i][j][k]);
ofgip[k]=ofgip[k]+0.001*0.000001*sgn[i][j][k]*vp[i][j][k]/bg[i][j][k];
}
fprintf(fw,"%7d oil:%10.2f water: %10.2f solutiongas: %10.2f free gas:%10.2f\n",
k,ooip[k]*rhosco*0.01602*158.988/1000.,owip[k]*158.988*rhoscw*0.01602/1000.,
odgip[k]*0.02832,ofgip[k]*0.02832);
}
tooip=0.0;
towip=0.0;
todgip=0.0;
tofgip=0.0;
for(k=1;k<=kk;k++)
{
tooip=tooip+ooip[k];
towip=towip+owip[k];
todgip=todgip+odgip[k];
tofgip=tofgip+ofgip[k];
}
fprintf(fw,"total: oil:%10.2f water: %10.2f solution gas:%10.2f free gas:%10.2f\n",
tooip*rhosco*158.988*0.01602/1000.,towip*rhoscw*158.988*0.01602/1000.,
todgip*0.02832,tofgip*0.02832);
stbo=scfo*d5615;
stbw=scfw*d5615;
mcfg=scfg*0.001;
mcfg1=scfg1*0.001;
stboi=stbo;
stbwi=stbw;
if(mcfgi<=1.e-7 && mcfgt<=1.e-7)mbeg=0.0;
initprtp();
/*---------------the end of volume -----------------------*/
fscanf(fr,"%s",words); /*------------ begin to read active data -----------------*/
fprintf(fw,"%s\n",words);
nmax=nn+1;
for(n=1;n<=nmax;n++)
{
Now=n;
if(ft<ftmax) goto Mark46;
if(numft==0)numft=1;
fscanf(fr,"%d%d%d%d%d%d%d%d%d",&iwlcng,&ichang,&iwlrep,&isumry,
&ipmap,&isomap,&iswmap,&isgmap,&ipbmap);
fscanf(fr,"%f%f%f%",&day,&dtmin,&dtmax);
if(feof(fr)!=0) goto M1006;
if(iwlcng==1) nodes();
if(iwlcng==-1) goto M1006;
delt=day;
ftmax=eti+ichang*delt;
Mark46:
if(n==1) { delt0=delt; goto Mark1050;}
if(dsmc<dsmax && dpmc<dpmax) delt=delt*fact1;
if(dsmc>dsmax || dpmc>dpmax)delt=delt*fact2;
if(delt<dtmin)delt=dtmin;
if(delt>dtmax)delt=dtmax;
if((eti+delt)>ftmax)delt=ftmax-eti;
Mark1050:
ft=eti+delt;
/*** itflag count the number of time step repitions. */
itflag=0;
/*** reentry point for repeated time step */
Mark1060:
div1=1.0/delt;
Mark105:
/**** esatblish rates & caculate bhfp(if pid is nonzero) */
if(nvqn!=0) qrate(sat,krot,krwt,krgt,pot,muot,bot,rsot,
pwt,muwt,bwt,rswt,pgt,mugt,bgt);
/**** caculate seven diagonal matrix for pressure solution */
solmat(div1,d288,sat,krot,krwt,krgt,pcowt,pcgot,pot,muot,rsot,
pwt,muwt,rswt,pgt,mugt);
/**** modify matrix elements for wells under implicit control. */
if(nvqn!=0) pratei(pot,bot,rsot,pwt,bwt,rswt,pgt,bgt);
/**** caculate new pressure distribution */
lsor();
/**** caculate implicit rates. */
if(nvqn!=0) prateo(pot,rsot,pwt,rswt);
/********* caculate new fluid saturations********/
scfo=0.0;
scfw=0.0;
scfg=0.0;
scfg1=0.0;
resvol=0.0;
for(k=1;k<=kk;k++)
for(j=1;j<=jj;j++)
for(i=1;i<=ii;i++)
{
ppn=pn[i][j][k];
pp=p[i][j][k];
bpt=pbot[i][j][k];
rso=intpvt(bpt,rslope,pot,rsot,mpot,pp,rso);
rsw=interp(pwt,rswt,mpwt,pp,rsw);
cr=interp(pgt,crt,mpgt,ppn,cr);
bpt=pbot[i][j][k];
bbo=intpvt(bpt,bslope,pot,bot,mpot,pp,bbo);
bbw=interp(pwt,bwt,mpwt,pp,bbw);
bbg=interp(pgt,bgt,mpgt,pp,bbg);
vpp=vp[i][j][k]*(1.0+cr*(p[i][j][k]-ppn));
resvol=resvol+vpp;
dp1=0.0;
dp2=0.0;
dp3=0.0;
dp4=0.0;
dp5=0.0;
dp6=0.0;
if((i-1)>0) dp1=p[i-1][j][k]-pp;
if((i+1)<=ii)dp2=p[i+1][j][k]-pp;
if((j-1)>0) dp3=p[i][j-1][k]-pp;
if((j+1)<=jj)dp4=p[i][j+1][k]-pp;
if((k-1)>0) dp5=p[i][j][k-1]-pp;
if((k+1)<=kk)dp6=p[i][j][k+1]-pp;
daodp=ow[i][j][k]*dp1+oe[i][j][k]*dp2+os[i][j][k]*dp3
+on[i][j][k]*dp4+ot[i][j][k]*dp5+ob[i][j][k]*dp6;
dawdp=ww[i][j][k]*dp1+we[i][j][k]*dp2+ws[i][j][k]*dp3
+wn[i][j][k]*dp4+wt[i][j][k]*dp5+wb[i][j][k]*dp6;
sw[i][j][k]=((dawdp+gwwt[i][j][k]-qw[i][j][k])*delt+vp[i][j][k]*
swn[i][j][k]/bw[i][j][k])*bbw/vpp;
so[i][j][k]=((daodp+gowt[i][j][k]-qo[i][j][k])*delt+vp[i][j][k]*
son[i][j][k]/bo[i][j][k])*bbo/vpp;
sg[i][j][k]=1.0-so[i][j][k]-sw[i][j][k];
if(sg[i][j][k]>0.0) goto M405;
sg[i][j][k]=0.0;
so[i][j][k]=1.0-sw[i][j][k];
M405:
vp[i][j][k]=vpp;
bo[i][j][k]=bbo;
bw[i][j][k]=bbw;
bg[i][j][k]=bbg;
ff1=so[i][j][k]/bo[i][j][k];
ff2=sw[i][j][k]/bw[i][j][k];
scfo=scfo+vp[i][j][k]*ff1;
scfw=scfw+vp[i][j][k]*ff2;
scfg=scfg+vp[i][j][k]*sg[i][j][k]/bg[i][j][k];
scfg1=scfg1+vp[i][j][k]*(rso*ff1+rsw*ff2);
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -