/* Copyright 2010 Statoil ASA. This file is part of The Open Porous Media project (OPM). OPM is free software: you can redistribute it and/or modify it under the terms of the GNU General Public License as published by the Free Software Foundation, either version 3 of the License, or (at your option) any later version. OPM is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details. You should have received a copy of the GNU General Public License along with OPM. If not, see . */ /** @file upscale_relperm.C @brief Upscales relative permeability as a fuction of water saturation assuming capillary equilibrium. Description: Reads in a lithofacies geometry in Eclipse format, reads in J(S_w) and relpermcurve(S_w) for each stone type, and calculates upscaled (three directions) relative permeability curves as a function of Sw. The relative permeability computation is based on - Capillary equilibrium, p_c is spatially invariant. - Optional gravitational effects. If gravity is not specified, gravity will be assumed to be zero. Units handling: - Assumes cornerpoint file reports lengths in cm. - Input surface tension is in dynes/cm - Input density is in g/cm^3 - The denominator \sigma * cos(\phi) in J-function scaling is what we call "surface tension". If angle dependency is to be included, calculate the "surface tension" yourself. - Outputted capillary pressure is in Pascals. Steps in the code: 1: Process command line options. 2: Read Eclipse file 3: Read relperm- and J-function for each stone-type. 4: Tesselate the grid (Sintef code) 5: Find minimum and maximum capillary pressure from the J-functions in each cell. 6: Upscale water saturation as a function of capillary pressure 7: Upscale single phase permeability. 8: Upscale phase permeability for capillary pressures that corresponds to a uniform saturation grid, and compute relative permeability. 9: Print output to screen and optionally to file. */ #include #include #include #if DUNE_VERSION_NEWER(DUNE_COMMON, 2, 3) #include #else #include #endif #include #include #include #include #include #include // for DBL_MAX/DBL_MIN #include #include #include #include #include #include #include #include #include #include #include using namespace Opm; using namespace std; static void usage() { cerr << "Usage: upscale_relperm stoneA.txt stoneB.txt ..." << endl << "where the options are:" << endl << " -bc -- which boundary conditions to use." << endl << " Possible values are f (fixed), l (linear)" << endl << " and p (periodic). Default f (fixed)." << endl << " -points -- Number of saturation points to upscale for." << endl << " Uniformly distributed within saturation endpoints." << endl << " Default 30." << endl << " -relPermCurve -- For isotropic input, the column number in the stone-files" << endl << " that represents the phase to be upscaled," << endl << " typically 2 (default) for water and 3 for oil." << endl << " -jFunctionCurve -- the column number in the stone-files that" << endl << " represent the Leverett J-function. Default 4." << endl << " -upscaleBothPhases -- If this is true, relperm for both phases will be upscaled" << endl << " and both will be outputted to Eclipse format. Default true." << endl << " For isotropic input, relPermCurves is assumed to be 2 and 3," << endl << " for anisotropic input, relPermCurves are assumed to be 3-5" << endl << " and 6-8 respectively for the two phases" << endl << " -gravity -- use 9.81 for standard gravity. Default zero. Unit m/s^2." << endl << " -surfaceTension -- Surface tension to use in J-function/Pc conversion." << endl << " Default 11 dynes/cm (oil-water systems). In absence of" << endl << " a correct value, the surface tension for gas-oil systems " << endl << " could be 22.5 dynes/cm." << endl << " -waterDensity -- density of water, only applicable to non-zero" << endl << " gravity, g/cm³. Default 1" << endl << " -oilDensity -- density of oil, only applicable to non-zero" << endl << " gravity, g/cm³. Default 0.6" << endl << " -output -- filename for where to write upscaled values." << endl << " If not supplied, output will only go to " << endl << " the terminal (standard out)." << endl << " -interpolate -- If supplied, the output data points will be" << endl << " interpolated using monotone cubic interpolation" << endl << " on a uniform grid with the specified number of" << endl << " points. Suggested value: 1000." << endl << " -maxPermContrast -- maximal permeability contrast in model." << endl << " Default 10^7" << endl << " -minPerm -- Minimum floating point value allowed for" << endl << " phase permeability in computations. If set to zero," << endl << " some models can end up singular. Default 10^-12" << endl << " -maxPerm -- Maximum floating point value allowed for" << endl << " permeability. " << endl << " Default 100000. Unit Millidarcy." << endl << " -fluids -- Either ow for oil/water systems or go for gas/oil systems. Default ow." << endl << " In case of go, the waterDensity option should be set to gas density" << endl << " Also remember to specify the correct surface tension" << endl << " -krowxswirr -- Oil relative permeability in x-direction at Swirr(from SWOF table)." << endl << " In case of oil/gas, this value is needed to ensure consistensy" << endl << " between SWOF and SGOF tables. Only has affect if fluids is set to go" << endl << " and upscaleBothPhases is true." << endl << " If not set, the point is not inserted into the final table." << endl << " -krowyswirr -- Oil relative permeability in y-direction at Swirr(from SWOF table). See krowxswirr." << endl << " -krowzswirr -- Oil relative permeability in z-direction at Swirr(from SWOF table). See krowxswirr." << endl << " -doEclipseCheck -- Default true. Check that input relperm curves includes relperms at critical" << endl << " saturation points, i.e. that krw(swcrit)=0 and krow(swmax) = 0 and similar for oil/gas." << endl << " -critRelpermThresh -- If minimum relperm values are less than this threshold, they are set to zero" << endl << " and will pass the EclipseCheck. Default 10^-6" << endl << "If only one stone-file is supplied, it is used for all stone-types defined" << endl << "in the geometry. If more than one, it corresponds to the SATNUM-values." << endl; // "minPoro" intentionally left undocumented // "saturationThreshold" also } static void usageandexit() { usage(); exit(1); } //! \brief Parse command line arguments into string map. //! \param[in,out] options The map of options. Should be filled with default values on entry. //! \param[in] varnum Number of arguments //! \param[in] vararg The arguments //! \param[in] verbose Whether or not to print parsed command line arguments //! \returns index of input file if positive, negated index to offending argument on failure. static int parseCommandLine(std::map& options, int varnum, char** vararg, bool verbose) { int argeclindex = 0; for (int argidx = 1; argidx < varnum; argidx += 2) { if (string(vararg[argidx]).substr(0,1) == "-") { string searchfor = string(vararg[argidx]).substr(1); // Chop off leading '-' /* Check if it is a match */ if (options.count(searchfor) == 1) { options[searchfor] = string(vararg[argidx+1]); if (verbose) cout << "Parsed command line option: " << searchfor << " := " << vararg[argidx+1] << endl; argeclindex = argidx + 2; } else return -argidx; } else { // if vararg[argidx] does not start in '-', // assume we have found the position of the Eclipse-file. argeclindex = argidx; break; // out of for-loop, } } return argeclindex; } //! \brief Return eclipse-style output filename. //! \param[in] opfname Base output file name. //! \param[in] comp Component (X, Y, Z). //! \param[in] sat Fluid system type. //! \return Eclipse-style filename for requested component/fluid system combination. static std::string getEclipseOutputFile(const std::string& opfname, char comp, char sat) { string fnbase = opfname.substr(0,opfname.find_last_of('.')); return fnbase + "-" +comp + ".S" + sat + "OF"; } //! \brief Write eclipse-style output file. //! \param[in] RelPermValues RelPerm values to write. //! \param[in] Satvalues Saturation values to write. //! \param[in] Pvalues Pressure values to write. //! \param[in] options Option structure. //! \param[in] component Component to write (0..2). //! \param[in] owsystem Fluid system type. template static void writeEclipseOutput(Lazy& RelPermValues, const std::vector& Satvalues, const std::vector& Pvalues, std::map& options, int component, bool owsystem) { std::stringstream swof; char sat = (owsystem?'W':'G'); char comp = 'x'+component; std::string krowstring = std::string("krow") + comp + "swirr"; double krowswirr = atof(options[krowstring].c_str()); const int outputprecision = atoi(options["outputprecision"].c_str()); const int fieldwidth = outputprecision + 8; // x-direction swof << "-- This file is based on the results in " << endl << "-- " << options["output"] << endl << "-- for relperm in " << comp << "-direction." << endl << "-- Pressure values (Pc) given in bars." << endl << "-- S" << (char)std::tolower(sat) << " Kr" << (char)std::tolower(sat) << comp << comp << " Kro" << (char)std::tolower(sat) << comp << comp << " Pc(bar)" << endl << "--S" << sat << "OF" << endl; if (krowswirr > 0) { swof << showpoint << setw(fieldwidth) << setprecision(outputprecision) << 0 << showpoint << setw(fieldwidth) << setprecision(outputprecision) << 0 << showpoint << setw(fieldwidth) << setprecision(outputprecision) << krowswirr << showpoint << setw(fieldwidth) << setprecision(outputprecision) << 0 << endl; } for (size_t i=0; i < Satvalues.size(); ++i) { swof << showpoint << setw(fieldwidth) << setprecision(outputprecision) << Satvalues[i] << showpoint << setw(fieldwidth) << setprecision(outputprecision) << RelPermValues[0][component][i] << showpoint << setw(fieldwidth) << setprecision(outputprecision) << RelPermValues[1][component][i] << showpoint << setw(fieldwidth) << setprecision(outputprecision) << Pvalues[i]/100000.0 << endl; } swof << "/" << endl; std::ofstream file; file.open(getEclipseOutputFile(options["output"], std::toupper(comp), sat), std::ios::out | std::ios::trunc); file << swof.str(); file.close(); } namespace { std::map defineOptions() { return { {"bc", "f"}, // Fixed boundary conditions {"points", "30"}, // Number of saturation points (uniformly distributed within saturation endpoints) {"relPermCurve", "2"}, // Which column in the rock types are upscaled {"upscaleBothPhases", "true"}, // Whether to upscale for both phases in the same run. Default true. {"jFunctionCurve", "4"}, // Which column in the rock type file is the J-function curve {"surfaceTension", "11"}, // Surface tension given in dynes/cm {"output", ""}, // If this is set, output goes to screen and to this file. {"gravity", "0.0"}, // default is no gravitational effects {"waterDensity", "1.0"}, // default density of water, only applicable to gravity {"oilDensity", "0.6"}, // ditto {"interpolate", "0"}, // default is not to interpolate {"maxpoints", "1000"}, // maximal number of saturation points. {"outputprecision", "4"}, // number of significant numbers to print {"maxPermContrast", "1e7"}, // maximum allowed contrast in each single-phase computation {"minPerm", "1e-12"}, // absolute minimum for allowed cell permeability {"maxPerm", "100000"}, // maximal allowed cell permeability {"minPoro", "0.0001"}, // this limit is necessary for pcmin/max computation {"saturationThreshold", "0.00001"}, // accuracy threshold for saturation, we ignore Pc values that // give so small contributions near endpoints. {"linsolver_tolerance", "1e-12"}, // residual tolerance for linear solver {"linsolver_verbosity", "0"}, // verbosity level for linear solver {"linsolver_max_iterations", "0"}, // Maximum number of iterations allow, specify 0 for default {"linsolver_type", "3"}, // Type of linear solver: 0 = ILU0/CG, 1 = AMG/CG, 2 KAMG/CG, 3 FAST_AMG/CG {"linsolver_prolongate_factor", "1.0"}, // Prolongation factor in AMG {"linsolver_smooth_steps", "1"}, // Number of smoothing steps in AMG {"fluids", "ow"}, // Whether upscaling for oil/water (ow) or gas/oil (go) {"krowxswirr", "-1"}, // Relative permeability in x direction of oil in corresponding oil/water system {"krowyswirr", "-1"}, // Relative permeability in y direction of oil in corresponding oil/water system {"krowzswirr", "-1"}, // Relative permeability in z direction of oil in corresponding oil/water system {"doEclipseCheck", "true"}, // Check if minimum relpermvalues in input are zero (specify critical saturations) {"critRelpermThresh", "1e-6"}, // Threshold for setting minimum relperm to 0 (thus specify critical saturations) }; } void invertCapillaryPressureIsotropic(const char* ROCKFILENAME, const int i, const int jFunctionCurve, const int relPermCurve, std::vector& JfunctionNames, Opm::RelPermUpscaleHelper& helper) { MonotCubicInterpolator Jtmp; try { Jtmp = MonotCubicInterpolator(ROCKFILENAME, 1, jFunctionCurve); } catch (const char* errormessage) { std::stringstream str; str << "Error: " << errormessage << endl << "Check filename and -jFunctionCurve" << endl; throw std::runtime_error(str.str()); } // Invert J-function, now we get saturation as a function of pressure: if (Jtmp.isStrictlyMonotone()) { helper.InvJfunctions.emplace_back(Jtmp.get_fVector(), Jtmp.get_xVector()); JfunctionNames.push_back(ROCKFILENAME); auto& krfunc = helper.Krfunctions[0]; if (helper.upscaleBothPhases) { krfunc[0].emplace_back(ROCKFILENAME, 1, 2); krfunc[1].emplace_back(ROCKFILENAME, 1, 3); } else { krfunc[0].emplace_back(ROCKFILENAME, 1, relPermCurve); } } else { std::stringstream str; str << "Error: Jfunction " << (i + 1) << " in rock file " << ROCKFILENAME << " was not invertible."; throw std::runtime_error(str.str()); } } void invertCapillaryPressureAnIsotropic(const char* ROCKFILENAME, const int i, const int /* jFunctionCurve */, const int /* relPermCurve */, std::vector& JfunctionNames, Opm::RelPermUpscaleHelper& helper) { MonotCubicInterpolator Pctmp; try { Pctmp = MonotCubicInterpolator(ROCKFILENAME, 2, 1); } catch (const char* errormessage) { std::stringstream str; str << "Error: " << errormessage << '\n' << "Check filename and columns 1 and 2 (Pc and " << helper.saturationstring << ')'; throw std::domain_error(str.str()); } // Invert Pc(Sw) curve into Sw(Pc): if (Pctmp.isStrictlyMonotone()) { helper.SwPcfunctions.emplace_back(Pctmp.get_fVector(), Pctmp.get_xVector()); JfunctionNames.push_back(ROCKFILENAME); auto& krfunc = helper.Krfunctions; krfunc[0][0].emplace_back(ROCKFILENAME, 2, 3); krfunc[1][0].emplace_back(ROCKFILENAME, 2, 4); krfunc[2][0].emplace_back(ROCKFILENAME, 2, 5); if (helper.upscaleBothPhases) { krfunc[0][1].emplace_back(ROCKFILENAME, 2, 6); krfunc[1][1].emplace_back(ROCKFILENAME, 2, 7); krfunc[2][1].emplace_back(ROCKFILENAME, 2, 8); } } else { std::stringstream str; str << "Error: Pc(" << helper.saturationstring << ") curve " << (i + 1) << " in rock file " << ROCKFILENAME << " was not invertible."; throw std::runtime_error(str.str()); } } void extractSaturationFunctions(const int varnum, char* vararg[], const int rockfileindex, const int stone_types, const int jFunctionCurve, const int relPermCurve, std::vector& JfunctionNames, Opm::RelPermUpscaleHelper& helper) { // If input is anisotropic, then we are in second mode with // different input file format auto invertPc = helper.anisotropic_input ? invertCapillaryPressureAnIsotropic : invertCapillaryPressureIsotropic; for (int i = 0 ; i < stone_types; ++i) { const char* ROCKFILENAME = vararg[(rockfileindex + stone_types == varnum) ? rockfileindex + i : rockfileindex]; // Check if rock file exists and is readable: { std::ifstream rockfile(ROCKFILENAME, std::ios::in); if (!rockfile) { std::stringstream str; str << "Error: Filename " << ROCKFILENAME << " not found or not readable."; throw std::runtime_error(str.str()); } } invertPc(ROCKFILENAME, i, jFunctionCurve, relPermCurve, JfunctionNames, helper); } } } int main(int varnum, char** vararg) try { // Variables used for timing/profiling: clock_t start, finish; double timeused = 0.0, timeused_tesselation = 0.0; double timeused_upscale_wallclock = 0.0; /****************************************************************************** * Step 1: * Process command line options */ Dune::MPIHelper& mpi=Dune::MPIHelper::instance(varnum, vararg); const int mpi_rank = mpi.rank(); #ifdef HAVE_MPI const int mpi_nodecount = mpi.size(); #endif if (varnum == 1) { /* If no arguments supplied ("upscale_relperm" is the first "argument") */ usage(); exit(1); } /* Populate options-map with default values */ auto options = defineOptions(); /* Check first if there is anything on the command line to look for */ if (varnum == 1) { if (mpi_rank == 0) cout << "Error: No Eclipsefile or stonefiles found on command line." << endl; usageandexit(); } /* 'argeclindex' is so that vararg[argeclindex] = the eclipse filename. */ int argeclindex = parseCommandLine(options, varnum, vararg, mpi_rank == 0); if (argeclindex < 0) { if (mpi_rank == 0) cout << "Option -" << vararg[-argeclindex] << " unrecognized." << endl; usageandexit(); } RelPermUpscaleHelper helper(mpi_rank, options); const auto owsystem = helper.saturationstring == "Sw"; // argeclindex should now point to the eclipse file static char* ECLIPSEFILENAME(vararg[argeclindex]); argeclindex += 1; // argeclindex jumps to next input argument, now it points to the stone files. // argeclindex now points to the first J-function. This index is not // to be touched now. static int rockfileindex = argeclindex; /* Check if at least one J-function is supplied on command line */ if (varnum <= rockfileindex) throw std::runtime_error("Error: No J-functions found on command line."); /* Check validity of boundary conditions chosen, and make booleans for boundary conditions, this allows more readable code later. */ helper.setupBoundaryConditions(); bool isFixed = helper.boundaryCondition == SinglePhaseUpscaler::Fixed, isLinear = helper.boundaryCondition == SinglePhaseUpscaler::Linear, isPeriodic = helper.boundaryCondition == SinglePhaseUpscaler::Periodic; // If this number is 1 or higher, the output will be interpolated, if not // the computed data is untouched. const int interpolationPoints = atoi(options["interpolate"].c_str()); const auto doInterpolate = interpolationPoints > 1; /*********************************************************************** * Step 2: * Load geometry and data from Eclipse file */ // Read data from the Eclipse file and // populate our vectors with data from the file // Test if filename exists and is readable start = clock(); auto deck = RelPermUpscaleHelper::parseEclipseFile(ECLIPSEFILENAME); finish = clock(); timeused = (double(finish) - double(start)) / CLOCKS_PER_SEC; if (helper.isMaster) { cout << " (" << timeused << " secs)" << endl; } { const double minPerm = atof(options["minPerm"].c_str()); const double maxPerm = atof(options["maxPerm"].c_str()); const double minPoro = atof(options["minPoro"].c_str()); helper.sanityCheckInput(deck, minPerm, maxPerm, minPoro); } /*************************************************************************** * Step 3: * Load relperm- and J-function-curves for the stone types. * We read columns from text-files, syntax allowed is determined * by MonotCubicInterpolator which actually opens and parses the * text files. * * If a standard eclipse data file is given as input, the data columns * should be: * Sw Krw Kro J-func * (In this case, the option -relPermCurve determines which of Krw or Kro is used) * * If output from this very program is given as input, then the data columns read * Pc Sw Krx Kry Krz * * (and the option -relPermCurve and -jFunctionCurve are ignored) * * How do we determine which mode of operation? * - If PERMY and PERMZ are present in grdecl-file, we are in the anisotropic mode * */ // Number of stone-types is max(satnums): // If there is only one J-function supplied on the command line, // use that for all stone types. const auto stone_types = static_cast(*(max_element(helper.satnums.begin(), helper.satnums.end()))); // This decides whether we are upscaling both phases in this run or only one helper.upscaleBothPhases = (options["upscaleBothPhases"] == "true"); helper.points = atoi(options["points"].c_str()); // Input for surfaceTension is dynes/cm // SI units are Joules/square metre => multiply with 10^-3 to obtain SI units. const double surfaceTension = atof(options["surfaceTension"].c_str()) * 1e-3; const double gravity = atof(options["gravity"].c_str()); const bool includeGravity = (std::fabs(gravity) > DBL_MIN); // true for non-zero gravity const int outputprecision = atoi(options["outputprecision"].c_str()); // Handle two command line input formats, either one J-function for all stone types // or one each. If there is only one stone type, both code blocks below are equivalent. if (! ((varnum == rockfileindex + stone_types) || (varnum == rockfileindex + 1 ))) { throw std::runtime_error("Error: Wrong number of stone-functions provided."); } // Placeholder for the names of the loaded J-functions. std::vector JfunctionNames; { // This decides whether we are upscaling water or oil relative permeability const int relPermCurve = atoi(options["relPermCurve"].c_str()); const int jFunctionCurve = atoi(options["jFunctionCurve"].c_str()); extractSaturationFunctions(varnum, vararg, rockfileindex, stone_types, jFunctionCurve, relPermCurve, JfunctionNames, helper); } // Check if input relperm curves satisfy ECLIPSE requirement of // specifying critical saturations if (helper.doEclipseCheck) { helper.checkCriticalSaturations(); } /***************************************************************************** * Step 4: * Generate tesselated grid: * This is a step needed for the later discretization code to figure out which * cells are connected to which. Each cornerpoint-cell is tesselated into 8 tetrahedrons. * * In case of non-zero gravity, calculate z-values of every cell: * 1) Compute height of model by averaging z-values of the top layer corners. * 2) Calculate density difference between phases in SI-units * 3) Go through each cell and find the z-values of the eight corners of the cell. * Set height of cell equal to average of z-values of the corners minus half of * model height. Now the cell height is relative to model centre. * Set pressure difference for the cell equal to density difference times gravity * constant times cell height times factor 10^-7 to obtain bars (same as p_c) */ timeused_tesselation = helper.tesselateGrid(deck); /* If gravity is to be included, calculate z-values of every cell: */ if (includeGravity) { helper.calculateCellPressureGradients(); } /****************************************************************************** * Step 5: * Go through each cell and calculate the minimum and * maximum capillary pressure possible in the cell, given poro, * perm and the J-function for the cell. This depends on the * J-function in that they represent all possible saturations, * ie. we do not want to extrapolate the J-functions (but we might * have to do that later in the computations). * * The user-supplied surface tension is ignored until * the final output of results. */ helper.calculateMinMaxCapillaryPressure(); /*************************************************************************** * Step 6: * Upscale capillary pressure function. * * This is upscaled in advance in order to be able to have uniformly distributed * saturation points for which upscaling is performed. * * Capillary pressure points are chosen heuristically in order to * ensure largest saturation interval between two saturation points * is 1/500 of the saturation interval. Monotone cubic interpolation * will be used afterwards for accessing the tabulated values. */ helper.upscaleCapillaryPressure(); /***************************************************************************** * Step 7: * Upscale single phase permeability * This uses the PERMX in the eclipse file as data, and upscales using * fixed boundary (no-flow) conditions * * In an MPI-environment, this is only done on the master node. */ helper.upscaleSinglePhasePermeability(); /***************************************************************** * Step 8: * * Loop through a given number of uniformly distributed saturation points * and upscale relative permeability for each of them. * a: Make vector of capillary pressure points corresponding to uniformly * distributed water saturation points between saturation endpoints. * b: Loop over capillary pressure points * 1) Loop over all cells to find the saturation value given the * capillary pressure found in (a). Given the saturation value, find the * phase permeability in the cell given input relperm curve and input * permeability values. * 2) Upscale phase permeability for the geometry. * c: Calculate relperm tensors from all the phase perm tensors. */ double avg_upscaling_time_pr_point; std::tie(timeused_upscale_wallclock, avg_upscaling_time_pr_point) = helper.upscalePermeability(mpi_rank); /* * Step 8c: Make relperm values from phaseperms * (only master node can do this) */ std::array>,2> RelPermValues; RelPermValues[0] = helper.getRelPerm(0); if (helper.upscaleBothPhases) { RelPermValues[1] = helper.getRelPerm(1); } /********************************************************************************* * Step 9 * * Output results to stdout and optionally to file. Note, we only output to * file if the '-outputWater'-option and/or '-outputOil' has been set, as this option is an * empty string by default. */ if (helper.isMaster) { stringstream outputtmp; // Print a table of all computed values: outputtmp << "######################################################################\n" << "# Results from upscaling relative permeability.\n" << "#\n"; #if defined(HAVE_MPI) && HAVE_MPI outputtmp << "# (MPI-version)\n"; #endif // HAVE_MPI { time_t now = std::time(NULL); outputtmp << "# Finished: " << asctime(localtime(&now)); } { utsname hostname; uname(&hostname); outputtmp << "# Hostname: " << hostname.nodename << endl; } outputtmp << "#" << '\n' << "# ECLIPSE file: " << ECLIPSEFILENAME << '\n' << "# cells: " << helper.tesselatedCells << '\n' << "# Pore volume: " << helper.poreVolume << '\n' << "# volume: " << helper.volume << '\n' << "# Porosity: " << (helper.poreVolume / helper.volume) << '\n' << "#\n"; if (! helper.anisotropic_input) { const auto& invJ = helper.InvJfunctions; for (int i = 0; i < stone_types ; ++i) { outputtmp << "# Stone " << (i + 1) << ": " << JfunctionNames[i] << " (" << invJ[i].getSize() << " points)\n"; } outputtmp << "# jFunctionCurve: " << options["jFunctionCurve"] << '\n'; if (! helper.upscaleBothPhases) { outputtmp << "# relPermCurve: " << options["relPermCurve"] << '\n'; } } else { // anisotropic input, not J-functions that are supplied on // command line (but vector JfunctionNames is still used) const auto& kr = helper.Krfunctions[0][0]; for (int i = 0; i < stone_types ; ++i) { outputtmp << "# Stone " << (i + 1) << ": " << JfunctionNames[i] << " (" << kr[i].getSize() << " points)\n"; } } outputtmp << "#" << '\n' << "# Timings: Tesselation: " << timeused_tesselation << " secs" << '\n' << "# Upscaling: " << timeused_upscale_wallclock << " secs"; #ifdef HAVE_MPI outputtmp << " (wallclock time)\n" << "# " << avg_upscaling_time_pr_point << " secs pr. saturation point" << '\n' << "# MPI-nodes: " << mpi_nodecount << '\n'; // Single phase upscaling time is included here, in possibly a hairy // way. const auto speedup = (avg_upscaling_time_pr_point * (helper.points + 1) + timeused_tesselation) / (timeused_upscale_wallclock + avg_upscaling_time_pr_point + timeused_tesselation); outputtmp << "# Speedup: " << speedup << ", efficiency: " << (speedup / mpi_nodecount) << endl; #else // !HAVE_MPI outputtmp << ", " << avg_upscaling_time_pr_point << " secs avg for " << helper.points << " runs" << endl; #endif // HAVE_MPI outputtmp << "#\n" << "# Options used:\n" << "# Boundary conditions: "; if (isFixed) { outputtmp << "Fixed (no-flow)\n"; } if (isPeriodic) { outputtmp << "Periodic\n"; } if (isLinear) { outputtmp << "Linear\n"; } outputtmp << "# points: " << options["points"] << '\n' << "# maxPermContrast: " << options["maxPermContrast"] << '\n' << "# minPerm: " << options["minPerm"] << " [mD]\n" << "# minPoro: " << options["minPoro"] << '\n' << "# surfaceTension: " << options["surfaceTension"] << " dynes/cm\n"; if (includeGravity) { outputtmp << "# gravity: " << options["gravity"] << " m/s²\n"; if (owsystem) { outputtmp << "# waterDensity: " << options["waterDensity"] << " g/cm³\n"; } else { outputtmp << "# gasDensity: " << options["waterDensity"] << " g/cm³\n"; } outputtmp << "# oilDensity: " << options["oilDensity"] << " g/cm³\n"; } else { outputtmp << "# gravity: 0\n"; } if (doInterpolate) { outputtmp << "# interpolate: " << options["interpolate"] << " points\n"; } { auto mD = [](const double k) { return Opm::unit::convert::to(k, Opm::prefix::milli*Opm::unit::darcy); }; outputtmp << "#\n" << "# Single phase permeability (mD)\n" << "# |Kxx Kxy Kxz| = " << mD(helper.permTensor(0,0)) << " " << mD(helper.permTensor(0,1)) << " " << mD(helper.permTensor(0,2)) << '\n'; outputtmp << "# |Kyx Kyy Kyz| = " << mD(helper.permTensor(1,0)) << " " << mD(helper.permTensor(1,1)) << " " << mD(helper.permTensor(1,2)) << '\n'; outputtmp << "# |Kzx Kzy Kzz| = " << mD(helper.permTensor(2,0)) << " " << mD(helper.permTensor(2,1)) << " " << mD(helper.permTensor(2,2)) << '\n'; outputtmp << "#\n"; } if (doInterpolate) { outputtmp << "# NB: Data points shown are interpolated.\n"; } outputtmp << "######################################################################\n"; if (helper.upscaleBothPhases) { string phase1, phase2; if (owsystem) { phase1 = "w"; } else { phase1 = "g"; } phase2 = "o"; if (isFixed) { outputtmp << "# Pc (Pa) " << helper.saturationstring << " Kr" << phase1 << "xx Kr" << phase1 << "yy Kr" << phase1 << "zz Kr" << phase2 << "xx Kr" << phase2 << "yy Kr" << phase2 << "zz\n"; } else if (isPeriodic || isLinear) { outputtmp << "# Pc (Pa) " << helper.saturationstring << " Kr" << phase1 << "xx Kr" << phase1 << "yy Kr" << phase1 << "zz Kr" << phase1 << "yz Kr" << phase1 << "xz Kr" << phase1 << "xy Kr" << phase1 << "zy Kr" << phase1 << "zx Kr" << phase1 << "yx Kr" << phase2 << "xx Kr" << phase2 << "yy Kr" << phase2 << "zz Kr" << phase2 << "yz Kr" << phase2 << "xz Kr" << phase2 << "xy Kr" << phase2 << "zy Kr" << phase2 << "zx Kr" << phase2 << "yx\n"; } } else { if (isFixed) { outputtmp << "# Pc (Pa) " << helper.saturationstring << " Krxx Kryy Krzz\n"; } else if (isPeriodic || isLinear) { outputtmp << "# Pc (Pa) " << helper.saturationstring << " Krxx Kryy Krzz " << "Kryz Krxz Krxy Krzy " << "Krzx Kryx\n"; } } auto Pvalues = helper.pressurePoints; // Multiply all pressures with the surface tension (potentially) // supplied at the command line. This multiplication has been // postponed to here to avoid division by zero and to avoid special // handling of negative capillary pressure in the code above. for (auto& press : Pvalues) { press *= surfaceTension; } auto Satvalues = helper.WaterSaturation; //.get_fVector(); // If user wants interpolated output, do monotone cubic interpolation // by modifying the data vectors that are to be printed. if (doInterpolate) { // Split interval between minumum and maximum saturation into // (interpolationPoints - 1) equally sized smaller intervals. auto SatvaluesInterp = std::vector{}; { const auto m = std::minmax_element(Satvalues.begin(), Satvalues.end()); const auto xmin = *m.first; const auto xmax = *m.second; // Make uniform grid in saturation axis SatvaluesInterp.reserve(interpolationPoints); auto interp = [xmin, xmax, interpolationPoints](const int i) -> double { const auto t = static_cast(i) / (interpolationPoints - 1); return xmin + t*(xmax - xmin); }; for (int i = 0; i < interpolationPoints; ++i) { SatvaluesInterp.push_back(interp(i)); } } // Now capillary pressure and computed relperm-values must be // viewed as functions of saturation, and then interpolated on // the uniform saturation grid. // Now overwrite existing Pvalues and relperm-data with // interpolated data: { const auto PvaluesVsSaturation = MonotCubicInterpolator(Satvalues, Pvalues); Pvalues.clear(); Pvalues.reserve(SatvaluesInterp.size()); for (const auto& s : SatvaluesInterp) { Pvalues.push_back(PvaluesVsSaturation.evaluate(s)); } } for (int voigtIdx = 0; voigtIdx < helper.tensorElementCount; ++voigtIdx) { auto& kr = RelPermValues[0][voigtIdx]; const auto RelPermVsSaturation = MonotCubicInterpolator(Satvalues, kr); kr.clear(); kr.reserve(SatvaluesInterp.size()); for (const auto& s : SatvaluesInterp) { kr.push_back(RelPermVsSaturation.evaluate(s)); } } if (helper.upscaleBothPhases) { for (int voigtIdx = 0; voigtIdx < helper.tensorElementCount; ++voigtIdx) { auto& kr = RelPermValues[1][voigtIdx]; const auto RelPermVsSaturation = MonotCubicInterpolator(Satvalues, kr); kr.clear(); kr.reserve(SatvaluesInterp.size()); for (const auto& s : SatvaluesInterp) { kr.push_back(RelPermVsSaturation.evaluate(s)); } } } // Now also overwrite Satvalues Satvalues = std::move(SatvaluesInterp); } // The code below does not care whether the data is interpolated or not. const auto fieldwidth = outputprecision + 8; for (decltype(Satvalues.size()) i = 0, n = Satvalues.size(); i < n; ++i) { outputtmp << showpoint << setw(fieldwidth) << setprecision(outputprecision) << Pvalues[i]; outputtmp << showpoint << setw(fieldwidth) << setprecision(outputprecision) << Satvalues[i]; for (int voigtIdx = 0; voigtIdx < helper.tensorElementCount; ++voigtIdx) { outputtmp << showpoint << setw(fieldwidth) << setprecision(outputprecision) << RelPermValues[0][voigtIdx][i]; } if (helper.upscaleBothPhases) { for (int voigtIdx = 0; voigtIdx < helper.tensorElementCount; ++voigtIdx) { outputtmp << showpoint << setw(fieldwidth) << setprecision(outputprecision) << RelPermValues[1][voigtIdx][i]; } } outputtmp << endl; } cout << outputtmp.str(); if (options["output"] != "") { cout << "Writing results to " << options["output"] << endl; ofstream outfile; outfile.open(options["output"].c_str(), ios::out | ios::trunc); outfile << outputtmp.str(); outfile.close(); } // If both phases are upscaled and output is specified, create SWOF // or SGOF files for ECLIPSE if (options["output"] != "" && helper.upscaleBothPhases) { cout << "Writing ECLIPSE compatible files to " << getEclipseOutputFile(options["output"], 'X', owsystem ? 'W' : 'G') << ", " << getEclipseOutputFile(options["output"], 'Y', owsystem ? 'W' : 'G') << " and " << getEclipseOutputFile(options["output"], 'Z', owsystem ? 'W' : 'G') << endl; for (int comp = 0; comp < 3; ++comp) { writeEclipseOutput(RelPermValues, Satvalues, Pvalues, options, comp, owsystem); } } } } catch (const std::exception& e) { std::cerr << e.what() << '\n'; usageandexit(); }