/* SODynSimoleo -- Oil reservoir simulator. */ /* Receives as input a reservoir data file consisting of: -> Size of a *square* reservoir; -> Water property functions descriptors; -> Oil property functions descriptors; -> General property functions descriptors; -> 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 #include #include #define PHASES 2 /* The number of fluids of the system: oil and water. */ #define DAYSECONDS 86400.0 /* Number of seconds of a day. */ #define YEARDAYS 365 /* Number of days of an year. */ #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 MAXPDIM 3 #define PLOTDIVS 64 #define STEPDIV 10 #define MAXDEPTHTREE 20 #define MIN_RELPERM 0.1 #define MAX_RELPERM 0.9 #define MAX_TREE_RANK 10 #define DIFF_LIMIT 0.0001 #define SOReservoir_FileFormat "2004-02-10" typedef struct Options { char *matName; /* Prefix for dot products *list*matrix filenames. */ 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. */ bool maximal; /* TRUE for maximal-support tents, FALSE for minimal ones. */ bool changeBas; /* TRUE changes the approximation basis dynamically. */ int refinedskip; /* Number of iterations with refined grids before new refinement. */ int fineskip; /* Number of iterations with fine grids. */ /* 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; 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 *res_side, int *minRank, int *maxRank, double *refinerate, double *initP, double *finP, double *init_press, double *init_sat, int *plot_skip, Resv_Func_vec *Mproperties, 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 *Mproperties, Resv_Func_vec *properties, SOFunction *porosity, SOFunction *pressure, SOFunction *saturation, SOFunction *dsaturation, SOFunction *Qflux ); 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 ); void ShatterNode ( SOGrid_Dim pDim, SOGrid_Node *p, SOGrid_Rank r, double *lo, double *hi, int minRank, int maxRank, double refinerate, int numwells, Resv_Well_vec wells ); void DynShatterNode ( SOGrid_Dim pDim, SOGrid_Node *p, SOGrid_Rank r, double *lo, double *hi, int minRank, int maxRank, double diff_limit, SOFunction *f ); void DynShatterNode2 ( SOGrid_Dim pDim, SOGrid_Node *p, SOGrid_Rank r, double *lo, double *hi, int minRank, int maxRank ); void SetGenDotMatrices ( ListMatEntry_vec *evaldot, ListMatEntry_vec *graddot, SOListMatrix Meval, SOListMatrix Mgrad, unsigned int offset, unsigned int maxindices ); void SetDotMatrices ( unsigned int dim, unsigned int offset, Basis bas, ListMatEntry_vec *evaldot, ListMatEntry_vec *graddot, ListMatEntry_vec *Hev_ents, ListMatEntry_vec *Hgr_ents ); void DynChangeFuncBasis ( SOGrid_Dim pDim, SOGrid_Tree *tree, double minRank, double maxRank, double satMin, double satMax, bool maximal, SOTent_vec *tents, Basis *bas, unsigned int *dim, unsigned int offset, ListMatEntry_vec *evaldot, ListMatEntry_vec *graddot, SOListMatrix *Hev, SOListMatrix *Hgr, SOFunction *pressure, SOFunction *Npressure, SOFunction *dsaturation, SOFunction *saturation, SOFunction *Qflux, SOFunction *porosity, Resv_Func_vec properties, double_vec B, double_vec A, double_vec Aini, double_vec P, double *TotH, double_vec b, int_vec *intersec ); void DynChangeFlux ( SOGrid_Dim pDim, Basis bas, double res_side, SOTent_vec tents, SOTent_vec new_tents, Resv_Well_vec wells, int numwells, unsigned int new_dim ); void ChangeBasis ( unsigned int offset, unsigned int dim, unsigned int new_dim, Basis bas, Basis new_bas, SOListMatrix Hev, ListMatEntry_vec *evaldot, double *TotH, double_vec b, double_vec coeffs ); void SaveResultStep ( SOGrid_Dim pDim, SOGrid_Tree *tree, SOFunction *pressure, double reserv_box, SOFunction *uf, FuncMap *NoMap, double time, int_vec *schedule, double_vec *avg_press, double_vec *oil_flux, double_vec *water_flux, double_vec *oil_vol_out, double_vec *water_vol_out, double_vec *appr_oil_flux, double_vec *appr_water_flux, double_vec *appr_oil_vol_out, double_vec *appr_water_vol_out, double water_flux_out, double oil_flux_out, double init_oil_vol, double water_volA, double water_vol, double oil_volA, double oil_vol, double tot_vol, double water_out, double water_in, double water_vol_in, double oil_out, double tStep, int plot_skip ); void SaveResults ( char *outName, char *init_time, char *end_time, int_vec schedule, double_vec avg_press, double_vec oil_flux, double_vec water_flux, double_vec oil_vol_out, double_vec water_vol_out, double_vec appr_oil_flux, double_vec appr_water_flux, double_vec appr_oil_vol_out, double_vec appr_water_vol_out ); void WriteArray( FILE *wr, char *descr, int_vec schedule, double_vec array ); void TestMatrixHSymmetry(SOMatrix H, unsigned int dim); void TestVecbSymmetry(double_vec b, unsigned int dim); int main(int argn, char **argc) { Options o = GetOptions(argn, argc); SOTentFunction *tf; Resv_Well_vec wells; Resv_Func_vec Mproperties = Resv_Func_vec_new(NUMMAIN); Resv_Func_vec properties = Resv_Func_vec_new(TOTPROP); int minRank, maxRank, i, j, l, plot_skip, outiter = 0, numwells, plotMeshDepth, refinedcount, finecount; int plotTotCount[TOTPROP+NUMMAIN], plot_count[TOTPROP+NUMMAIN]; double init_press, init_sat, time, res_side, initP[MAXPDIM], finP[MAXPDIM], maxP[MAXPDIM]; double reserv_box, refinerate, tent_min[MAXPDIM], tent_max[MAXPDIM]; double pressMaxVar, pressMax, pressMin, satMax, satMin, dsatMax, dsatMin, qfluxMax, 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}; /* Reads the reservoir data file: */ SOGrid_Dim pDim; char *reservoirName = txtcat(o.reservoirName, ".dat"); FILE *rd = open_read(reservoirName); assert((rd != NULL), "Needs a reservoir file descriptor."); ReservoirDataRead( rd, &pDim, &res_side, &minRank, &maxRank, &refinerate, initP, finP, &init_press, &init_sat, &plot_skip, &Mproperties, &properties, &numwells, &wells); /* Builds the initial {SOGrid} tree. */ SOGrid_Tree *tree = SOGrid_Tree_new(pDim); SOGrid_Tree *maxtree = SOGrid_Tree_new(pDim); /* Used for reservoir volumes calculation. */ double lo[pDim], hi[pDim]; for(i = 0; i < pDim; i++){ lo[i] = 0.0; hi[i] = (double)1.0/pow(2.0, (double)i/pDim); } /* Finds the lower rank grid to build the initial functions. */ ShatterNode (pDim, tree->root, 0, lo, hi, minRank, minRank, refinerate, numwells, wells); /* Finds the highest rank grid to use with integrals. */ ShatterNode (pDim, maxtree->root, 0, lo, hi, MAX_TREE_RANK, MAX_TREE_RANK, refinerate, numwells, wells); /* Gets basis and reservoir domain dimension */ SOTent_vec tents, new_tents; Basis bas; unsigned int dim, maxdim = (1 << maxRank), max_cells = (1 << pDim); if (o.maximal) tents = SOTent_maximal_basis(tree); else tents = SOTent_minimal_basis(tree); bas = SOFunctionRef_vec_new(tents.nel); dim = bas.nel; char *basName = txtcat(o.outName, ".bas"); for (i = 0; i < tents.nel; i++) { SOTent t = tents.el[i]; bas.el[i] = (SOFunction *)SOTentFunction_Make(pDim, t.cx, t.r); } /* 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 *porosity; /* Porosity function. */ SOFunction *Qflux; /* Flux by volume function. */ double_vec Aini = double_vec_new(maxdim); /* Coeffs for the initial oil pressure function. */ double_vec A = double_vec_new(maxdim); /* Coeffs for the oil pressure approximation function. */ double_vec Anew = double_vec_new(maxdim); /* Coeffs for the new oil pressure approximation function. */ double_vec dB = double_vec_new(maxdim); /* Coeffs for the differential water saturation approximation function. */ double_vec B = double_vec_new(maxdim); /* Coeffs for the water saturation approximation function. */ double_vec dAdBnew = double_vec_new(2*maxdim); /* Coeffs of new solution, A(t) and dB(t). */ double_vec b = double_vec_new(2*maxdim); /* Right-hand side vector, C_w - M_w and C_o - M_o. */ double_vec P = double_vec_new(maxdim); /* Coeffs of the porosity function. */ double_vec Q = double_vec_new(maxdim); /* Coeffs of the flux by volume function. */ /* Result values. */ int tintervals = floor(o.tFin / o.tStep) + 100; char *init_time, *end_time; int_vec schedule = int_vec_new(tintervals); schedule.nel = 0; double_vec avg_press = double_vec_new(tintervals); avg_press.nel = 0; double_vec oil_flux = double_vec_new(tintervals); oil_flux.nel = 0; double_vec water_flux = double_vec_new(tintervals); water_flux.nel = 0; double_vec oil_vol_out = double_vec_new(tintervals); oil_vol_out.nel = 0; double_vec water_vol_out = double_vec_new(tintervals); water_vol_out.nel = 0; double_vec appr_oil_flux = double_vec_new(tintervals); appr_oil_flux.nel = 0; double_vec appr_water_flux = double_vec_new(tintervals); appr_water_flux.nel = 0; double_vec appr_oil_vol_out = double_vec_new(tintervals); appr_oil_vol_out.nel = 0; double_vec appr_water_vol_out = double_vec_new(tintervals); appr_water_vol_out.nel = 0; reserv_box = 1.0; for(i = 0; i < pDim; i++) reserv_box *= res_side; /* Builds the initial pressure {A} and saturation {B} functions. */ for(i = 0; i < maxdim; i++) { Aini.el[i] = 0.0; A.el[i] = 0.0; Anew.el[i] = 0.0; B.el[i] = 0.0; dB.el[i] = 0.0; P.el[i] = 0.0; Q.el[i] = 0.0; } if (o.maximal) { SOFunction *init_sat_max = SOApprox_ReadFunction("init_sat.fun"); SOLinCombFunction *init_sat_max_lc = SOLinCombFunction_Cast(init_sat_max); SOFunction *init_press_max = SOApprox_ReadFunction("init_press.fun"); SOLinCombFunction *init_press_max_lc = SOLinCombFunction_Cast(init_press_max); for(i = 0; i < dim; i++) { Aini.el[i] = init_press_max_lc->d->coef.el[i]; P.el[i] = 1.0; B.el[i] = init_sat_max_lc->d->coef.el[i]; } } else for(i = 0; i < dim; i++) { Aini.el[i] = init_press; P.el[i] = 1.0; B.el[i] = init_sat; } pressMax = init_press; satMax = init_sat; pressure = SOApprox_BuildIterFunction(bas, A, basName, pDim, 1); Npressure = SOApprox_BuildIterFunction(bas, Anew, basName, pDim, 1); saturation = SOApprox_BuildIterFunction(bas, B, basName, pDim, 1); dsaturation = SOApprox_BuildIterFunction(bas, dB, basName, pDim, 1); porosity = SOApprox_BuildIterFunction(bas, P, basName, pDim, 1); Qflux = SOApprox_BuildIterFunction(bas, Q, basName, pDim, 1); /* pressure = SOApprox_BuildFunction(bas, Aini, "init_press", basName, dsaturation); */ /* saturation = SOApprox_BuildFunction(bas, B, "init_sat", basName, dsaturation); */ /* 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(&Mproperties, &properties, porosity, pressure, saturation, dsaturation, Qflux); FuncMap NoMap = IdentityMap(1, pDim); for(i = 0; i < NUMMAIN; i++) Mproperties.el[i].FMap = &NoMap; /* Gets the pre-calculated dot product matrices. */ double *TotH = malloc(sizeof(double) * (2*maxdim) * (2*maxdim + 1)); SOMatrix H; SOListMatrix Hev, Hgr, Meval, Mgrad; H.ents = MatEntry_vec_new((2 * maxdim) * (2 * maxdim)); H.rows = (2 * maxdim); H.cols = (2 * maxdim); /* and get all dot product combinations between minRank to maxRank basis. */ Meval = SOApprox_ReadListMatrix(txtcat(o.matName, "-ev.lmat")); Mgrad = SOApprox_ReadListMatrix(txtcat(o.matName, "-gr.lmat")); unsigned int offset = (1 << minRank), maxindices; /* and may have only up to the total maxRank basis (grid all refined). */ Hev.ents = ListMatEntry_vec_new(maxdim * maxdim); Hev.rows = maxdim; Hev.cols = maxdim; Hev.length = max_cells; Hgr.ents = ListMatEntry_vec_new(maxdim * maxdim); Hgr.rows = maxdim; Hgr.cols = maxdim; Hgr.length = max_cells; /* and keep and in a direct access organization. */ maxindices = 0; for(i = minRank; i <= maxRank; i++) maxindices += (1 << i); if (o.maximal) { offset = (1 << pDim); maxindices = 0; for(i = pDim; i <= maxRank; i++) maxindices += (1 << i); } ListMatEntry_vec evaldot[maxindices], graddot[maxindices]; SetGenDotMatrices(evaldot, graddot, Meval, Mgrad, offset, maxindices); /* Gets the dot product values for the main matrices and . */ SetDotMatrices( dim, offset, bas, evaldot, graddot, &(Hev.ents), &(Hgr.ents) ); /* Builds the intersection array for speeding up linear combinations. */ int_vec intersec[maxdim]; for(i = 0; i < maxdim; i++) intersec[i] = int_vec_new(maxdim/4); /* Rearranges the main grid and its derived functions. */ if( tree->root->c[LO] != NULL && tree->root->c[HI] != NULL ) { SOGrid_subtree_free(tree->root->c[LO]); tree->root->c[LO] = NULL; SOGrid_subtree_free(tree->root->c[HI]); tree->root->c[HI] = NULL; } ShatterNode (pDim, tree->root, 0, lo, hi, o.changeBas?minRank:maxRank, maxRank, refinerate, numwells, wells); DynChangeFuncBasis( pDim, tree, o.changeBas?minRank:maxRank, maxRank, 0.0, init_sat, o.maximal, &new_tents, &bas, &dim, offset, evaldot, graddot, &Hev, &Hgr, pressure, Npressure, dsaturation, saturation, Qflux, porosity, properties, B, A, Aini, P, TotH, b, intersec ); free((SOTent *)tents.el); tents.el = new_tents.el; tents.nel = new_tents.nel; /* 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]; SOApprox_GetBasisIndeces(bas, pDim, numwells, wellind, wellpos); for(i = 0; i < numwells; i++) wells.el[i].ind = wellind[i]; double tot_flux = 0; for(i = 0; i < numwells; i++) if(wells.el[i].type == 0) { tot_flux += wells.el[i].tflux; tf = (SOTentFunction *)bas.el[wells.el[i].ind]; SOGrid_cell_coords(tf->d->index, tf->d->rank, pDim, tent_min, tent_max); for(j = 0; j < pDim; j++) /* Canonical geometry. */ { tent_min[j] = tent_min[j] / pow(2.0, (double)j/pDim); tent_max[j] = tent_max[j] / pow(2.0, (double)j/pDim); } wells.el[i].vol = 1.0; for(j = 0; j < pDim; j++) wells.el[i].vol *= 2.0 * (tent_max[j] - tent_min[j]) * res_side; wells.el[i].vol = (double)wells.el[i].vol / (1 << pDim); wells.el[i].qflux = (double)wells.el[i].tflux / wells.el[i].vol; Q.el[wells.el[i].ind] = wells.el[i].qflux; wells.el[i].wflux = wells.el[i].tflux; fprintf(stderr, " Well %d, INJ> flux(WTR): %.16g, tent: %d, vol: %.16g, qflux: %.16g\n\n", i, wells.el[i].wflux, wells.el[i].ind, wells.el[i].vol, wells.el[i].qflux); } else { tf = (SOTentFunction *)bas.el[wells.el[i].ind]; SOGrid_cell_coords(tf->d->index, tf->d->rank, pDim, tent_min, tent_max); for(j = 0; j < pDim; j++) /* Canonical geometry. */ { tent_min[j] = tent_min[j] / pow(2.0, (double)j/pDim); tent_max[j] = tent_max[j] / pow(2.0, (double)j/pDim); } wells.el[i].vol = 1.0; for(j = 0; j < pDim; j++) wells.el[i].vol *= 2.0 * (tent_max[j] - tent_min[j]) * res_side; wells.el[i].vol = (double)wells.el[i].vol / (1 << pDim); wells.el[i].qflux = (double)wells.el[i].tflux / wells.el[i].vol; Q.el[wells.el[i].ind] = -wells.el[i].qflux; wells.el[i].wflux = 0.0; //-(double)wells.el[i].tflux * init_sat; wells.el[i].oflux = -(double)wells.el[i].tflux; //* (1.0 - init_sat); fprintf(stderr, " Well %d, EXT> flux(OIL): %.16g, flux(WTR): %.16g, tent: %d, vol: %.16g, qflux: %.16g \n\n", i, wells.el[i].oflux, wells.el[i].wflux, wells.el[i].ind, wells.el[i].vol, wells.el[i].qflux); fprintf(stderr, " Well CAPACITY> (wells.el[%d].flux: %.16g, tstep: %.16g)-> vol: %.16g well_vol: %.16g \n\n", wells.el[i].ind, wells.el[i].tflux, o.tStep, wells.el[i].tflux * o.tStep, wells.el[i].vol); } qfluxMax = fabs(wells.el[0].qflux); /* Plot initialization. */ plotMeshDepth = (int)ceil(2*log(o.plt.figSize/o.plt.meshSize)/log(2)); for(i = 0; i < TOTPROP; i++) { plot_count[i] = 0; if(properties.el[i].plot) OpenPlot(&(properties.el[i]), &o, &(plotTotCount[i])); } for(i = 0; i < NUMMAIN; i++) { plot_count[TOTPROP+i] = 0; OpenPlot(&(Mproperties.el[i]), &o, &(plotTotCount[TOTPROP+i])); } /* Estimates the average time for simulation. */ double tot_vol, init_oil_vol, oil_vol, oil_volA, water_vol, water_volA; double oil_flux_out, water_flux_out, water_vol_in = 0.0; /* Gets the total volume. */ SOProcFunction *unit = SOProcFunction_FromName("unit", pDim, 1); SOFunction *uf = SOFunction_Cast(unit); SOFunction_RawDot( porosity, FMpor, uf, &NoMap, NULL, maxtree, &tot_vol ); tot_vol *= reserv_box; TotalOilVol( saturation, &NoMap, porosity, FMpor, maxtree, pDim, &tot_vol, &water_vol, &oil_vol, &reserv_box ); init_oil_vol = oil_vol; fprintf(stderr, " RESERVOIR Total Vol> %.16g \n\n", tot_vol); fprintf(stderr, " INITIAL OIL Vol> %.16g \n\n", oil_vol); fprintf(stderr, " INITIAL WATER Vol> %.16g \n\n", water_vol); fprintf(stderr, " TOTAL (IN/OUT) FLUX> %.16g m^3/s or %.16g m^3/day \n\n", tot_flux, tot_flux * DAYSECONDS); double tot_est_time = (double) init_oil_vol/tot_flux; double oil_out = 0.0, water_out = 0.0, water_in = 0.0; /* Save initial result values. */ SaveResultStep( pDim, maxtree, pressure, reserv_box, uf, &NoMap, time, &schedule, &avg_press, &oil_flux, &water_flux, &oil_vol_out, &water_vol_out, &appr_oil_flux, &appr_water_flux, &appr_oil_vol_out, &appr_water_vol_out, 0.0, wells.el[0].tflux, init_oil_vol, water_vol, water_vol, oil_vol, oil_vol, tot_vol, 0.0, 0.0, 0.0, 0.0, o.tStep, plot_skip ); fprintf(stderr, "\n\n ********* ESTIMATED TOTAL TIME FOR SIMULATION: %d YEARS and %d DAYS ************\n\n", (int)(tot_est_time/DAYSECONDS)/YEARDAYS, (int)(tot_est_time/DAYSECONDS)%YEARDAYS); init_time = today(); time = 0.0; refinedcount = 0; finecount = 0; while (time <= o.tFin && outiter < 222) { printf("\n=== time = %f DAYS, outiter = %d ===============================\n", time/DAYSECONDS, outiter); /* Given the saturation field {B} for this instant, computes the pressure {A} and water saturation changes {dB}. */ int iter = 0; while (TRUE) { H.ents.nel = 0; H.rows = 2 * dim; H.cols = 2 * dim; SOApprox_SimoGetLeftMatrix /* Obtains the system's matrix: */ ( bas, &H, Hev, Hgr, res_side, saturation, pressure, porosity, FMpor, FMrelperm, FMvisc, FMform, FMabperm, intersec ); SOApprox_SimoGetRightHandSide /* Obtains the system's right-hand side: */ ( bas, Hev, reserv_box, wells, saturation, pressure, 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?"); } /* ########## DEBUG #################### */ /* TestMatrixHSymmetry(H, dim); */ /* TestVecbSymmetry(b, dim); */ /* TestVecbSymmetry(dAdBnew, dim); */ /* for(i = 0; i < 2 * dim; i++) printf(" DEBUG: dAdBnew[%d]-> %f \n", i, dAdBnew.el[i]); */ /* for(i = 0; i < dim; i++) printf(" DEBUG PRESSAO: A[%d]-> %.16g \n", i, A.el[i]); */ /* for(i = 0; i < 2 * dim; i++) printf(" DEBUG b: b[%d]-> %f \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++) { Anew.el[i] = Aini.el[i] + dAdBnew.el[i]; dB.el[i] = dAdBnew.el[i + dim]; } SOApprox_GetLimitValues(pressure, maxP, &pressMax, &pressMin); fprintf(stderr, "PRESSURE at TIME %f, ITER %d >> MAX: %.16g MIN: %.16g \n\n", time, iter, pressMax, pressMin); SOApprox_GetLimitValues(dsaturation, maxP, &dsatMax, &dsatMin); /* Prepare for next iteration: */ fprintf(stderr, "Distance change: %f to ", d); d = SOApprox_UpdateSolution(Anew, A); fprintf(stderr, " %f \n", d); norm = rn_norm(dim, A.el); /* Checks pressure convergence. */ fprintf(stderr, "RelTol: %f, AbsTol: %f, Norm: %f, pressMaxVar: %f \n", o.relTol, o.absTol, norm, pressMaxVar); if ((iter > 0) && ((d <= o.relTol*norm) && (d <= o.absTol))) { break; } if (iter > o.maxIter) { 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); for(i = 0; i < dim; i++)A.el[i] -= Aini.el[i]; PlotState( &(Mproperties.el[0]), pDim, tree, initP, finP, &o, &plotMeshDepth, "RUNNING", time, iter, outiter, &pressMax, &(plot_count[TOTPROP]), &(plotTotCount[TOTPROP]), &water_out, &oil_out, &init_oil_vol, &oil_vol, &oil_volA); for(i = 0; i < dim; i++)A.el[i] += Aini.el[i]; PlotState( &(Mproperties.el[1]), pDim, tree, initP, finP, &o, &plotMeshDepth, "RUNNING", time, iter, outiter, &satMax, &(plot_count[TOTPROP+1]), &(plotTotCount[TOTPROP+1]), &water_out, &oil_out, &init_oil_vol, &oil_vol, &oil_volA); PlotState( &(Mproperties.el[2]), pDim, tree, initP, finP, &o, &plotMeshDepth, "RUNNING", time, iter, outiter, &dsatMax, &(plot_count[TOTPROP+2]), &(plotTotCount[TOTPROP+2]), &water_out, &oil_out, &init_oil_vol, &oil_vol, &oil_volA); PlotState( &(Mproperties.el[3]), pDim, tree, initP, finP, &o, &plotMeshDepth, "RUNNING", time, iter, outiter, &qfluxMax, &(plot_count[TOTPROP+3]), &(plotTotCount[TOTPROP+3]), &water_out, &oil_out, &init_oil_vol, &oil_vol, &oil_volA); } iter++; } for(i = 0; i < dim; i++) B.el[i] = B.el[i] + dB.el[i] * o.tStep; SOApprox_GetLimitValues(saturation, maxP, &satMax, &satMin); fprintf(stderr, "SATURATION at TIME %f>> MAX: %.16g MIN: %.16g \n", time, satMax, satMin); oil_flux_out = 0.0; water_flux_out = 0.0; for(i = 0; i < dim; i++) Q.el[i] = 0.0; /* Producer well flux is adapted: */ for(i = 0; i < numwells; i++) if(wells.el[i].type == 1) { double sa, po, kra, kro, ma, mo, ba, bo, moba, mobo; /* Counts the flux and volume taken out. */ water_flux_out += - wells.el[i].wflux; water_out += - wells.el[i].wflux * o.tStep; oil_flux_out += - wells.el[i].oflux; oil_out += - wells.el[i].oflux * o.tStep; SOGrid_Index celli; moba = 0.0; mobo = 0.0; for(l = 0; l < max_cells; l++) { tf = (SOTentFunction *)bas.el[wells.el[i].ind]; celli = SOTent_get_brick(tf->d->index, pDim, tf->d->rank, l); SOGrid_cell_coords(celli, tf->d->rank, pDim, tent_min, tent_max); for(j = 0; j < pDim; j++) /* Canonical geometry. */ { tent_min[j] = tent_min[j] / pow(2.0, (double)j/pDim); tent_max[j] = tent_max[j] / pow(2.0, (double)j/pDim); } for(j = 0; j < pDim; j++) tent_max[j] = (tent_max[j] + tent_min[j]) / 2.0; saturation->m->eval(saturation, tent_max, &sa); pressure->m->eval(pressure, tent_max, &po); FMrelperm[WATER]->map(&sa, tent_max, &kra); FMrelperm[OIL]->map(&sa, tent_max, &kro); FMvisc[WATER]->map(&po, tent_max, &ma); FMvisc[OIL]->map(&po, tent_max, &mo); FMform[WATER]->map(&po, tent_max, &ba); FMform[OIL]->map(&po, tent_max, &bo); moba += (double)kra / (ma * ba) * (double)(wells.el[i].vol / max_cells); mobo += (double)kro / (mo * bo) * (double)(wells.el[i].vol / max_cells); } /* Resets the water and oil flux. */ if(moba == 0.0) wells.el[i].wflux = 0.0; else wells.el[i].wflux = (double) - wells.el[i].tflux / ((double) mobo / moba + 1.0); if(wells.el[i].wflux > 0.0) wells.el[i].wflux = 0.0; wells.el[i].oflux = - (wells.el[i].tflux + wells.el[i].wflux); if(wells.el[i].oflux > 0.0) wells.el[i].oflux = 0.0; //Q.el[wells.el[i].ind] = (double)wells.el[i].wflux / wells.el[i].vol; Q.el[wells.el[i].ind] = (double)(wells.el[i].wflux + wells.el[i].oflux) / wells.el[i].vol; fprintf(stderr, " DBG-> flux_OIL: %.16g, flux_WATER: %.16g, SOMA: %.16g \n\n", wells.el[i].oflux, wells.el[i].wflux, wells.el[i].oflux+wells.el[i].wflux); } else { water_vol_in += wells.el[i].wflux * o.tStep; water_in += wells.el[i].wflux * o.tStep; Q.el[wells.el[i].ind] = (double)(wells.el[i].wflux + wells.el[i].oflux) / wells.el[i].vol; } if(!((outiter-1) % plot_skip) && (outiter > 0))//(!(outiter % plot_skip)) { /* Gets the new oil total volume.*/ oil_volA = oil_vol; water_volA = water_vol; TotalOilVol( saturation, &NoMap, porosity, FMpor, maxtree, pDim, &tot_vol, &water_vol, &oil_vol, &reserv_box ); fprintf(stderr, " >>> DBG outiter: %d, oil_volA: %.16g, oil_vol: %.16g, water_volA: %.16g, water_vol: %.16g \n\n", outiter, oil_volA, oil_vol, water_volA, water_vol); 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); for(i = 0; i < dim; i++)A.el[i] -= Aini.el[i]; PlotState( &(Mproperties.el[0]), pDim, tree, initP, finP, &o, &plotMeshDepth, "TIME_STEP", time, iter, outiter, &pressMax, &(plot_count[TOTPROP]), &(plotTotCount[TOTPROP]), &water_out, &oil_out, &init_oil_vol, &oil_vol, &oil_volA); for(i = 0; i < dim; i++)A.el[i] += Aini.el[i]; PlotState( &(Mproperties.el[1]), pDim, tree, initP, finP, &o, &plotMeshDepth, "TIME_STEP", time, iter, outiter, &satMax, &(plot_count[TOTPROP+1]), &(plotTotCount[TOTPROP+1]), &water_out, &oil_out, &init_oil_vol, &oil_vol, &oil_volA); PlotState( &(Mproperties.el[2]), pDim, tree, initP, finP, &o, &plotMeshDepth, "TIME_STEP", time, iter, outiter, &dsatMax, &(plot_count[TOTPROP+2]), &(plotTotCount[TOTPROP+2]), &water_out, &oil_out, &init_oil_vol, &oil_vol, &oil_volA); PlotState( &(Mproperties.el[3]), pDim, tree, initP, finP, &o, &plotMeshDepth, "TIME_STEP", time, iter, outiter, &qfluxMax, &(plot_count[TOTPROP+3]), &(plotTotCount[TOTPROP+3]), &water_out, &oil_out, &init_oil_vol, &oil_vol, &oil_volA); /* Save result values. */ SaveResultStep( pDim, maxtree, pressure, reserv_box, uf, &NoMap, time+o.tStep, &schedule, &avg_press, &oil_flux, &water_flux, &oil_vol_out, &water_vol_out, &appr_oil_flux, &appr_water_flux, &appr_oil_vol_out, &appr_water_vol_out, water_flux_out, oil_flux_out, init_oil_vol, water_volA, water_vol, oil_volA, oil_vol, tot_vol, water_out, water_in, water_vol_in, oil_out, o.tStep, plot_skip ); water_vol_in = 0.0; } /* Dynamic basis adaptation - changes to the finest grid. */ if(o.changeBas && finecount == 0 && refinedcount == o.refinedskip) { DynChangeFuncBasis( pDim, tree, maxRank, maxRank, satMin, satMax, o.maximal, &new_tents, &bas, &dim, offset, evaldot, graddot, &Hev, &Hgr, pressure, Npressure, dsaturation, saturation, Qflux, porosity, properties, B, A, Aini, P, TotH, b, intersec ); DynChangeFlux( pDim, bas, res_side, tents, new_tents, wells, numwells, dim ); free((SOTent *)tents.el); tents.el = new_tents.el; tents.nel = new_tents.nel; finecount++; } /* Dynamic basis adaptation - changes back to the refined grid. */ if(o.changeBas && finecount == o.fineskip) { if( tree->root->c[LO] != NULL && tree->root->c[HI] != NULL ) { SOGrid_subtree_free(tree->root->c[LO]); tree->root->c[LO] = NULL; SOGrid_subtree_free(tree->root->c[HI]); tree->root->c[HI] = NULL; } ShatterNode (pDim, tree->root, 0, lo, hi, minRank, maxRank, refinerate, numwells, wells); DynChangeFuncBasis( pDim, tree, minRank, maxRank, satMin, satMax, o.maximal, &new_tents, &bas, &dim, offset, evaldot, graddot, &Hev, &Hgr, pressure, Npressure, dsaturation, saturation, Qflux, porosity, properties, B, A, Aini, P, TotH, b, intersec ); DynChangeFlux( pDim, bas, res_side, tents, new_tents, wells, numwells, dim ); free((SOTent *)tents.el); tents.el = new_tents.el; tents.nel = new_tents.nel; finecount = 0; refinedcount = 0; } if(finecount == 0) refinedcount++; else finecount++; outiter++; iter = 0; time += o.tStep; } end_time = today(); SaveResults( o.outName, init_time, end_time, schedule, avg_press, oil_flux, water_flux, oil_vol_out, water_vol_out, appr_oil_flux, appr_water_flux, appr_oil_vol_out, appr_water_vol_out ); for(i = 0; i < TOTPROP; i++) if(properties.el[i].plot) SOPlot_CloseStream(properties.el[i].ps); for(i = 0; i < NUMMAIN; i++) SOPlot_CloseStream(Mproperties.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, "SODynSimoleo \\\n"); PPUSAGE(pp, " -reservoirName \\\n"); PPUSAGE(pp, " -matName NAME \\\n"); PPUSAGE(pp, " -outName NAME \\\n"); PPUSAGE(pp, " -tFin TOTAL_TIME \\\n"); PPUSAGE(pp, " -tStep TIME_INTERVAL \\\n"); PPUSAGE(pp, " [ -maximal | -minimal ] \\\n"); PPUSAGE(pp, " [ -changeBas ] \\\n"); PPUSAGE(pp, " [ -refinedskip NUM ] \\\n"); PPUSAGE(pp, " [ -fineskip NUM ] \\\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, " [ -plotAll ] [ -plotFinal ] \\\n"); PPUSAGE(pp, SOPlotParams_FunctionHelp " \n"); SOParams_GetKeyword(pp, "-reservoirName"); o.reservoirName = SOParams_GetNext(pp); SOParams_GetKeyword(pp, "-matName"); o.matName = SOParams_GetNext(pp); SOParams_GetKeyword(pp, "-outName"); o.outName = SOParams_GetNext(pp); if (SOParams_KeywordPresent(pp, "-maximal")) { o.maximal = TRUE; } else if (SOParams_KeywordPresent(pp, "-minimal")) { o.maximal = FALSE; } else { SOParams_Error(pp, "must specify \"-maximal\" or \"-minimal\""); } o.changeBas = SOParams_KeywordPresent(pp, "-changeBas"); if (SOParams_KeywordPresent(pp, "-refinedskip")) { o.refinedskip = SOParams_GetNextInt(pp, 0, INT_MAX); } if (SOParams_KeywordPresent(pp, "-fineskip")) { o.fineskip = SOParams_GetNextInt(pp, 0, INT_MAX); } 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); o.tFin *= DAYSECONDS; SOParams_GetKeyword(pp, "-tStep"); o.tStep = SOParams_GetNextDouble(pp, 0.0, 1.0e8); o.tStep *= DAYSECONDS; /* 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 *res_side, int *minRank, int *maxRank, double *refinerate, double *initP, double *finP, double *init_press, double *init_sat, int *plot_skip, Resv_Func_vec *Mproperties, Resv_Func_vec *properties, int *numwells, Resv_Well_vec *wells ) { int i; char *tmpstr; filefmt_read_header(rd, "SOReservoir", SOReservoir_FileFormat); *pDim = nget_int(rd, "domain_dim"); fget_eol(rd); nget_name_eq(rd, "res_side"); *res_side = fget_double(rd); fget_eol(rd); *minRank = nget_int(rd, "minRank"); fget_eol(rd); *maxRank = nget_int(rd, "maxRank"); fget_eol(rd); *refinerate = nget_double(rd, "refinerate"); fget_eol(rd); nget_name_eq(rd, "plot_points"); for (i = 0; i < *pDim; i++){ initP[i] = fget_double(rd); }//initP[i] = initP[i] / pow(2.0, (double)i / *pDim); } for (i = 0; i < *pDim; i++){ finP[i] = fget_double(rd); }//finP[i] = finP[i] / pow(2.0, (double)i / *pDim); } 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, &(Mproperties->el[0]), "PRESS" ); ResvFuncRead( rd, &(Mproperties->el[1]), "SAT" ); ResvFuncRead( rd, &(Mproperties->el[2]), "DSAT" ); ResvFuncRead( rd, &(Mproperties->el[3]), "QFLUX" ); 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]), *pDim ); 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 pDim ) { int i; well->type = nget_int(rd, "type"); fget_eol(rd); nget_name_eq(rd, "position"); well->pos = malloc(pDim * sizeof(double)); assert(well->pos != NULL, "No memory for well->pos allocation."); for (i = 0; i < pDim; i++) { well->pos[i] = fget_double(rd); }//well->pos[i] = well->pos[i] / pow(2.0, (double)i/pDim); } fget_eol(rd); well->tflux = nget_double(rd, "flux"); } void SetPropFunctions ( Resv_Func_vec *Mproperties, Resv_Func_vec *properties, SOFunction *porosity, SOFunction *pressure, SOFunction *saturation, SOFunction *dsaturation, SOFunction *Qflux ) { /* ++++++++++++++++++++++++++++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 = porosity; /* POROSITY */ /*++++++++++++++++++++++++++++MAIN_PROPERTIES++++++++++++++++*/ Mproperties->el[0].F = pressure; /* PRESSURE */ Mproperties->el[1].F = saturation; /* SATURATION */ Mproperties->el[2].F = dsaturation; /* DIFFERENTIAL_SATURATION */ Mproperties->el[3].F = Qflux; /* FLUX BY VOLUME (QFLUX) */ } void OpenPlot ( Resv_Func *RF, Options *o, int *plotTotCount ) { char paperSize[50]; int pos = 0, hCount = 1, vCount = 1; strcpy(paperSize,""); 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, RF->ps_eps?0:7, 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, SQRTHALF}; printf("\n+++ Plot %s +++++++++++++++\n", RF->name); printf("iter: %d, outiter: %d, fMax: %.3g, 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:%.2f 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, !RF->ps_eps ? tmpstr : (char*)NULL ); (*plot_count)++; } if(!((*plot_count) % *plotTotCount)) { sprintf(tmpstr,"%08d", *plot_count); SOPlot_AddPage(RF->ps, tmpstr); } if(fMax == NULL) { sprintf(tmpstr,"%s - %s \n iter: %d - outiter: %d \n time:%.3g days \n 2D_PlotMax: %.3g 2D_PlotMin: %.3g \n", RF->name, tag, iter, outiter, time/DAYSECONDS, PlotMax, PlotMin); sprintf(tmpstr2," init oil vol: %.3g oil vol: %.3g \n water out: %.3g oil out: %.3g \n vol_diff: %.3g oil_diff: %.3g", *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:%.3g days fMax:%.3g \n 2D_PlotMax: %.3g 2D_PlotMin: %.3g \n", RF->name, tag, iter, outiter, time/DAYSECONDS, *fMax, PlotMax, PlotMin); sprintf(tmpstr2," init oil vol: %.3g oil vol: %.3g \n water out: %.3g oil out: %.3g \n vol_diff: %.3g oil_diff: %.3g", *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; double limit_val = RF->fRange; //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), +(limit_val), //-(limit_val * (double)10/8), +(limit_val * (double)10/8), limit_val, PLOTDIVS, 0.20, &PlotMin, &PlotMax, !RF->ps_eps ? tmpstr : (char*)NULL ); (*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", *plot_count); SOPlot_AddPage(RF->ps, tmpstr); } fprintf(stderr, "\n++++++++++++++++++++++++\n"); } void ShatterNode ( SOGrid_Dim pDim, SOGrid_Node *p, SOGrid_Rank r, double *lo, double *hi, int minRank, int maxRank, double refinerate, int numwells, Resv_Well_vec wells ) { int i, j; bool split; double d2, lowd2, sz, szi, ci, lo_child[pDim], hi_child[pDim]; double dist = 1.0, idist; lowd2 = 1.0; sz = hi[0] - lo[0]; for(i = 0; i < pDim; i++) { szi = hi[i] - lo[i]; if (szi > sz) { sz = szi; } } for(j = 0; j < numwells; j++) { d2 = 0.0; for(i = 0; i < pDim; i++) { ci = (double)(lo[i] + hi[i])/2.0 - wells.el[j].pos[i]; d2 += ci*ci; } //if(d2 < lowd2) lowd2 = (wells.el[j].type?d2:d2*0.6); if(d2 < lowd2) lowd2 = d2; //if(d2 < lowd2 && wells.el[j].type == 1) lowd2 = d2; } split = (sz > (double) lowd2 * refinerate); for(i = 0; i < pDim; i++) { if(lo[i] < 0.5) idist = lo[i]; if(hi[i] > 0.5) idist = ((double)1.0/pow(2.0, (double)i/pDim)) - hi[i]; if(idist < dist) dist = idist; } split = (split || dist < (double)1.0 / (1 << (maxRank / pDim))); assert((maxRank >= minRank), "Wrong values for grid depths: minRank > maxRank."); if (r < minRank) split = TRUE; if (r >= maxRank) split = FALSE; if (split) { SOGrid_Node_split(p); /* Computes low and high coordinates of low child. */ for (i = 0; i < pDim; i++) { lo_child[i] = lo[i]; hi_child[i] = (i == (r%pDim) ? (lo[i]+hi[i])/2 : hi[i]); } ShatterNode(pDim, p->c[LO], r+1, lo_child, hi_child, minRank, maxRank, refinerate, numwells, wells); /* Compute low and high coordinates of high child. */ for (i = 0; i < pDim; i++) { lo_child[i] = (i == (r%pDim) ? (lo[i]+hi[i])/2 : lo[i]); hi_child[i] = hi[i]; } ShatterNode(pDim, p->c[HI], r+1, lo_child, hi_child, minRank, maxRank, refinerate, numwells, wells); } } void DynShatterNode ( SOGrid_Dim pDim, SOGrid_Node *p, SOGrid_Rank r, double *lo, double *hi, //double *maxhi, int minRank, int maxRank, double diff_limit, SOFunction *f ) { bool split;//, border = FALSE; int i, j; double mid_val, maxdiff=0, diff, lo_child[pDim], hi_child[pDim]; //, abs_mid_val; double pnt_mid[pDim], pnt_lo[pDim], pnt_hi[pDim], pnt_lo_val, pnt_hi_val; for(i = 0; i < pDim; i++) pnt_mid[i] = (double) (hi[i] + lo[i]) / 2.0; f->m->eval(f, pnt_mid, &mid_val); for(j = 0; j < pDim; j++) { for(i = 0; i < pDim; i++) if(i != j){ pnt_lo[i] = pnt_mid[i]; pnt_hi[i] = pnt_mid[i]; } pnt_lo[j] = lo[j]; f->m->eval(f, pnt_lo, &pnt_lo_val); pnt_hi[j] = hi[j]; f->m->eval(f, pnt_hi, &pnt_hi_val); diff = fabs(mid_val - (pnt_lo_val + pnt_hi_val)/2.0); if(diff > maxdiff) maxdiff = diff; /* if(lo[j] == 0.0 || hi[j] == maxhi[j]){ border = TRUE; abs_mid_val = fabs(mid_val); } */ } split = (maxdiff > diff_limit);// || (border && abs_mid_val > 2.0 * diff_limit)); if (r < minRank){split = TRUE;} if (r >= maxRank){split = FALSE;} if (split) { if( p->c[LO] == NULL && p->c[HI] == NULL )SOGrid_Node_split(p); /* Compute low and high coordinates of low child. */ for (i = 0; i < pDim; i++) { lo_child[i] = lo[i]; hi_child[i] = (i == (r%pDim) ? (lo[i]+hi[i])/2 : hi[i]); } DynShatterNode(pDim, p->c[LO], r+1, lo_child, hi_child, minRank, maxRank, diff_limit, f); /* ShatterNode(pDim, p->c[LO], r+1, lo_child, hi_child, maxhi, minRank, maxRank, diff_limit, f); */ /* Compute low and high coordinates of high child. */ for (i = 0; i < pDim; i++) { lo_child[i] = (i == (r%pDim) ? (lo[i]+hi[i])/2 : lo[i]); hi_child[i] = hi[i]; } DynShatterNode(pDim, p->c[HI], r+1, lo_child, hi_child, minRank, maxRank, diff_limit, f); /* ShatterNode(pDim, p->c[HI], r+1, lo_child, hi_child, maxhi, minRank, maxRank, diff_limit, f); */ } /* if(unite) */ /* { */ /* if( p->c[LO] != NULL && p->c[HI] != NULL ) */ /* { SOGrid_subtree_free(p->c[LO]); p->c[LO] = NULL; */ /* SOGrid_subtree_free(p->c[HI]); p->c[HI] = NULL; */ /* } */ /* } */ } void DynShatterNode2 ( SOGrid_Dim pDim, SOGrid_Node *p, SOGrid_Rank r, double *lo, double *hi, int minRank, int maxRank ) { int i; double lo_child[pDim], hi_child[pDim]; if(r < maxRank) { if( p->c[LO] == NULL && p->c[HI] == NULL ) { SOGrid_Node_split(p); if(r < (maxRank-1)) { SOGrid_Node_split(p->c[LO]); SOGrid_Node_split(p->c[HI]); } } else { /* Compute low and high coordinates of low child. */ for (i = 0; i < pDim; i++) { lo_child[i] = lo[i]; hi_child[i] = (i == (r%pDim) ? (lo[i]+hi[i])/2 : hi[i]); } DynShatterNode2(pDim, p->c[LO], r+1, lo_child, hi_child, minRank, maxRank); /* Compute low and high coordinates of high child. */ for (i = 0; i < pDim; i++) { lo_child[i] = (i == (r%pDim) ? (lo[i]+hi[i])/2 : lo[i]); hi_child[i] = hi[i]; } DynShatterNode2(pDim, p->c[HI], r+1, lo_child, hi_child, minRank, maxRank); } } } void SetGenDotMatrices ( ListMatEntry_vec *evaldot, ListMatEntry_vec *graddot, SOListMatrix Meval, SOListMatrix Mgrad, unsigned int offset, unsigned int maxindices ) { unsigned int i, j, index, indexi; for(i = 0; i < maxindices; i++) { evaldot[i] = ListMatEntry_vec_new(Meval.cols); evaldot[i].nel = 0; graddot[i] = ListMatEntry_vec_new(Mgrad.cols); graddot[i].nel = 0; } i = 0; while(i < Meval.ents.nel) { index = Meval.ents.el[i].row; while(Meval.ents.el[i].row == index) { indexi = Meval.ents.el[i].row - offset; evaldot[indexi].el[evaldot[indexi].nel].row = Meval.ents.el[i].row; evaldot[indexi].el[evaldot[indexi].nel].col = Meval.ents.el[i].col; evaldot[indexi].el[evaldot[indexi].nel].va = (double *)malloc(sizeof(double)*Meval.length); for(j = 0; j < Meval.length; j++) evaldot[indexi].el[evaldot[indexi].nel].va[j] = Meval.ents.el[i].va[j]; evaldot[indexi].nel++; i++; } } i = 0; while(i < Mgrad.ents.nel) { index = Mgrad.ents.el[i].row; while(Mgrad.ents.el[i].row == index) { indexi = Mgrad.ents.el[i].row - offset; graddot[indexi].el[graddot[indexi].nel].row = Mgrad.ents.el[i].row; graddot[indexi].el[graddot[indexi].nel].col = Mgrad.ents.el[i].col; graddot[indexi].el[graddot[indexi].nel].va = (double *)malloc(sizeof(double)*Mgrad.length); for(j = 0; j < Mgrad.length; j++) graddot[indexi].el[graddot[indexi].nel].va[j] = Mgrad.ents.el[i].va[j]; graddot[indexi].nel++; i++; } } } void SetDotMatrices ( unsigned int dim, unsigned int offset, Basis bas, ListMatEntry_vec *evaldot, ListMatEntry_vec *graddot, ListMatEntry_vec *Hev_ents, ListMatEntry_vec *Hgr_ents ) { int i, j, k, hi_ind, low_ind; SOTentFunction *itf, *jtf; Hev_ents->nel = 0; Hgr_ents->nel = 0; for(i = 0; i < dim; i++) for(j = 0; j < dim; j++) { itf = (SOTentFunction *)bas.el[i]; jtf = (SOTentFunction *)bas.el[j]; if(itf->d->index > jtf->d->index){ hi_ind = itf->d->index; low_ind = jtf->d->index; } else { hi_ind = jtf->d->index; low_ind = itf->d->index; } k = 0; while(k < evaldot[hi_ind - offset].nel) { if(evaldot[hi_ind - offset].el[k].col == low_ind) { Hev_ents->el[Hev_ents->nel].row = i; Hev_ents->el[Hev_ents->nel].col = j; Hev_ents->el[Hev_ents->nel].va = (double *)evaldot[hi_ind - offset].el[k].va; Hev_ents->nel++; } k++; } k = 0; while(k < graddot[hi_ind - offset].nel) { if(graddot[hi_ind - offset].el[k].col == low_ind) { Hgr_ents->el[Hgr_ents->nel].row = i; Hgr_ents->el[Hgr_ents->nel].col = j; Hgr_ents->el[Hgr_ents->nel].va = (double *)graddot[hi_ind - offset].el[k].va; Hgr_ents->nel++; } k++; } } } void ChangeBasis ( unsigned int offset, unsigned int dim, unsigned int new_dim, Basis bas, Basis new_bas, SOListMatrix Hev, ListMatEntry_vec *evaldot, double *TotH, double_vec b, double_vec coeffs ) { int i, j, l, k, hi_ind, low_ind, max_cells; SOTentFunction *itf, *jtf; jtf = (SOTentFunction *)bas.el[0]; max_cells = 1 << (jtf->d->fn.pDim); for(i = 0; i < new_dim; i++) { itf = (SOTentFunction *)new_bas.el[i]; b.el[i] = 0.0; for(j = 0; j < dim; j++) { jtf = (SOTentFunction *)bas.el[j]; if(itf->d->index > jtf->d->index){ hi_ind = itf->d->index; low_ind = jtf->d->index; } else { hi_ind = jtf->d->index; low_ind = itf->d->index; } k = 0; while(k < evaldot[hi_ind - offset].nel) { if(evaldot[hi_ind - offset].el[k].col == low_ind) for(l = 0; l < max_cells; l++) b.el[i] += evaldot[hi_ind - offset].el[k].va[l] * coeffs.el[j]; k++; } } } SOApprox_ListGaussElim(TotH, Hev, b, coeffs, FALSE); } void SaveResultStep ( SOGrid_Dim pDim, SOGrid_Tree *tree, SOFunction *pressure, double reserv_box, SOFunction *uf, FuncMap *NoMap, double time, int_vec *schedule, double_vec *avg_press, double_vec *oil_flux, double_vec *water_flux, double_vec *oil_vol_out, double_vec *water_vol_out, double_vec *appr_oil_flux, double_vec *appr_water_flux, double_vec *appr_oil_vol_out, double_vec *appr_water_vol_out, double water_flux_out, double oil_flux_out, double init_oil_vol, double water_volA, double water_vol, double oil_volA, double oil_vol, double tot_vol, double water_out, double water_in, double water_vol_in, double oil_out, double tStep, int plot_skip ) { double total_press = 0; schedule->el[schedule->nel] = (int)time/DAYSECONDS; schedule->nel++; /* Gets the total pressure. */ SOFunction_RawDot( pressure, NoMap, uf, NoMap, NULL, tree, &total_press ); printf("################ RESULTS day %d ########################### \n", schedule->el[schedule->nel - 1]); avg_press->el[avg_press->nel] = (double)total_press * reserv_box / tot_vol; avg_press->nel++; printf(" AVG_PRESS on day %d -> %.16g \n", schedule->el[schedule->nel - 1], avg_press->el[avg_press->nel - 1]); oil_flux->el[oil_flux->nel] = (double)(oil_volA - oil_vol) / (tStep * plot_skip); oil_flux->nel++; printf(" OIL_FLUX on day %d -> %.16g \n", schedule->el[schedule->nel - 1], oil_flux->el[oil_flux->nel - 1]); appr_oil_flux->el[appr_oil_flux->nel] = oil_flux_out; appr_oil_flux->nel++; printf(" APPR_OIL_FLUX on day %d -> %.16g \n", schedule->el[schedule->nel - 1], appr_oil_flux->el[appr_oil_flux->nel - 1]); water_flux->el[water_flux->nel] = (double)(water_vol_in + water_volA - water_vol) / (tStep * plot_skip); water_flux->nel++; printf(" WATER_FLUX on day %d -> %.16g \n", schedule->el[schedule->nel - 1], water_flux->el[water_flux->nel - 1]); appr_water_flux->el[appr_water_flux->nel] = water_flux_out; appr_water_flux->nel++; printf(" APPR_WATER_FLUX on day %d -> %.16g \n", schedule->el[schedule->nel - 1], appr_water_flux->el[appr_water_flux->nel - 1]); oil_vol_out->el[oil_vol_out->nel] = init_oil_vol - oil_vol; oil_vol_out->nel++; printf(" OIL_OUT on day %d -> %.16g \n", schedule->el[schedule->nel - 1], oil_vol_out->el[oil_vol_out->nel - 1]); appr_oil_vol_out->el[appr_oil_vol_out->nel] = oil_out; appr_oil_vol_out->nel++; printf(" APPR_OIL_OUT on day %d -> %.16g \n", schedule->el[schedule->nel - 1], appr_oil_vol_out->el[appr_oil_vol_out->nel - 1]); water_vol_out->el[water_vol_out->nel] = water_in + (tot_vol - init_oil_vol) - water_vol; water_vol_out->nel++; printf(" WATER_OUT on day %d -> %.16g \n", schedule->el[schedule->nel - 1], water_vol_out->el[water_vol_out->nel - 1]); appr_water_vol_out->el[appr_water_vol_out->nel] = water_out; appr_water_vol_out->nel++; printf(" APPR_WATER_OUT on day %d -> %.16g \n", schedule->el[schedule->nel - 1], appr_water_vol_out->el[appr_water_vol_out->nel - 1]); printf("########################################### \n\n"); } void SaveResults ( char *outName, char *init_time, char *end_time, int_vec schedule, double_vec avg_press, double_vec oil_flux, double_vec water_flux, double_vec oil_vol_out, double_vec water_vol_out, double_vec appr_oil_flux, double_vec appr_water_flux, double_vec appr_oil_vol_out, double_vec appr_water_vol_out ) { FILE *wr = open_write(txtcat(outName, ".log")); fprintf(wr, "## Results from simulation %s\n\n", outName); fprintf(wr, "Started at: %s\n", init_time); fprintf(wr, "Finished at: %s\n\n", end_time); WriteArray(wr, "AVG_PRESS", schedule, avg_press); WriteArray(wr, "OIL_FLUX", schedule, oil_flux); WriteArray(wr, "APPR_OIL_FLUX", schedule, appr_oil_flux); WriteArray(wr, "WATER_FLUX", schedule, water_flux); WriteArray(wr, "APPR_WATER_FLUX", schedule, appr_water_flux); WriteArray(wr, "OIL_OUT", schedule, oil_vol_out); WriteArray(wr, "APPR_OIL_OUT", schedule, appr_oil_vol_out); WriteArray(wr, "WATER_OUT", schedule, water_vol_out); WriteArray(wr, "APPR_WATER_OUT", schedule, appr_water_vol_out); fclose(wr); SO1DPlot_Array(avg_press, txtcat("AVG_PRESS", outName), 12000000.0, 14000000.0, 10.0); SO1DPlot_Array(oil_flux, txtcat("OIL_FLUX", outName), 0.006, 0.010, 10.0); SO1DPlot_Array(appr_oil_flux, txtcat("APPR_OIL_FLUX", outName), 0.006, 0.010, 10.0); SO1DPlot_Array(water_flux, txtcat("WATER_FLUX", outName), -0.04, 0.0008, 10.0); SO1DPlot_Array(appr_water_flux, txtcat("APPR_WATER_FLUX", outName), -0.04, 0.0008, 10.0); SO1DPlot_Array(oil_vol_out, txtcat("OIL_OUT", outName), 2000.0, 750000.0, 10.0); SO1DPlot_Array(appr_oil_vol_out, txtcat("APPR_OIL_OUT", outName), 2000.0, 750000.0, 10.0); SO1DPlot_Array(water_vol_out, txtcat("WATER_OUT", outName), 100.0, 9500.0, 10.0); SO1DPlot_Array(appr_water_vol_out, txtcat("APPR_WATER_OUT", outName), 100.0, 9500.0, 10.0); } void WriteArray( FILE *wr, char *descr, int_vec schedule, double_vec array ) { int i; fprintf(wr, ">> %s:\n\n", descr); for(i = 0; i < array.nel; i++) fprintf(wr, " day: %d, value: %.16g \n", schedule.el[i], array.el[i]); fprintf(wr, "\n\n"); } void DynChangeFuncBasis ( SOGrid_Dim pDim, SOGrid_Tree *tree, double minRank, double maxRank, double satMin, double satMax, bool maximal, SOTent_vec *tents, Basis *bas, unsigned int *dim, unsigned int offset, ListMatEntry_vec *evaldot, ListMatEntry_vec *graddot, SOListMatrix *Hev, SOListMatrix *Hgr, SOFunction *pressure, SOFunction *Npressure, SOFunction *dsaturation, SOFunction *saturation, SOFunction *Qflux, SOFunction *porosity, Resv_Func_vec properties, double_vec B, double_vec A, double_vec Aini, double_vec P, double *TotH, double_vec b, int_vec *intersec ) { int i; double lo[pDim], hi[pDim]; SOTent_vec new_tents; Basis new_bas; unsigned int new_dim; SOLinCombFunction *lcf; for(i = 0; i < pDim; i++){ lo[i] = 0.0; hi[i] = (double)1.0/pow(2.0, (double)i/pDim); } if(minRank == maxRank) DynShatterNode2( pDim, tree->root, 0, lo, hi, minRank, maxRank ); else DynShatterNode( pDim, tree->root, 0, lo, hi, minRank, maxRank, DIFF_LIMIT, saturation ); /* Gets new basis and reservoir domain dimension */ if (maximal) new_tents = SOTent_maximal_basis(tree); else new_tents = SOTent_minimal_basis(tree); new_bas = SOFunctionRef_vec_new(new_tents.nel); new_dim = new_bas.nel; printf(" DYN BASIS CHANGE>> new_bas.nel: %d \n\n\n", new_bas.nel); for (i = 0; i < new_tents.nel; i++) { SOTent t = new_tents.el[i]; new_bas.el[i] = (SOFunction *)SOTentFunction_Make(pDim, t.cx, t.r); } /* Gets the new dot product values for the main matrices and . */ SetDotMatrices( new_dim, offset, new_bas, evaldot, graddot, &(Hev->ents), &(Hgr->ents) ); Hev->rows = new_dim; Hev->cols = new_dim; Hgr->rows = new_dim; Hgr->cols = new_dim; for(i = 0; i < new_dim; i++) intersec[i].nel = 0; for(i = 0; i < Hev->ents.nel; i++) { intersec[Hev->ents.el[i].row].el[intersec[Hev->ents.el[i].row].nel] = Hev->ents.el[i].col; intersec[Hev->ents.el[i].row].nel++; } // for(i = 0; i < *dim; i++)printf(" DEBUG>>> B.el[%d]: %.16g \n\n", i, B.el[i]); ChangeBasis( offset, *dim, new_dim, *bas, new_bas, *Hev, evaldot, TotH, b, B ); // for(i = 0; i < new_dim; i++)printf(" DEBUG NEW>>> B.el[%d]: %.16g \n\n", i, B.el[i]); // for(i = 0; i < *dim; i++)printf(" DEBUG>>> A.el[%d]: %.16g \n\n", i, A.el[i]); ChangeBasis( offset, *dim, new_dim, *bas, new_bas, *Hev, evaldot, TotH, b, A ); // for(i = 0; i < new_dim; i++)printf(" DEBUG NEW>>> A.el[%d]: %.16g \n\n", i, A.el[i]); // for(i = 0; i < *dim; i++)printf(" DEBUG>>> Aini.el[%d]: %.16g \n\n", i, Aini.el[i]); ChangeBasis( offset, *dim, new_dim, *bas, new_bas, *Hev, evaldot, TotH, b, Aini ); // for(i = 0; i < new_dim; i++)printf(" DEBUG NEW>>> Aini.el[%d]: %.16g \n\n", i, Aini.el[i]); if(strcmp(properties.el[TOTPROP - 1].descr, "Fpor_NoMap") == 0) ChangeBasis( offset, *dim, new_dim, *bas, new_bas, *Hev, evaldot, TotH, b, P ); lcf = SOLinCombFunction_Cast(pressure); lcf->d->bas = new_bas; lcf = SOLinCombFunction_Cast(Npressure); lcf->d->bas = new_bas; lcf = SOLinCombFunction_Cast(saturation); lcf->d->bas = new_bas; lcf = SOLinCombFunction_Cast(dsaturation); lcf->d->bas = new_bas; lcf = SOLinCombFunction_Cast(Qflux); lcf->d->bas = new_bas; lcf = SOLinCombFunction_Cast(porosity); lcf->d->bas = new_bas; free((SOFunctionRef *)bas->el); *dim = new_dim; *tents = new_tents; *bas = new_bas; } void DynChangeFlux ( SOGrid_Dim pDim, Basis bas, double res_side, SOTent_vec tents, SOTent_vec new_tents, Resv_Well_vec wells, int numwells, unsigned int new_dim ) { int i, j, new_ind[numwells]; for(i = 0; i < numwells; i++) { new_ind[i] = -1; printf(" DEBUG WELL CGH 1>>>> bas ind: %d , ind: %d \n\n", (int)wells.el[i].ind, (int)tents.el[wells.el[i].ind].cx); for(j = 0; j < new_dim; j++) if(tents.el[wells.el[i].ind].cx == new_tents.el[j].cx) new_ind[i] = j; wells.el[i].ind = new_ind[i]; if(new_ind[i] == -1) { SOApprox_GetBasisIndeces(bas, pDim, 1, &(wells.el[i].ind), wells.el[i].pos); SOTentFunction *tf = (SOTentFunction *)bas.el[wells.el[i].ind]; double tent_min[pDim], tent_max[pDim]; SOGrid_cell_coords(tf->d->index, tf->d->rank, pDim, tent_min, tent_max); for(j = 0; j < pDim; j++) /* Canonical geometry. */ { tent_min[j] = tent_min[j] / pow(2.0, (double)j/pDim); tent_max[j] = tent_max[j] / pow(2.0, (double)j/pDim); } wells.el[i].vol = 1.0; for(j = 0; j < pDim; j++) wells.el[i].vol *= 2.0 * (tent_max[j] - tent_min[j]) * res_side; wells.el[i].vol = (double)wells.el[i].vol / (1 << pDim); wells.el[i].qflux = (double)wells.el[i].tflux / wells.el[i].vol; } printf(" DEBUG WELL CGH 2>>>> bas ind: %d , ind: %d \n\n", (int)wells.el[i].ind, (int)new_tents.el[wells.el[i].ind].cx); } } void TestMatrixHSymmetry(SOMatrix H, unsigned int dim) { int i, j; printf(" \n\n ****** TestMatrixHSymmetry ...\n\n"); for(i = 0; i < H.ents.nel; i++) if(H.ents.el[i].row < dim) for(j = 0; j < H.ents.nel; j++) { if(H.ents.el[i].col < dim && H.ents.el[j].row == H.ents.el[i].col && H.ents.el[j].col == H.ents.el[i].row) if(H.ents.el[i].va != H.ents.el[j].va) printf(" H Symm>> ERROR!! -> %.16g != -> %.16g \n\n", H.ents.el[i].row, H.ents.el[i].col, H.ents.el[i].va, H.ents.el[j].row, H.ents.el[j].col, H.ents.el[j].va); if(H.ents.el[j].row == H.ents.el[i].col-dim && H.ents.el[j].col == H.ents.el[i].row+dim) if(H.ents.el[i].va != H.ents.el[j].va) printf(" H Symm>> ERROR!! -> %.16g != -> %.16g \n\n", H.ents.el[i].row, H.ents.el[i].col, H.ents.el[i].va, H.ents.el[j].row, H.ents.el[j].col, H.ents.el[j].va); } else for(j = 0; j < H.ents.nel; j++) { if(H.ents.el[j].row == H.ents.el[i].col+dim && H.ents.el[j].col == H.ents.el[i].row-dim) if(H.ents.el[i].va != H.ents.el[j].va) printf(" H Symm>> ERROR!! -> %.16g != -> %.16g \n\n", H.ents.el[i].row, H.ents.el[i].col, H.ents.el[i].va, H.ents.el[j].row, H.ents.el[j].col, H.ents.el[j].va); if(H.ents.el[j].row >= dim && H.ents.el[j].row == H.ents.el[i].col && H.ents.el[j].col == H.ents.el[i].row) if(H.ents.el[i].va != H.ents.el[j].va) printf(" H Symm>> ERROR!! -> %.16g != -> %.16g \n\n", H.ents.el[i].row, H.ents.el[i].col, H.ents.el[i].va, H.ents.el[j].row, H.ents.el[j].col, H.ents.el[j].va); } for(i = 0; i < H.ents.nel; i++) for(j = 0; j < H.ents.nel; j++) if(H.ents.el[i].row < (int)dim/2 && H.ents.el[i].col > (int)dim/2 && H.ents.el[i].col < dim) if(H.ents.el[j].row == dim - H.ents.el[i].row - 1 && H.ents.el[j].col == dim - H.ents.el[i].col - 1) if(fabs(H.ents.el[i].va - H.ents.el[j].va) > 10e-9) printf(" H Symm>> ERROR 2!! -> %.16g != -> %.16g \n\n", H.ents.el[i].row, H.ents.el[i].col, H.ents.el[i].va, H.ents.el[j].row, H.ents.el[j].col, H.ents.el[j].va); for(i = 0; i < H.ents.nel; i++) for(j = 0; j < H.ents.nel; j++) if(H.ents.el[i].row < (int)dim/2 && (H.ents.el[i].col-dim) > (int)dim/2 && (H.ents.el[i].col-dim) < dim) if(H.ents.el[j].row == dim - H.ents.el[i].row - 1 && (H.ents.el[j].col-dim) == dim - (H.ents.el[i].col-dim) - 1) if(fabs(H.ents.el[i].va - H.ents.el[j].va) > 10e-9) printf(" H Symm>> ERROR 3!! -> %.16g != -> %.16g \n\n", H.ents.el[i].row, H.ents.el[i].col, H.ents.el[i].va, H.ents.el[j].row, H.ents.el[j].col, H.ents.el[j].va); for(i = 0; i < H.ents.nel; i++) for(j = 0; j < H.ents.nel; j++) if((H.ents.el[i].row-dim) < (int)dim/2 && H.ents.el[i].col > (int)dim/2 && H.ents.el[i].col < dim) if((H.ents.el[j].row-dim) == dim - (H.ents.el[i].row-dim) - 1 && H.ents.el[j].col == dim - H.ents.el[i].col - 1) if(fabs(H.ents.el[i].va - H.ents.el[j].va) > 10e-9) printf(" H Symm>> ERROR 4!! -> %.16g != -> %.16g \n\n", H.ents.el[i].row, H.ents.el[i].col, H.ents.el[i].va, H.ents.el[j].row, H.ents.el[j].col, H.ents.el[j].va); for(i = 0; i < H.ents.nel; i++) for(j = 0; j < H.ents.nel; j++) if((H.ents.el[i].row-dim) < (int)dim/2 && (H.ents.el[i].col-dim) > (int)dim/2 && (H.ents.el[i].col-dim) < dim) if((H.ents.el[j].row-dim) == dim - (H.ents.el[i].row-dim) - 1 && (H.ents.el[j].col-dim) == dim - (H.ents.el[i].col-dim) - 1) if(fabs(H.ents.el[i].va - H.ents.el[j].va) > 10e-9) printf(" H Symm>> ERROR 5!! -> %.16g != -> %.16g \n\n", H.ents.el[i].row, H.ents.el[i].col, H.ents.el[i].va, H.ents.el[j].row, H.ents.el[j].col, H.ents.el[j].va); } void TestVecbSymmetry(double_vec b, unsigned int dim) { int i; printf(" \n\n ****** TestVecbSymmetry ...\n\n"); for(i = 0; i < (int)dim/2; i++) if(fabs(b.el[i] - b.el[dim - i - 1]) > 10e-8 || fabs(b.el[dim + i] - b.el[dim + dim - i - 1]) > 10e-8) printf(" b Symm>> ERROR!! -> <%.8g, %.8g> -> <%.8g, %.8g> \n", i, dim-i-1, b.el[i], b.el[dim-i-1], dim+i, dim+dim-i-1, b.el[dim+i], b.el[dim+dim-i-1]); }