/* SOSimoleo -- Oil reservoir simulator. */ /* Last edited on 2003-09-15 17:28:31 by ra014852 */ /* Receives as input a reservoir data file consisting of: -> Size of each reservoir side (its dimensions); -> Water property functions descriptors {a string + coeff double array}; -> Oil property functions descriptors {a string + coeff double array}; -> General property functions descriptors {a string + coeff double array}; -> Well descriptors {IN/OUT well def. + flux + position double array}. */ #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #define PHASES 2 /* The number of fluids of the system: oil and water. */ #define DAYSECONDS 86400 /* Number of seconds of a day. */ #define NUMPROP 4 /* Number of phase properties. */ #define NUMGENPROP 2 /* Number of general reservoir properties. */ #define TOTPROP PHASES*NUMPROP+NUMGENPROP /* Total number of properties. */ #define NUMMAIN 4 /* Number of approximated main functions. */ #define WATER 0 #define OIL 1 #define MAXDIM 3 #define PLOTDIVS 80 #define STEPDIV 10 #define MAXDEPTHTREE 20 #define MIN_RELPERM 0.1 #define MAX_RELPERM 0.9 #define SOReservoir_FileFormat "2004-02-10" typedef struct Options { char *basisName; /* Prefix for basis filenames (minus extension). */ char *matName; /* Prefix for matrix filenames (default {basisName}). */ char *reservoirName; /* Prefix of the reservoir input file (minus extension). */ char *outName; /* Prefix for output filenames. */ nat maxIter; /* Maximum iterations. */ double tFin; /* Total simulation time (in seconds). */ double tStep; /* Time interval between global iterations (in seconds). */ double relTol; /* Max change in {a} coeffs between iterations. */ double absTol; /* Max function difference between iterations. */ nat gaussOrder; /* Gaussian quadrature order for dot products. */ /* Method used to solve the linear system: */ bool gaussElim; /* TRUE uses the generic Gauss elimination method. */ bool gaussSeidel; /* TRUE uses Gauss-Seidel iteration. */ bool conjGrad; /* TRUE uses the conjugate gradient method. */ double omega; /* Relaxation factor for Gauss-Seidel. */ /* Plotting options: */ bool plotFinal; /* TRUE to plot final solution and error. */ bool plotAll; /* TRUE to plot every iteration and error. */ PlotOptions plt; /* Plot format and style options. */ } Options; typedef struct Resv_Func { char *name; /* Name of the reservoir function. */ char *descr; /* Descriptor of the reservoir function. */ int degree; /* For polynomial approximated functions, its degree. */ double *coeff; /* For polynomial approximated functions, its coefficients. */ SOFunction *F; /* Pressure, Saturation or their differentials. */ FuncMap *FMap; /* The mapping function. May be null. */ int plot; /* Binary verification for plotting: **1: 1D, *1*: 2D, 1**: 3D. */ double fRange; /* Plotting range, from -fRange to fRange. */ double fStep; /* Plotting step. */ int changrng; /* Changes automatically fRange and fStep. */ int ps_eps; /* O: Ps, 1: Eps plot file format. */ int plot_times; /* Number of images to be plotted. */ SOPlot_PSStream *ps; /* Plotting stream descriptor. */ } Resv_Func; typedef struct Resv_Well { int type; /* Water Influx (0) or Oil/Water Outflux (1). */ double *pos; /* Point coordenates of the well. */ double flux; /* Well flux. */ } Resv_Well; typedef struct Resv_Func_vec { nat nel; Resv_Func *el; } Resv_Func_vec; typedef struct Resv_Well_vec { nat nel; Resv_Well *el; } Resv_Well_vec; Options GetOptions(int argn, char **argc); Resv_Func_vec Resv_Func_vec_new(nat nel); Resv_Well_vec Resv_Well_vec_new(nat nel); void ReservoirDataRead ( FILE *rd, SOGrid_Dim pDim, double_vec *res_sizes, double *initP, double *finP, double *init_press, double *init_sat, int *plot_skip, Resv_Func_vec *properties, int *numwells, Resv_Well_vec *wells ); void ResvFuncRead( FILE *rd, Resv_Func *property, char *propname ); void ResvWellRead( FILE *rd, Resv_Well *well, int resdim ); void SetPropFunctions ( Resv_Func_vec *properties, SOFunction *pressure, SOFunction *saturation, SOFunction *dpressure, SOFunction *dsaturation ); void OpenPlot ( Resv_Func *RF, Options *o, int *plotTotCount ); void TotalOilVol ( SOFunction *saturation, FuncMap *NoMap, SOFunction *porosity, FuncMap *FMpor, SOGrid_Tree *tree, SOGrid_Dim pDim, double *tot_vol, double *water_vol, double *oil_vol, double *reserv_box ); void PlotState ( Resv_Func *RF, SOGrid_Dim pDim, SOGrid_Tree *tree, double *initP, double *finP, Options *o, int *plotMeshDepth, char *tag, double time, int iter, int outiter, double *fMax, int *plot_count, int *plotTotCount, double *water_out, double *oil_out, double *init_oil_vol, double *oil_vol, double *oil_volA ); int main(int argn, char **argc) { Options o = GetOptions(argn, argc); SOTentFunction *tf; Resv_Well_vec wells; double_vec res_sizes, flux[PHASES]; Resv_Func_vec properties = Resv_Func_vec_new(TOTPROP+NUMMAIN); int i, j, plot_skip, outiter = 0, numwells, plotMeshDepth; int plotTotCount[TOTPROP+NUMMAIN], plot_count[TOTPROP+NUMMAIN], num_out_wells; double init_press, init_sat, time, initP[MAXDIM], finP[MAXDIM]; double extra_flux, extra, well_vol, local_sat, reserv_box, init_influx; double dpressMaxVar = 0.0, pressMax, dpressMax = 0.0, satMax, satMin, dsatMax = 0.0, dsatMin, d, norm; SOSimoMap_Data Frelperm[PHASES], Fvisc[PHASES], Fform[PHASES], Fgravden[PHASES], Fabperm, Fpor; FuncMap *FMrelperm[PHASES], *FMvisc[PHASES], *FMform[PHASES], *FMgravden[PHASES], *FMabperm, *FMpor; SOSimoMap_Data *Maps[TOTPROP] = {&(Frelperm[WATER]), &(Fvisc[WATER]), &(Fform[WATER]), &(Fgravden[WATER]), &(Frelperm[OIL]), &(Fvisc[OIL]), &(Fform[OIL]), &(Fgravden[OIL]), &Fabperm, &Fpor}; FuncMap **FMaps[TOTPROP] = {&(FMrelperm[WATER]), &(FMvisc[WATER]), &(FMform[WATER]), &(FMgravden[WATER]), &(FMrelperm[OIL]), &(FMvisc[OIL]), &(FMform[OIL]), &(FMgravden[OIL]), &FMabperm, &FMpor}; /* Gets basis and reservoir domain dimension */ char *basName = txtcat(o.basisName, ".bas"); Basis bas = SOApprox_ReadBasis(basName); unsigned int dim = bas.nel; SOGrid_Dim pDim = bas.el[0]->d->pDim; /* Sets approximate functions and their weight arrays: */ SOFunction *saturation; /* Current Water saturation approximation function. */ SOFunction *dsaturation; /* Current differential Water saturation approximation function. */ SOFunction *pressure; /* Current Oil pressure approximation function. */ SOFunction *Npressure; /* New Oil pressure approximation function. */ SOFunction *dpressure; /* Current Oil pressure difference approximation function. */ SOFunction *Ndpressure; /* New Oil pressure difference approximation function. */ SOFunction *porosity; /* Porosity function. */ double_vec A = double_vec_new(dim); /* Coeffs for the Oil pressure approximation function. */ double_vec Anew = double_vec_new(dim); /* Coeffs for the new Oil pressure approximation function. */ double_vec dA = double_vec_new(dim); /* Coeffs for the Oil pressure difference approximation function. */ double_vec dAnew = double_vec_new(dim); /* Coeffs for the new Oil pressure difference approximation function. */ double_vec dB = double_vec_new(dim); /* Coeffs for the differential Water saturation approximation function. */ double_vec B = double_vec_new(dim); /* Coeffs for the Water saturation approximation function. */ double_vec dAdBnew = double_vec_new(2*dim); /* Coeffs of new solution, A(t) and dB(t). */ double_vec b = double_vec_new(2*dim); /* Right-hand side vector, C_w - M_w and C_o - M_o. */ double_vec PorCoeff = double_vec_new(dim); /* Coeffs of the porosity function. */ /* Gets the {SOGrid} tree */ char treeFile[100]; strcpy(treeFile, o.basisName); i = strlen(treeFile); treeFile[i - 4] = '.'; treeFile[i - 3] = 't'; treeFile[i - 2] = 'r'; treeFile[i - 1] = 'e'; treeFile[i] = 'e'; treeFile[i + 1] = '\0'; FILE *rd = open_read(treeFile); SOGrid_Tree *tree = SOGrid_Tree_read(rd); /* Reads the reservoir data file: */ char *reservoirName = txtcat(o.reservoirName, ".dat"); rd = open_read(reservoirName); assert((rd != NULL), "Needs a reservoir file descriptor."); ReservoirDataRead( rd, pDim, &res_sizes, initP, finP, &init_press, &init_sat, &plot_skip, &properties, &numwells, &wells); reserv_box = 1.0; for(i = 0; i < pDim; i++) reserv_box *= res_sizes.el[i]; /* Builds the initial pressure {A} and saturation {B} functions. */ flux[WATER] = double_vec_new(dim); flux[OIL] = double_vec_new(dim); for(i = 0; i < dim; i++) { A.el[i] = init_press; Anew.el[i] = init_press; dA.el[i] = 0.0; B.el[i] = init_sat; dB.el[i] = 0.0; flux[WATER].el[i] = 0.0; flux[OIL].el[i] = 0.0; PorCoeff.el[i] = 1.0; } pressMax = init_press; satMax = init_sat; pressure = SOApprox_BuildIterFunction(bas, A, basName, pDim, 1); Npressure = SOApprox_BuildIterFunction(bas, Anew, basName, pDim, 1); dpressure = SOApprox_BuildIterFunction(bas, dA, basName, pDim, 1); Ndpressure = SOApprox_BuildIterFunction(bas, dAnew, basName, pDim, 1); saturation = SOApprox_BuildIterFunction(bas, B, basName, pDim, 1); dsaturation = SOApprox_BuildIterFunction(bas, dB, basName, pDim, 1); porosity = SOApprox_BuildIterFunction(bas, PorCoeff, basName, pDim, 1); /* Sets well positions and initial values. */ int wellind[numwells]; double wellpos[numwells * pDim]; for(i = 0; i < numwells * pDim; i++) wellpos[i] = wells.el[(int)i/pDim].pos[i%pDim]; /* Gets the pre-calculated dot product matrices. */ double *TotH = malloc(sizeof(double) * (2*dim) * (2*dim + 1)); SOMatrix Hev = SOMatrix_Null(dim,dim); SOMatrix Hgr = SOMatrix_Null(dim,dim); SOMatrix H; H.ents = MatEntry_vec_new((2 * dim) * (2 * dim)); H.rows = (2 * dim); H.cols = (2 * dim); Hev = SOApprox_ReadMatrix(txtcat(o.matName, "-ev.mat")); Hgr = SOApprox_ReadMatrix(txtcat(o.matName, "-gr.mat")); SOApprox_GetBasisIndeces(bas, pDim, numwells, wellind, wellpos); double min_range[pDim], max_range[pDim], tot_flux = 0; double_vec qflux = double_vec_new(numwells); num_out_wells = 0; for(i = 0; i < numwells; i++) if(wells.el[i].type == 0) { tot_flux += wells.el[i].flux; flux[WATER].el[wellind[i]] = wells.el[i].flux; printf(" DEBUG INJ >>>> flux[WATER].el[%d]: %.16g \n\n\n", wellind[i], flux[WATER].el[wellind[i]]); } else { num_out_wells ++; well_vol = 1.0; tf = (SOTentFunction *)bas.el[wellind[i]]; SOGrid_cell_coords(tf->d->index, tf->d->rank, pDim, min_range, max_range); for(j = 0; j < pDim; j++) well_vol *= 2.0 * (max_range[j] - min_range[j]) * res_sizes.el[j]; qflux.el[i] = (double)wells.el[i].flux * (double)(1 << pDim) / well_vol; printf("\n TOTAL EXT >>>>( wells.el[%d].flux: %.16g tstep: %.16g )-> vol: %.16g well_vol: %.16g \n\n", wellind[i], wells.el[i].flux, o.tStep, (double)(wells.el[i].flux * o.tStep), well_vol); flux[WATER].el[wellind[i]] = -(double)wells.el[i].flux * init_sat; flux[OIL].el[wellind[i]] = -(double)wells.el[i].flux * (1.0 - init_sat); printf(" DEBUG EXT >>>> flux[WATER].el[%d]: %.16g \n\n\n", wellind[i], flux[WATER].el[wellind[i]]); printf(" DEBUG EXT >>>> flux[OIL].el[%d]: %.16g qflux.el[%d]: %.16g \n\n\n", wellind[i], flux[OIL].el[wellind[i]], i, qflux.el[i]); } /* Sets the property functions maps. */ for(i = 0; i < TOTPROP; i++) { *(Maps[i]) = SOSimoMap_FromName(properties.el[i].descr, properties.el[i].degree, properties.el[i].coeff, pDim); *(FMaps[i]) = &((*(Maps[i])).FMap); properties.el[i].FMap = *(FMaps[i]); } SetPropFunctions(&properties, pressure, saturation, dpressure, dsaturation); FuncMap NoMap = IdentityMap(1, pDim); properties.el[TOTPROP+NUMMAIN - 5].F = porosity; /* Porosity. */ properties.el[TOTPROP+NUMMAIN - 4].FMap = &NoMap; properties.el[TOTPROP+NUMMAIN - 3].FMap = &NoMap; properties.el[TOTPROP+NUMMAIN - 2].FMap = &NoMap; properties.el[TOTPROP+NUMMAIN - 1].FMap = &NoMap; /* Plot initialization. */ plotMeshDepth = (int)ceil(2*log(o.plt.figSize/o.plt.meshSize)/log(2)); for(i = 0; i < TOTPROP+NUMMAIN; i++) { plot_count[i] = 0; if(properties.el[i].plot || i >= TOTPROP) OpenPlot(&(properties.el[i]), &o, &(plotTotCount[i])); } /* Estimates the average time for simulation. */ double tot_vol, init_oil_vol, oil_vol, oil_volA, water_vol, water_volA; SOProcFunction *unit = SOProcFunction_FromName("unit", pDim, 1); SOFunction *uf = SOFunction_Cast(unit); /* Gets the total volume. */ SOFunction_RawDot( porosity, &NoMap, uf, &NoMap, NULL, tree, &tot_vol ); tot_vol *= reserv_box; TotalOilVol( saturation, &NoMap, porosity, FMpor, tree, pDim, &tot_vol, &water_vol, &oil_vol, &reserv_box ); init_oil_vol = oil_vol; double tot_est_time = (double) init_oil_vol/tot_flux; double oil_out = 0.0, water_out = 0.0; printf("\n\n ************ ESTIMATED TOTAL TIME FOR SIMULATION: %.16g DAYS ********************\n\n", tot_est_time/DAYSECONDS); time = 0.0; while (time <= o.tFin) { fprintf(stderr, "\n=== time = %f DAYS ===============================\n", time/DAYSECONDS); /* Given the saturation field {B} for this instant, computes the pressure field difference {dA} and water saturation change {dB}. */ int iter = 0; while (TRUE) { H.ents.nel = 0; SOApprox_SimoGetLeftMatrix /* Obtains the system's matrix: */ ( bas, &H, Hev, Hgr, reserv_box, saturation, Npressure, porosity, FMpor, FMrelperm, FMvisc, FMform, FMabperm, FALSE ); SOApprox_SimoGetRightHandSide /* Obtains the system's right-hand side: */ ( bas, tree, reserv_box, Hev, saturation, Npressure, porosity, FMpor, flux, res_sizes.el, FMrelperm, FMvisc, FMform, FMgravden, FMabperm, b ); /* Solve the system obtaining the new guesses {dAdBnew}: */ if (o.gaussElim) { SOApprox_GaussElim(TotH, H, b, dAdBnew, FALSE); } else { assert(FALSE , "no solution method?"); } /* for(i = 0; i < dAdBnew.nel; i++) printf(" DEBUG: dAdBnew[%d]-> %.16g \n", i, dAdBnew.el[i]); */ /* for(i = 0; i < A.nel; i++) printf(" DEBUG PRESSAO: A[%d]-> %.16g \n", i, A.el[i]); */ /* for(i = 0; i < b.nel; i++) printf(" DEBUG B: b[%d]-> %.16g \n", i, b.el[i]); */ /* for(i = 0; i < H.ents.nel; i++) printf(" DEBUG: H.ents.el[%d].va-> %.16g \n", i, H.ents.el[i].va); */ /* Extracts the saturation and pressure changes: */ for(i = 0; i < dim; i++) { dAnew.el[i] = dAdBnew.el[i]; dB.el[i] = dAdBnew.el[i + dim]; } printf("\n OIL PRESSURE(iter:%d): \n", iter); SOApprox_PrintMaxErrorValues(Ndpressure, dpressure, &dpressMax, &dpressMaxVar); SOApprox_GetLimitValues(dsaturation, &dsatMax, &dsatMin); /* Prepare for next iteration: */ printf("Distance change: %f to ", d); d = SOApprox_UpdateSolution(dAnew, dA); printf(" %f \n", d); norm = rn_norm(dim, dA.el); for(i = 0; i < dim; i++) Anew.el[i] = A.el[i] + dA.el[i]; /* Checks pressure convergence. */ printf("RelTol: %f, AbsTol: %f, Norm: %f, dpressMaxVar: %f \n", o.relTol, o.absTol, norm, dpressMaxVar); if ((iter > 0) && ((d <= o.relTol*norm) && (d <= o.absTol))) { break; } if ((o.plotAll) && !(iter % plot_skip)) { for(i = 0; i < TOTPROP; i++) if(properties.el[i].plot) PlotState( &(properties.el[i]), pDim, tree, initP, finP, &o, &plotMeshDepth, "RUNNING", time, iter, outiter, NULL, &(plot_count[i]), &(plotTotCount[i]), &water_out, &oil_out, &init_oil_vol, &oil_vol, &oil_volA); PlotState( &(properties.el[TOTPROP+NUMMAIN-4]), pDim, tree, initP, finP, &o, &plotMeshDepth, "RUNNING", time, iter, outiter, &pressMax, &(plot_count[i]), &(plotTotCount[i]), &water_out, &oil_out, &init_oil_vol, &oil_vol, &oil_volA); PlotState( &(properties.el[TOTPROP+NUMMAIN-3]), pDim, tree, initP, finP, &o, &plotMeshDepth, "RUNNING", time, iter, outiter, &satMax, &(plot_count[i]), &(plotTotCount[i]), &water_out, &oil_out, &init_oil_vol, &oil_vol, &oil_volA); PlotState( &(properties.el[TOTPROP+NUMMAIN-2]), pDim, tree, initP, finP, &o, &plotMeshDepth, "RUNNING", time, iter, outiter, &dpressMax, &(plot_count[i]), &(plotTotCount[i]), &water_out, &oil_out, &init_oil_vol, &oil_vol, &oil_volA); PlotState( &(properties.el[TOTPROP+NUMMAIN-1]), pDim, tree, initP, finP, &o, &plotMeshDepth, "RUNNING", time, iter, outiter, &dsatMax, &(plot_count[i]), &(plotTotCount[i]), &water_out, &oil_out, &init_oil_vol, &oil_vol, &oil_volA); } iter++; } /* Now the new saturation is generated. */ for(i = 0; i < dim; i++){B.el[i] = B.el[i] + dB.el[i] * o.tStep; A.el[i] = Anew.el[i];} SOApprox_GetLimitValues(saturation, &satMax, &satMin); fprintf(stderr, "Maximum SATURATION, at TIME %f: %16.12f \n", time, satMax); /* extra = 0.0; */ /* if(satMax > 0.9) */ /* { extra = satMax - 0.9; */ /* if(extra > 0.2) extra = 0.2; */ /* } */ /* if(satMin < 0.0) */ /* { if(satMin < -extra) extra = fabs(satMin); */ /* if(extra > 0.2) extra = 0.2; */ /* } */ /* extra_flux = (0.2 - extra)/0.2; */ /* if(outiter == 0) {extra_flux = 1.0;} */ /* Gets the new oil total volume.*/ oil_volA = oil_vol; water_volA = water_vol; TotalOilVol( saturation, &NoMap, porosity, FMpor, tree, pDim, &tot_vol, &water_vol, &oil_vol, &reserv_box ); printf("<> oil_volA: %.16g water_volA: %.16g \n", oil_volA, water_volA); printf("<> oil_vol: %.16g water_vol: %.16g \n", oil_vol, water_vol); printf("<> oil_vol - oil_volA: %.16g water_vol - water_volA: %.16g \n\n", oil_vol - oil_volA, water_vol - water_volA); if(outiter == 0) {init_influx = water_vol - water_volA; extra_flux = 1.0;} else extra_flux = fabs((double)(water_vol - water_volA) / init_influx); if(extra_flux < 0.5)extra_flux = 1.0; printf("\n\n>>DBG: EXTRA_FLUX: %.16g, EXTRA: %.16f, outiter: %d, satMax: %.16g satMin: %.16g \n\n", extra_flux, extra, outiter, satMax, satMin); /* Producer well flux is adapted: */ for(i = 0; i < numwells; i++) if(wells.el[i].type == 1) { /* Counts the flux taken out. */ water_out += - flux[WATER].el[wellind[i]] * o.tStep; oil_out += - flux[OIL].el[wellind[i]] * o.tStep; printf("<> flux[OIL] * tstep: %.16g (oil_volA - oil_vol): %.16g \n\n", - flux[OIL].el[wellind[i]] * o.tStep, oil_volA - oil_vol); local_sat = 0.0; for(j = 0; j < Hev.ents.nel; j++) if(Hev.ents.el[j].row == wellind[i]) { local_sat += Hev.ents.el[j].va * B.el[Hev.ents.el[j].col]; /* printf(" DEBUG HEV 1 >> j: %d, row: %d, col: %d val: %.16g B.el[%d]: %.16g local_sat: %.16g \n\n", */ /* j, Hev.ents.el[j].row, Hev.ents.el[j].col, Hev.ents.el[j].va, */ /* Hev.ents.el[j].col, B.el[Hev.ents.el[j].col], local_sat); */ } /* Resets water flux. */ flux[WATER].el[wellind[i]] = - qflux.el[i] * reserv_box * local_sat; printf("<> WATER qflux.el[%d]: %.16g reserv_box: %.16g local_sat: %.16g \n\n", i, qflux.el[i], reserv_box, local_sat); local_sat = 0.0; for(j = 0; j < Hev.ents.nel; j++) if(Hev.ents.el[j].row == wellind[i]) local_sat += Hev.ents.el[j].va * (1.0 - B.el[Hev.ents.el[j].col]); /* Resets oil flux. */ flux[OIL].el[wellind[i]] = - qflux.el[i] * reserv_box * local_sat; printf("<> OIL qflux.el[%d]: %.16g reserv_box: %.16g local_sat: %.16g \n\n", i, qflux.el[i], reserv_box, local_sat); /* Flux correction. */ /* if(flux[OIL].el[wellind[i]] > 0.0) flux[OIL].el[wellind[i]] = 0.0; */ /* if(flux[WATER].el[wellind[i]] > 0.0) flux[WATER].el[wellind[i]] = 0.0; */ flux[OIL].el[wellind[i]] *= extra_flux; flux[WATER].el[wellind[i]] *= extra_flux; /* printf("<> oil_volA: %.16g oil_vol: %.16g water_vol: %.16g \n", oil_volA, oil_vol, water_vol); */ /* printf("<> flux[WATER].el[%d]: %.16g \n", wellind[i], flux[WATER].el[wellind[i]]); */ /* printf("<> flux[OIL].el[%d]: %.16g \n", wellind[i], flux[OIL].el[wellind[i]]); */ /* printf("<> flux[OIL]+flux[WATER]: %.16g \n", flux[OIL].el[wellind[i]] + flux[WATER].el[wellind[i]]); */ /* /\* Relative permeability adjustment. *\/ */ /* if(flux[WATER].el[wellind[i]] / (flux[WATER].el[wellind[i]]+flux[OIL].el[wellind[i]]) < MIN_RELPERM) */ /* { */ /* printf("<>: %.16g \n", */ /* flux[WATER].el[wellind[i]] / (flux[WATER].el[wellind[i]]+flux[OIL].el[wellind[i]])); */ /* flux[OIL].el[wellind[i]] += flux[WATER].el[wellind[i]]; */ /* flux[WATER].el[wellind[i]] = 0.0; */ /* printf("<> flux[WATER].el[%d]: %.16g \n", wellind[i], flux[WATER].el[wellind[i]]); */ /* printf("<> flux[OIL].el[%d]: %.16g \n", wellind[i], flux[OIL].el[wellind[i]]); */ /* } */ } else flux[WATER].el[wellind[i]] = wells.el[i].flux * extra_flux; if (!(outiter % plot_skip)) { for(i = 0; i < TOTPROP; i++) if(properties.el[i].plot) PlotState( &(properties.el[i]), pDim, tree, initP, finP, &o, &plotMeshDepth, "TIME_STEP", time, iter, outiter, NULL, &(plot_count[i]), &(plotTotCount[i]), &water_out, &oil_out, &init_oil_vol, &oil_vol, &oil_volA); PlotState( &(properties.el[TOTPROP+NUMMAIN-4]), pDim, tree, initP, finP, &o, &plotMeshDepth, "TIME_STEP", time, iter, outiter, &pressMax, &(plot_count[i]), &(plotTotCount[i]), &water_out, &oil_out, &init_oil_vol, &oil_vol, &oil_volA); PlotState( &(properties.el[TOTPROP+NUMMAIN-3]), pDim, tree, initP, finP, &o, &plotMeshDepth, "TIME_STEP", time, iter, outiter, &satMax, &(plot_count[i]), &(plotTotCount[i]), &water_out, &oil_out, &init_oil_vol, &oil_vol, &oil_volA); PlotState( &(properties.el[TOTPROP+NUMMAIN-2]), pDim, tree, initP, finP, &o, &plotMeshDepth, "TIME_STEP", time, iter, outiter, &dpressMax, &(plot_count[i]), &(plotTotCount[i]), &water_out, &oil_out, &init_oil_vol, &oil_vol, &oil_volA); PlotState( &(properties.el[TOTPROP+NUMMAIN-1]), pDim, tree, initP, finP, &o, &plotMeshDepth, "TIME_STEP", time, iter, outiter, &dsatMax, &(plot_count[i]), &(plotTotCount[i]), &water_out, &oil_out, &init_oil_vol, &oil_vol, &oil_volA); } outiter++; iter = 0; time += o.tStep; } for(i = 0; i < TOTPROP+NUMMAIN; i++) if(properties.el[i].plot || i >= TOTPROP) SOPlot_CloseStream(properties.el[i].ps); return 0; } #define PPUSAGE SOParams_SetUsage Options GetOptions(int argn, char **argc) { Options o; SOParams_T *pp = SOParams_NewT(stderr, argn, argc); PPUSAGE(pp, "SOSimoleo \\\n"); PPUSAGE(pp, " -basisName NAME [ -matName NAME ] \\\n"); PPUSAGE(pp, " -reservoirName \\\n"); PPUSAGE(pp, " -tFin TOTAL_TIME \\\n"); PPUSAGE(pp, " -tStep TIME_INTERVAL \\\n"); PPUSAGE(pp, " [ -gaussSeidel [-omega NUM] | \\\n"); PPUSAGE(pp, " -conjGrad | \\\n"); PPUSAGE(pp, " -gaussElim ] \\\n"); PPUSAGE(pp, " [ -maxIter NUM ] [ -relTol NUM ] [ -absTol NUM ] \\\n"); PPUSAGE(pp, " [ -gaussOrder NUM ] \\\n"); PPUSAGE(pp, " -outName NAME \\\n"); PPUSAGE(pp, " [ -plotAll ] [ -plotFinal ] \\\n"); PPUSAGE(pp, SOPlotParams_FunctionHelp " \n"); SOParams_GetKeyword(pp, "-basisName"); o.basisName = SOParams_GetNext(pp); if (SOParams_KeywordPresent(pp, "-matName")) { o.matName = SOParams_GetNext(pp); } else { o.matName = o.basisName; } SOParams_GetKeyword(pp, "-reservoirName"); o.reservoirName = SOParams_GetNext(pp); SOParams_GetKeyword(pp, "-outName"); o.outName = SOParams_GetNext(pp); o.gaussSeidel = SOParams_KeywordPresent(pp, "-gaussSeidel"); o.conjGrad = SOParams_KeywordPresent(pp, "-conjGrad"); o.gaussElim = SOParams_KeywordPresent(pp, "-gaussElim"); if (! (o.gaussSeidel || o.conjGrad || o.gaussElim)) { SOParams_Error(pp, "must specify a linear solution method"); } else if ((o.gaussSeidel + o.conjGrad + o.gaussElim) > 1) { SOParams_Error(pp, "please specify only one linear solution method"); } if (o.gaussSeidel) { if (SOParams_KeywordPresent(pp, "-omega")) { o.omega = SOParams_GetNextDouble(pp, 0.0, MAXDOUBLE); } else { o.omega = 1.0; } } if (SOParams_KeywordPresent(pp, "-maxIter")) { o.maxIter = SOParams_GetNextInt(pp, 0, INT_MAX); } else { o.maxIter = 20; } if (SOParams_KeywordPresent(pp, "-relTol")) { o.relTol = SOParams_GetNextDouble(pp, 0.0, MAXDOUBLE); } else { o.relTol = 0.0; } if (SOParams_KeywordPresent(pp, "-absTol")) { o.absTol = SOParams_GetNextDouble(pp, 0.0, MAXDOUBLE); } else { o.absTol = 0.0; } SOParams_GetKeyword(pp, "-gaussOrder"); o.gaussOrder = SOParams_GetNextInt(pp, 1, INT_MAX); SOParams_GetKeyword(pp, "-tFin"); o.tFin = SOParams_GetNextDouble(pp, 0.0, 1.0e12); SOParams_GetKeyword(pp, "-tStep"); o.tStep = SOParams_GetNextDouble(pp, 0.0, 1.0e8); /* Plotting options: */ o.plotAll = SOParams_KeywordPresent(pp, "-plotAll"); o.plotFinal = SOParams_KeywordPresent(pp, "-plotFinal"); o.plt = SOPlotParams_FunctionParse(pp); SOParams_Finish(pp); return o; } Resv_Func_vec Resv_Func_vec_new(nat nel) { /* This is not a macro only because gcc does not allow cast of struct: */ vec_t v = vec_new(nel, sizeof(Resv_Func)); Resv_Func_vec r; r.nel = v.nel; r.el = (Resv_Func*)v.el; return r; } Resv_Well_vec Resv_Well_vec_new(nat nel) { /* This is not a macro only because gcc does not allow cast of struct: */ vec_t v = vec_new(nel, sizeof(Resv_Well)); Resv_Well_vec r; r.nel = v.nel; r.el = (Resv_Well*)v.el; return r; } void ReservoirDataRead ( FILE *rd, SOGrid_Dim pDim, double_vec *res_sizes, double *initP, double *finP, double *init_press, double *init_sat, int *plot_skip, Resv_Func_vec *properties, int *numwells, Resv_Well_vec *wells ) { int i, resdim; char *tmpstr; filefmt_read_header(rd, "SOReservoir", SOReservoir_FileFormat); resdim = nget_int(rd, "domain_dim"); fget_eol(rd); assert(resdim == pDim, "bad dimension"); *res_sizes = double_vec_new(resdim); nget_name_eq(rd, "sizes_width"); for (i = 0; i < resdim; i++) res_sizes->el[i] = fget_double(rd); fget_eol(rd); nget_name_eq(rd, "plot_points"); for (i = 0; i < resdim; i++) initP[i] = fget_double(rd); for (i = 0; i < resdim; i++) finP[i] = fget_double(rd); fget_eol(rd); *init_press = nget_double(rd, "init_press"); fget_eol(rd); *init_sat = nget_double(rd, "init_sat"); fget_eol(rd); *plot_skip = nget_int(rd, "plot_skip"); fget_eol(rd); fget_skip_formatting_chars(rd); tmpstr = nget_string(rd, "section"); assert(!strcmp(tmpstr, "++++++++++++++++++++++++++++WATER_PROPERTIES++++++++++++++++"), "WATER specifications not found."); free(tmpstr); fget_skip_formatting_chars(rd); ResvFuncRead( rd, &(properties->el[0]), "wREL_PERM" ); ResvFuncRead( rd, &(properties->el[1]), "wVISC" ); ResvFuncRead( rd, &(properties->el[2]), "wFORM_VOL" ); ResvFuncRead( rd, &(properties->el[3]), "wGRAV_DEN" ); tmpstr = nget_string(rd, "section"); assert(!strcmp(tmpstr, "++++++++++++++++++++++++++++OIL_PROPERTIES++++++++++++++++"), "OIL specifications not found."); free(tmpstr); fget_skip_formatting_chars(rd); ResvFuncRead( rd, &(properties->el[4]), "oREL_PERM" ); ResvFuncRead( rd, &(properties->el[5]), "oVISC" ); ResvFuncRead( rd, &(properties->el[6]), "oFORM_VOL" ); ResvFuncRead( rd, &(properties->el[7]), "oGRAV_DEN" ); tmpstr = nget_string(rd, "section"); assert(!strcmp(tmpstr, "++++++++++++++++++++++++++++GENERAL_PROPERTIES++++++++++++++++"), "GENERAL specifications not found."); free(tmpstr); fget_skip_formatting_chars(rd); ResvFuncRead( rd, &(properties->el[8]), "AB_PERM" ); ResvFuncRead( rd, &(properties->el[9]), "POR" ); tmpstr = nget_string(rd, "section"); assert(!strcmp(tmpstr, "++++++++++++++++++++++++++++MAIN_PROPERTIES++++++++++++++++"), "MAIN specifications not found."); free(tmpstr); fget_skip_formatting_chars(rd); ResvFuncRead( rd, &(properties->el[10]), "PRESS" ); ResvFuncRead( rd, &(properties->el[11]), "SAT" ); ResvFuncRead( rd, &(properties->el[12]), "DPRESS" ); ResvFuncRead( rd, &(properties->el[13]), "DSAT" ); tmpstr = nget_string(rd, "section"); assert(!strcmp(tmpstr, "++++++++++++++++++++++++++++WELLS++++++++++++++++"), "WELLS specifications not found."); free(tmpstr); fget_eol(rd); *numwells = nget_int(rd, "numwells"); fget_skip_formatting_chars(rd); *wells = Resv_Well_vec_new(*numwells); for (i = 0; i < *numwells; i++) { ResvWellRead( rd, &(wells->el[i]), resdim ); fget_skip_formatting_chars(rd); } filefmt_read_footer(rd, "SOReservoir"); } void ResvFuncRead( FILE *rd, Resv_Func *property, char *propname ) { int i; property->name = nget_string(rd, "name"); fget_eol(rd); assert(!strcmp(propname, property->name),"Wrong property."); property->descr = nget_string(rd, "descr"); fget_eol(rd); property->degree = nget_int(rd, "degree"); property->coeff = malloc((property->degree + 1) * sizeof(double)); assert(property->coeff != NULL, "No memory for property->coeff allocation."); for (i = 0; i <= property->degree; i++) property->coeff[i] = fget_double(rd); fget_eol(rd); property->plot = nget_int(rd, "plot"); fget_eol(rd); property->fRange = nget_double(rd, "fRange"); fget_eol(rd); property->fStep = nget_double(rd, "fStep"); fget_eol(rd); property->changrng = nget_int(rd, "changrng"); fget_eol(rd); property->ps_eps = nget_int(rd, "ps_eps"); fget_eol(rd); property->plot_times = nget_int(rd, "plot_times"); fget_eol(rd); fget_skip_formatting_chars(rd); } void ResvWellRead( FILE *rd, Resv_Well *well, int resdim ) { int i; well->type = nget_int(rd, "type"); fget_eol(rd); nget_name_eq(rd, "position"); well->pos = malloc(resdim * sizeof(double)); assert(well->pos != NULL, "No memory for well->pos allocation."); for (i = 0; i < resdim; i++) { well->pos[i] = fget_double(rd); } fget_eol(rd); well->flux = nget_double(rd, "flux"); } void SetPropFunctions ( Resv_Func_vec *properties, SOFunction *pressure, SOFunction *saturation, SOFunction *dpressure, SOFunction *dsaturation ) { /* ++++++++++++++++++++++++++++WATER_PROPERTIES++++++++++++++++ */ properties->el[0].F = saturation; /* RELATIVE_PERMEABILITY */ properties->el[1].F = pressure; /* VISCOSITY */ properties->el[2].F = pressure; /* FORMATION_VOLUME */ properties->el[3].F = pressure; /* GRAVITATIONAL_DENSITY */ /* ++++++++++++++++++++++++++++OIL_PROPERTIES++++++++++++++++ */ properties->el[4].F = saturation; /* RELATIVE_PERMEABILITY */ properties->el[5].F = pressure; /* VISCOSITY */ properties->el[6].F = pressure; /* FORMATION_VOLUME */ properties->el[7].F = pressure; /* GRAVITATIONAL_DENSITY */ /*++++++++++++++++++++++++++++GENERAL_PROPERTIES++++++++++++++++*/ properties->el[8].F = pressure; /* ABSOLUTE_PERMEABILITY */ properties->el[9].F = pressure; /* POROSITY */ /*++++++++++++++++++++++++++++MAIN_PROPERTIES++++++++++++++++*/ properties->el[10].F = pressure; /* PRESSURE */ properties->el[11].F = saturation; /* SATURATION */ properties->el[12].F = dpressure; /* DIFFERENTIAL_PRESSURE */ properties->el[13].F = dsaturation; /* DIFFERENTIAL_SATURATION */ } void OpenPlot ( Resv_Func *RF, Options *o, int *plotTotCount ) { char paperSize[50]; int pos = 0, hCount = 1, vCount = 1; if(!RF->ps_eps) { strcpy(paperSize,"a3"); pos = 2; hCount = 2; vCount = 3; } *plotTotCount = hCount * vCount; paperSize[pos] = '\0'; RF->ps = SOPlot_OpenStream(RF->ps_eps, txtcat(RF->name, o->outName), paperSize, 350.0, 350.0 * SQRTHALF); SOPlot_ChangeLayout(RF->ps, paperSize, 300.0, 300.0 * SQRTHALF, 1, hCount, vCount); } void TotalOilVol ( SOFunction *saturation, FuncMap *NoMap, SOFunction *porosity, FuncMap *FMpor, SOGrid_Tree *tree, SOGrid_Dim pDim, double *tot_vol, double *water_vol, double *oil_vol, double *reserv_box ) { *water_vol = 0, *oil_vol = 0; /* Gets the water volume. */ SOFunction_RawDot( saturation, NoMap, porosity, FMpor, NULL, tree, water_vol); *water_vol *= *reserv_box; *oil_vol = *tot_vol - *water_vol; } void PlotState ( Resv_Func *RF, SOGrid_Dim pDim, SOGrid_Tree *tree, double *initP, double *finP, Options *o, int *plotMeshDepth, char *tag, double time, int iter, int outiter, double *fMax, int *plot_count, int *plotTotCount, double *water_out, double *oil_out, double *init_oil_vol, double *oil_vol, double *oil_volA ) { bool split; char tmpstr[350], tmpstr2[150]; double PlotMin = 10e+20, PlotMax = 0, minR[2] = {0.0, 0.0}, maxR[2] = {1.0, 1.0}; // fprintf(stderr, "\n+++ Plot %s +++++++++++++++\n", RF->name); printf("\n+++ Plot %s +++++++++++++++\n", RF->name); printf("iter: %d, outiter: %d, fMax: %f, fRange: %.16g, fStep: %.16g \n", iter, outiter, fMax == NULL ? -1 : *fMax, RF->fRange, RF->fStep); sprintf(tmpstr,"%s - %s \n iter: %d - outiter: %d \n time:%f days", RF->name, tag, iter, outiter, time/DAYSECONDS); if(RF->plot & 2) { SO2DPlot_AddFunctionPicture( RF->ps, RF->F, RF->FMap, tree, minR, maxR, -(RF->fRange), +(RF->fRange), RF->fStep, o->plt.isolineWidth, *plotMeshDepth - 6, 2, *plotMeshDepth, o->plt.gridWidth, (! o->plt.noGrid), MAXDEPTHTREE, (! o->plt.noBands), (! o->plt.noIsolines), &PlotMin, &PlotMax, tmpstr ); (*plot_count)++; } if(fMax == NULL) { sprintf(tmpstr,"%s - %s \n iter: %d - outiter: %d \n time:%f days \n 2D_PlotMax: %f 2D_PlotMin: %f \n", RF->name, tag, iter, outiter, time/DAYSECONDS, PlotMax, PlotMin); sprintf(tmpstr2," init oil vol: %f oil vol: %f \n water out: %f oil out: %f \n vol_diff: %f oil_diff: %f", *init_oil_vol, *oil_vol, *water_out, *oil_out, *init_oil_vol - (*oil_vol + *oil_out), *oil_vol - *oil_volA); strcpy(&(tmpstr[strlen(tmpstr)]), tmpstr2); } else { sprintf(tmpstr,"%s - %s \n iter: %d - outiter: %d \n time:%f days fMax:%f \n 2D_PlotMax: %f 2D_PlotMin: %f \n", RF->name, tag, iter, outiter, time/DAYSECONDS, *fMax, PlotMax, PlotMin); sprintf(tmpstr2," init oil vol: %f oil vol: %f \n water out: %f oil out: %f \n vol_diff: %f oil_diff: %f", *init_oil_vol, *oil_vol, *water_out, *oil_out, *init_oil_vol - (*oil_vol + *oil_out), *oil_vol - *oil_volA); strcpy(&(tmpstr[strlen(tmpstr)]), tmpstr2); } if(RF->plot & 1) { double limit_val = 1.0; if(strcmp(RF->name, "SAT")||strcmp(RF->name, "oREL_PERM") ||strcmp(RF->name, "wREL_PERM")) limit_val = RF->fRange; SO1DPlot_AddFunctionPicture( RF->ps, RF->F, RF->FMap, pDim, initP, finP, -(limit_val * (double)10/8), +(limit_val * (double)10/8), limit_val, PLOTDIVS, o->plt.isolineWidth, &PlotMin, &PlotMax, tmpstr ); (*plot_count)++; } split = fabs((RF->fStep * STEPDIV) - PlotMax) > (RF->fStep/2); if(RF->changrng && split) { if((iter != 0 || outiter != 0) && (RF->fRange <= PlotMax * 100000)) { if((double)PlotMax / STEPDIV < RF->fRange) RF->fStep = (double)PlotMax / STEPDIV; printf(" NEW> fRange: %.16g, PlotMax: %.16g, fStep: %.16g \n\n", RF->fRange, PlotMax, RF->fStep); } if((iter + outiter == 1) && (PlotMax + PlotMin)) { printf("Range change: PlotMax: %f, PlotMin: %f, fRange: %f \n", PlotMax, PlotMin, RF->fRange); RF->fRange = fabs(PlotMax) > fabs(PlotMin) ? fabs(PlotMax) : fabs(PlotMin); RF->fStep = (double)RF->fRange / STEPDIV; } } if(!((*plot_count) % *plotTotCount)) { sprintf(tmpstr,"%08d", iter); SOPlot_AddPage(RF->ps, tmpstr); } fprintf(stderr, "\n++++++++++++++++++++++++\n"); }