diff --git a/CMakeLists_files.cmake b/CMakeLists_files.cmake index 04d6fb6..292e65d 100644 --- a/CMakeLists_files.cmake +++ b/CMakeLists_files.cmake @@ -64,10 +64,10 @@ list (APPEND TEST_DATA_FILES tests/input_data/grids/stone2.txt tests/input_data/grids/stoneAniso.txt tests/input_data/grids/stonefile_benchmark.txt - tests/input_data/reference_solutions/upscale_perm_BCp_PeriodicTilted.txt - tests/input_data/reference_solutions/upscale_perm_BCflp_27cellsAniso.txt - tests/input_data/reference_solutions/upscale_perm_BCflp_27cellsIso.txt - tests/input_data/reference_solutions/upscale_perm_BCfl_EightCells.txt + tests/input_data/reference_solutions/upscale_perm_BCp_PeriodicTilted.txt + tests/input_data/reference_solutions/upscale_perm_BCflp_27cellsAniso.txt + tests/input_data/reference_solutions/upscale_perm_BCflp_27cellsIso.txt + tests/input_data/reference_solutions/upscale_perm_BCfl_EightCells.txt tests/input_data/reference_solutions/upscale_perm_BCflp_Hummocky.txt tests/input_data/reference_solutions/upscale_relperm_BCf_pts20_surfTens11_stonefile_benchmark_stonefile_benchmark_benchmark_tiny_grid.txt tests/input_data/reference_solutions/upscale_relperm_BCf_pts30_surfTens11_stone1_stone1_EightCells.txt diff --git a/benchmarks/upscale_relperm_benchmark.cpp b/benchmarks/upscale_relperm_benchmark.cpp index 8761539..4f3129d 100644 --- a/benchmarks/upscale_relperm_benchmark.cpp +++ b/benchmarks/upscale_relperm_benchmark.cpp @@ -27,7 +27,7 @@ The main machinery is the same, but with some changes: - Input data (grid file, rock file and reference solution) is not provided from command line, - but is built in at compiler time by embedding hexadecimal (1 byte) input data files. See + but is built in at compiler time by embedding hexadecimal (1 byte) input data files. See README for further documentation. - Other command line options are not supported. - The construction of deck and stone data is changed due to change in input routine. diff --git a/examples/aniso_implicitcap_test.cpp b/examples/aniso_implicitcap_test.cpp index 36a0146..5c4db47 100644 --- a/examples/aniso_implicitcap_test.cpp +++ b/examples/aniso_implicitcap_test.cpp @@ -81,5 +81,3 @@ catch (const std::exception &e) { std::cerr << "Program threw an exception: " << e.what() << "\n"; throw; } - - diff --git a/examples/cpchop.cpp b/examples/cpchop.cpp index d169f8f..2169d0c 100644 --- a/examples/cpchop.cpp +++ b/examples/cpchop.cpp @@ -109,7 +109,7 @@ try Opm::unit::square(Opm::unit::meter)); // Input for surfaceTension is dynes/cm, SI units are Joules/square metre - const double surfaceTension = param.getDefault("surfaceTension", 11.0) * 1e-3; // multiply with 10^-3 to obtain SI units + const double surfaceTension = param.getDefault("surfaceTension", 11.0) * 1e-3; // multiply with 10^-3 to obtain SI units bool endpoints = param.getDefault("endpoints", false); // whether to upscale saturation endpoints bool cappres = param.getDefault("cappres", false); // whether to upscale capillary pressure @@ -158,18 +158,18 @@ try if (!rock_stream) { OPM_THROW(std::runtime_error, "Could not open file " + rockfilename); } - + if (! anisorocks) { //Isotropic input rocks (Sw Krw Kro J) Opm::MonotCubicInterpolator Jtmp; try { - Jtmp = Opm::MonotCubicInterpolator(rockname, 1, 4); + Jtmp = Opm::MonotCubicInterpolator(rockname, 1, 4); } catch (const char * errormessage) { std::cerr << "Error: " << errormessage << std::endl; std::cerr << "Check filename" << std::endl; exit(1); } - + // Invert J-function, now we get saturation as a function of pressure: if (Jtmp.isStrictlyMonotone()) { InvJfunctions.push_back(Opm::MonotCubicInterpolator(Jtmp.get_fVector(), Jtmp.get_xVector())); @@ -178,7 +178,7 @@ try std::cerr << "Error: Jfunction " << i+1 << " in rock file " << rockname << " was not invertible." << std::endl; exit(1); } - + jfuncendpoints_[i][0] = Jtmp.getMinimumX().second; jfuncendpoints_[i][1] = Jtmp.getMaximumX().second; rocksatendpoints_[i][0] = Jtmp.getMinimumX().first; @@ -213,7 +213,7 @@ try } rocksatendpoints_[i][0] = Pctmp.getMinimumX().first; rocksatendpoints_[i][1] = Pctmp.getMaximumX().first; - } + } } } @@ -270,7 +270,7 @@ try std::cout << "***** WARNING: Unused parameters: *****\n"; param.displayUsage(); } - + // Note that end is included in interval for uniform_int. std::uniform_int_distribution<> disti(imin, imax - ilen); std::uniform_int_distribution<> distj(jmin, jmax - jlen); @@ -295,9 +295,9 @@ try std::vector > pcvalues; std::vector dipangs, azimuths; - // Initialize a matrix for subsample satnum volumes. + // Initialize a matrix for subsample satnum volumes. // Outer index is subsample index, inner index is SATNUM-value - std::vector > rockvolumes; + std::vector > rockvolumes; int maxSatnum = 0; // This value is determined from the chopped cells. int finished_subsamples = 0; // keep explicit count of successful subsamples @@ -323,7 +323,7 @@ try if (upscale) { auto subdeck = ch.subDeck(); Opm::SinglePhaseUpscaler upscaler; - + upscaler.init(subdeck, bctype, minpermSI, residual_tolerance, linsolver_verbosity, linsolver_type, false); @@ -334,7 +334,7 @@ try porosities.push_back(upscaler.upscalePorosity()); if (ch.hasNTG()) { netporosities.push_back(upscaler.upscaleNetPorosity()); - ntgs.push_back(upscaler.upscaleNTG()); + ntgs.push_back(upscaler.upscaleNTG()); } if (ch.hasSWCR()) { swcrs.push_back(upscaler.upscaleSWCR(ch.hasNTG())); @@ -355,7 +355,7 @@ try // Create single-phase upscaling object to get poro and perm values from the grid auto subdeck = ch.subDeck(); std::vector perms = subdeck["PERMX"].back().getRawDoubleData(); - Opm::SinglePhaseUpscaler upscaler; + Opm::SinglePhaseUpscaler upscaler; upscaler.init(subdeck, bctype, minpermSI, residual_tolerance, linsolver_verbosity, linsolver_type, false); std::vector satnums = subdeck["SATNUM"].back().getIntData(); @@ -407,7 +407,7 @@ try } // Total porevolume and total volume -> upscaled porosity: - double poreVolume = std::accumulate(cellPoreVolumes.begin(), + double poreVolume = std::accumulate(cellPoreVolumes.begin(), cellPoreVolumes.end(), 0.0); double Swir = Swirvolume/poreVolume; @@ -421,7 +421,7 @@ try double Ptestvalue; while (largestSaturationInterval > (Swor-Swir)/double(nsatpoints)) { if (Pcmax == Pcmin) { - // This is a dummy situation, we go through once and then + // This is a dummy situation, we go through once and then // we are finished (this will be triggered by zero permeability) Ptestvalue = Pcmin; largestSaturationInterval = 0; @@ -457,12 +457,12 @@ try double waterSaturationCell = 0.0; if (satnums[cell_idx] > 0) { // handle "no rock" cells with satnum zero double PtestvalueCell; - + PtestvalueCell = Ptestvalue; - - if (!anisorocks) { + + if (!anisorocks) { double Jvalue = sqrt(perms[cell_idx] * milliDarcyToSqMetre /poros[cell_idx]) * PtestvalueCell / surfaceTension; - waterSaturationCell + waterSaturationCell = InvJfunctions[int(satnums[cell_idx])-1].evaluate(Jvalue); } else { // anisotropic_input, then we do not do J-function-scaling @@ -476,7 +476,7 @@ try WaterSaturationVsCapPressure.chopFlatEndpoints(saturationThreshold); std::vector wattest = WaterSaturationVsCapPressure.get_fVector(); std::vector cprtest = WaterSaturationVsCapPressure.get_xVector(); - Opm::MonotCubicInterpolator CapPressureVsWaterSaturation(WaterSaturationVsCapPressure.get_fVector(), + Opm::MonotCubicInterpolator CapPressureVsWaterSaturation(WaterSaturationVsCapPressure.get_fVector(), WaterSaturationVsCapPressure.get_xVector()); std::vector pcs; for (int satp=0; satp mincellvolume) { - std::pair xydip = gridinspector.cellDips(i, j, k); + std::pair xydip = gridinspector.cellDips(i, j, k); xdips_subsample.push_back(xydip.first); ydips_subsample.push_back(xydip.second); } @@ -515,7 +515,7 @@ try // double azimuth = atan(xydip.first/xydip.second); // double dip = acos(1.0/sqrt(pow(xydip.first,2.0)+pow(xydip.second,2.0)+1.0)); // dips_subsample.push_back( xydip.first ); - // azims_subsample.push_back(atan(xydip.first/xydip.second)); + // azims_subsample.push_back(atan(xydip.first/xydip.second)); // Average xdips and ydips double xdipaverage = accumulate(xdips_subsample.begin(), xdips_subsample.end(), 0.0)/xdips_subsample.size(); @@ -525,7 +525,7 @@ try double azimuth = atan(xdipaverage/ydipaverage)+azimuthdisplacement; double dip = acos(1.0/sqrt(pow(xdipaverage,2.0)+pow(ydipaverage,2.0)+1.0)); dipangs.push_back(dip); - azimuths.push_back(azimuth); + azimuths.push_back(azimuth); } if (satnumvolumes) { @@ -565,7 +565,7 @@ try std::cerr << "Warning: Upscaling chopped subsample nr. " << sample << " failed, proceeding to next subsample\n"; } } - + // Make stream of output data, to be outputted to screen and optionally to file std::stringstream outputtmp; @@ -589,8 +589,8 @@ try outputtmp << "# (auto) seed: " << autoseed << std::endl; } else { - outputtmp << "# (manual) seed: " << userseed << std::endl; - } + outputtmp << "# (manual) seed: " << userseed << std::endl; + } outputtmp << "################################################################################################" << std::endl; outputtmp << "# id"; if (upscale) { @@ -628,7 +628,7 @@ try } outputtmp << " Swir Swor"; if (cappres) { - outputtmp << " Pc(Swir) Pc2 Pc3 Pc4 Pc(Swor)"; + outputtmp << " Pc(Swir) Pc2 Pc3 Pc4 Pc(Swor)"; } } if (dips) { @@ -659,7 +659,7 @@ try outputtmp << std::showpoint << std::setw(fieldwidth) << std::setprecision(outputprecision) << permyzs[sample-1] << '\t' << std::showpoint << std::setw(fieldwidth) << std::setprecision(outputprecision) << permxzs[sample-1] << '\t' << - std::showpoint << std::setw(fieldwidth) << std::setprecision(outputprecision) << permxys[sample-1] << '\t'; + std::showpoint << std::setw(fieldwidth) << std::setprecision(outputprecision) << permxys[sample-1] << '\t'; } if (ch.hasNTG()) { outputtmp << diff --git a/examples/cpchop_depthtrend.cpp b/examples/cpchop_depthtrend.cpp index 6cadca2..b8248bc 100644 --- a/examples/cpchop_depthtrend.cpp +++ b/examples/cpchop_depthtrend.cpp @@ -48,12 +48,12 @@ #include /** - This program is a variant of cpchop. Instead of subsampling randomly, + This program is a variant of cpchop. Instead of subsampling randomly, it picks subsamples downwards in a model. It is specifically designed for extracting laterally extensive subsamples in specified intervals - downwards through a model, in order to analyze depth trends in + downwards through a model, in order to analyze depth trends in porosity in particular - + Default is to pick a new subsample every meter (zresolution=1) The subsample height might be smaller or larger than zresolution. */ @@ -105,22 +105,22 @@ try double residual_tolerance = param.getDefault("residual_tolerance", 1e-8); double linsolver_verbosity = param.getDefault("linsolver_verbosity", 0); double linsolver_type = param.getDefault("linsolver_type", 1); - + // Check for unused parameters (potential typos). if (param.anyUnused()) { std::cout << "***** WARNING: Unused parameters: *****\n"; param.displayUsage(); } - - // Check that we do not have any user input + + // Check that we do not have any user input // that goes outside the coordinates described in // the cornerpoint file (runtime-exception will be thrown in case of error) - ch.verifyInscribedShoebox(imin, ilen, imax, + ch.verifyInscribedShoebox(imin, ilen, imax, jmin, jlen, jmax, zmin, zlen, zmax); std::mt19937 gen; - + // Seed the random number generators with the current time, unless specified on command line // Warning: Current code does not allow 0 for the seed!! if (userseed == 0) { @@ -129,14 +129,14 @@ try else { gen.seed(userseed); } - + // Note that end is included in interval for uniform_int. std::uniform_int_distribution<> disti(imin, imax - ilen); std::uniform_int_distribution<> distj(jmin, jmax - jlen); auto ri = [&disti, &gen] { return disti(gen); }; auto rj = [&distj, &gen] { return distj(gen); }; - + // Storage for results std::vector zstarts; std::vector porosities; @@ -144,7 +144,7 @@ try std::vector permys; std::vector permzs; - + /* z_start is the topmost point of the subsample to extract */ for (double zstart = 0.0; zstart <= zmax-zlen; zstart += zresolution) { /* Horizontally, we pick by random, even though default behaviour is @@ -169,35 +169,35 @@ try Opm::SinglePhaseUpscaler upscaler; upscaler.init(subdeck, Opm::SinglePhaseUpscaler::Fixed, minpermSI, residual_tolerance, linsolver_verbosity, linsolver_type, false); - + Opm::SinglePhaseUpscaler::permtensor_t upscaled_K = upscaler.upscaleSinglePhase(); upscaled_K *= (1.0/(Opm::prefix::milli*Opm::unit::darcy)); - - + + zstarts.push_back(zstart); porosities.push_back(upscaler.upscalePorosity()); permxs.push_back(upscaled_K(0,0)); permys.push_back(upscaled_K(1,1)); permzs.push_back(upscaled_K(2,2)); - - } + + } } catch (...) { std::cerr << "Warning: Upscaling chopped subsample at z=" << zstart << "failed, proceeding to next subsample\n"; } } - + if (upscale) { - + // Make stream of output data, to be outputted to screen and optionally to file std::stringstream outputtmp; - + outputtmp << "################################################################################################" << std::endl; outputtmp << "# Results from depth trend analysis on subsamples" << std::endl; outputtmp << "#" << std::endl; time_t now = time(NULL); outputtmp << "# Finished: " << asctime(localtime(&now)); - + utsname hostname; uname(&hostname); outputtmp << "# Hostname: " << hostname.nodename << std::endl; outputtmp << "#" << std::endl; @@ -209,10 +209,10 @@ try outputtmp << "# zresolution: " << zresolution << std::endl; outputtmp << "################################################################################################" << std::endl; outputtmp << "# zstart porosity permx permy permz" << std::endl; - + const int fieldwidth = outputprecision + 8; for (size_t sample = 1; sample <= porosities.size(); ++sample) { - outputtmp << + outputtmp << std::showpoint << std::setw(fieldwidth) << std::setprecision(outputprecision) << zstarts[sample-1] << '\t' << std::showpoint << std::setw(fieldwidth) << std::setprecision(outputprecision) << porosities[sample-1] << '\t' << std::showpoint << std::setw(fieldwidth) << std::setprecision(outputprecision) << permxs[sample-1] << '\t' << @@ -220,17 +220,17 @@ try std::showpoint << std::setw(fieldwidth) << std::setprecision(outputprecision) << permzs[sample-1] << '\t' << std::endl; } - + if (resultfile != "") { std::cout << "Writing results to " << resultfile << std::endl; std::ofstream outfile; outfile.open(resultfile.c_str(), std::ios::out | std::ios::trunc); outfile << outputtmp.str(); - outfile.close(); + outfile.close(); } - - - + + + std::cout << outputtmp.str(); } } @@ -238,4 +238,3 @@ catch (const std::exception &e) { std::cerr << "Program threw an exception: " << e.what() << "\n"; throw; } - diff --git a/examples/cpregularize.cpp b/examples/cpregularize.cpp index 2f9fdd3..6a79688 100644 --- a/examples/cpregularize.cpp +++ b/examples/cpregularize.cpp @@ -28,7 +28,7 @@ your asked-for horizontal resolution does not divide the initial number of pillars in x/y, you will not obtain a fully regular grid, but still easier numerically. - + - Be careful with non-flat top and bottom boundary. */ @@ -76,7 +76,7 @@ try } Dune::MPIHelper::instance(argc, argv); - + Opm::ParameterGroup param(argc, argv); std::string gridfilename = param.get("gridfilename"); Opm::CornerPointChopper ch(gridfilename); @@ -105,15 +105,15 @@ try param.displayUsage(); } - // Check that we do not have any user input + // Check that we do not have any user input // that goes outside the coordinates described in // the cornerpoint file (runtime-exception will be thrown in case of error) // (ilen, jlen and zlen set to zero, does not apply here) - ch.verifyInscribedShoebox(imin, 0, imax, + ch.verifyInscribedShoebox(imin, 0, imax, jmin, 0, jmax, zmin, 0, zmax); - + // Storage for properties for regularized cells std::vector poro; std::vector permx; @@ -159,7 +159,7 @@ try zcorn_c.push_back(zmin + zidx_c * (zmax-zmin)/zres); } zcorn_c.push_back(zmax); - + // Run through the new regular grid to find its properties @@ -183,13 +183,13 @@ try Opm::SinglePhaseUpscaler upscaler; upscaler.init(subdeck, Opm::SinglePhaseUpscaler::Fixed, minpermSI, residual_tolerance, linsolver_verbosity, linsolver_type, false); - + Opm::SinglePhaseUpscaler::permtensor_t upscaled_K = upscaler.upscaleSinglePhase(); upscaled_K *= (1.0/(Opm::prefix::milli*Opm::unit::darcy)); poro.push_back(upscaler.upscalePorosity()); permx.push_back(upscaled_K(0,0)); permy.push_back(upscaled_K(1,1)); - permz.push_back(upscaled_K(2,2)); + permz.push_back(upscaled_K(2,2)); } catch (...) { std::cout << "Warning: Upscaling for cell failed to convert, values set to zero\n"; @@ -207,7 +207,7 @@ try std::cerr << "Could not open file " << resultgrid << "\n"; throw std::runtime_error("Could not open output file."); } - out << "SPECGRID\n" << ires << ' ' << jres << ' ' << zres + out << "SPECGRID\n" << ires << ' ' << jres << ' ' << zres << " 1 F\n/\n\n"; out << "COORD\n"; @@ -223,12 +223,12 @@ try Write ZCORN, that is the Z-coordinates along the pillars, specifying the eight corners of each cell. Each corner is specified for each cell, even though it is the same corner that is used in other - cells. + cells. We loop over corners in each grid cell, directions: z, y, x (x innermost). The code here *IS* redundant, but the grid is also very redundant for a grid that is really regular.. - */ + */ out << "ZCORN\n"; double zlen = zmax-zmin; for (int zidx=0; zidx < zres; ++zidx) { @@ -252,36 +252,34 @@ try } } out << "/\n\n"; - + out << "PORO\n"; for (size_t idx=0; idx < (size_t)poro.size(); ++idx) { out << poro[idx] << std::endl; } out << "/\n\n"; - + out << "PERMX\n"; for (size_t idx=0; idx < (size_t)permx.size(); ++idx) { out << permx[idx] << std::endl; } out << "/\n\n"; - + out << "PERMY\n\n"; for (size_t idx=0; idx < (size_t)permy.size(); ++idx) { out << permy[idx] << std::endl; } out << "/\n\n"; - + out << "PERMZ\n\n"; for (size_t idx=0; idx < (size_t)permz.size(); ++idx) { out << permz[idx] << std::endl; } out << "/\n"; - + out.close(); } catch (const std::exception &e) { std::cerr << "Program threw an exception: " << e.what() << "\n"; throw; } - - diff --git a/examples/exp_variogram.cpp b/examples/exp_variogram.cpp index e6da314..9bae571 100644 --- a/examples/exp_variogram.cpp +++ b/examples/exp_variogram.cpp @@ -19,7 +19,7 @@ */ /* - This program computes data for an experimental + This program computes data for an experimental variogram from a cornerpoint geometry with properties. This works by choosing pairs of volumes chosen @@ -94,7 +94,7 @@ try double residual_tolerance = param.getDefault("residual_tolerance", 1e-8); int linsolver_verbosity = param.getDefault("linsolver_verbosity", 0); int linsolver_type = param.getDefault("linsolver_type", 1); - + // Check for unused parameters (potential typos). if (param.anyUnused()) { std::cout << "***** WARNING: Unused parameters: *****\n"; @@ -115,7 +115,7 @@ try exit(1); } - // Check user supplied variogram direction, either horizontal or vertical + // Check user supplied variogram direction, either horizontal or vertical enum variogram_directions { undefined, horizontal, vertical }; variogram_directions variogram_direction = undefined; std::string distancemetric; @@ -133,17 +133,17 @@ try std::cerr << "Error: variogram direction is undefined, user supplied '" << direction << "'.\n"; exit(1); } - - - // Check that we do not have any user input + + + // Check that we do not have any user input // that goes outside the coordinates described in // the cornerpoint file (runtime-exception will be thrown in case of error) ch.verifyInscribedShoebox(imin, ilen, imax, jmin, jlen, jmax, zmin, zlen, zmax); - + std::mt19937 gen; - + // Seed the random number generators with the current time, unless specified on command line // Warning: Current code does not allow 0 for the seed!! if (userseed == 0) { @@ -152,7 +152,7 @@ try else { gen.seed(userseed); } - + // Note that end is included in interval for uniform_int. std::uniform_int_distribution<> disti(imin, imax - ilen); @@ -161,14 +161,14 @@ try auto ri = [&disti, &gen] { return disti(gen); }; auto rj = [&distj, &gen] { return distj(gen); }; auto rz = [&distz, &gen] { return distz(gen); }; - + // Storage for results std::vector distances; std::vector porodiffs; std::vector permxdiffs; std::vector permydiffs; std::vector permzdiffs; - + for (int pair = 1; pair <= pairs; ++pair) { int istart_1 = ri(); int jstart_1 = rj(); @@ -195,9 +195,9 @@ try istart_2 = istart_1; jstart_2 = jstart_1; zstart_2 = rz(); - } + } ch.chop(istart_2, istart_2 + ilen, jstart_2, jstart_2 + jlen, zstart_2, zstart_2 + zlen, false); - + auto subdeck_2 = ch.subDeck(); Opm::SinglePhaseUpscaler upscaler_2; upscaler_2.init(subdeck_2, Opm::SinglePhaseUpscaler::Fixed, minpermSI, @@ -217,16 +217,16 @@ try permydiffs.push_back(fabs(upscaled_K_2(1,1) - upscaled_K_1(1,1))); permzdiffs.push_back(fabs(upscaled_K_2(2,2) - upscaled_K_1(2,2))); } - + // Make stream of output data, to be outputted to screen and optionally to file std::stringstream outputtmp; - + outputtmp << "################################################################################################" << std::endl; outputtmp << "# Data for experimental variogram" << std::endl; outputtmp << "#" << std::endl; time_t now = time(NULL); outputtmp << "# Finished: " << asctime(localtime(&now)); - + utsname hostname; uname(&hostname); outputtmp << "# Hostname: " << hostname.nodename << std::endl; outputtmp << "#" << std::endl; @@ -239,7 +239,7 @@ try outputtmp << "# pairs: " << pairs << std::endl; outputtmp << "################################################################################################" << std::endl; outputtmp << "# distance (" << distancemetric << ") porositydiff permxdiff permydiff permzdiff" << std::endl; - + const int fieldwidth = outputprecision + 8; for (int pair = 1; pair <= pairs; ++pair) { outputtmp << std::showpoint << std::setw(fieldwidth) << std::setprecision(outputprecision) << distances[pair-1] << '\t' << @@ -249,21 +249,20 @@ try std::showpoint << std::setw(fieldwidth) << std::setprecision(outputprecision) << permzdiffs[pair-1] << '\t' << std::endl; } - + if (resultfile != "") { std::cout << "Writing results to " << resultfile << std::endl; std::ofstream outfile; outfile.open(resultfile.c_str(), std::ios::out | std::ios::trunc); outfile << outputtmp.str(); - outfile.close(); + outfile.close(); } - - - + + + std::cout << outputtmp.str(); } catch (const std::exception &e) { std::cerr << "Program threw an exception: " << e.what() << "\n"; throw; } - diff --git a/examples/grdecldips.cpp b/examples/grdecldips.cpp index 894c89b..c8e62a1 100644 --- a/examples/grdecldips.cpp +++ b/examples/grdecldips.cpp @@ -52,17 +52,17 @@ using namespace std; -int main(int argc, char** argv) try { - if (argc == 1) { // If no arguments supplied +int main(int argc, char** argv) try { + if (argc == 1) { // If no arguments supplied cout << "Usage: grdecldips gridfilename=foo.grdecl [mincellvolume=1e-8] " << endl; cout << " [listallcells=false] [output=filename.txt]" << endl; exit(1); - } + } Dune::MPIHelper::instance(argc, argv); Opm::ParameterGroup param(argc, argv); - + std::string gridfilename = param.get("gridfilename"); double minCellVolume = param.getDefault("mincellvolume", 1e-8); bool listallcells = param.getDefault("listallcells", false); @@ -76,13 +76,13 @@ int main(int argc, char** argv) try { auto deck = Opm::RelPermUpscaleHelper::parseEclipseFile(gridfilename); Opm::EclipseGridInspector gridinspector(deck); - + // Check that we have the information we need from the eclipse file, we will check PERM-fields later - if (! (deck.hasKeyword("SPECGRID") && deck.hasKeyword("COORD") && deck.hasKeyword("ZCORN"))) { - cerr << "Error: Did not find SPECGRID, COORD and ZCORN in Eclipse file " << gridfilename << endl; - exit(1); + if (! (deck.hasKeyword("SPECGRID") && deck.hasKeyword("COORD") && deck.hasKeyword("ZCORN"))) { + cerr << "Error: Did not find SPECGRID, COORD and ZCORN in Eclipse file " << gridfilename << endl; + exit(1); } - + /*************************** * Find dips for every cell. */ @@ -94,11 +94,11 @@ int main(int argc, char** argv) try { griddims[2] = specgridRecord.getItem("NZ").get< int >(0); vector xdips, ydips, cellvolumes; vector cellidxs_i, cellidxs_j, cellidxs_k; - - int ignoredCellCount = 0; + + int ignoredCellCount = 0; for (int k=0; k < griddims[2]; ++k) { for (int j=0; j < griddims[1]; ++j) { - for (int i=0; i < griddims[0]; ++i) { + for (int i=0; i < griddims[0]; ++i) { double cellVolume = gridinspector.cellVolumeVerticalPillars(i, j, k); if (cellVolume > minCellVolume) { std::pair xydip = gridinspector.cellDips(i, j, k); @@ -115,13 +115,13 @@ int main(int argc, char** argv) try { } } } - - // Average xdips and ydips + + // Average xdips and ydips double xdipaverage = accumulate(xdips.begin(), xdips.end(), 0.0)/xdips.size(); double ydipaverage = accumulate(ydips.begin(), ydips.end(), 0.0)/ydips.size(); - + stringstream outputtmp; - + // Print a table of all computed values: outputtmp << "###############################################################################" << endl; outputtmp << "# Results from upscaling dips."<< endl; @@ -129,7 +129,7 @@ int main(int argc, char** argv) try { time_t now = time(NULL); outputtmp << "# Finished: " << asctime(localtime(&now)); utsname hostname; uname(&hostname); - + outputtmp << "#" << endl; outputtmp << "# Eclipse file: " << gridfilename << endl; outputtmp << "#" << endl; @@ -143,13 +143,13 @@ int main(int argc, char** argv) try { outputtmp << "# i j k xdip ydip cellvolume" << endl; for (unsigned int i=0; i < xdips.size(); ++i) { outputtmp << cellidxs_i[i] << " " << cellidxs_j[i] << " " << cellidxs_k[i] << "\t" << xdips[i] << "\t" << ydips[i] << "\t" << cellvolumes[i] << endl; - } + } } - else { + else { outputtmp << "x_dip_average " << xdipaverage << endl; outputtmp << "y_dip_average " << ydipaverage << endl; - } - + } + cout << endl << outputtmp.str(); if (outputfilename != "") { cout << "Writing results to " << outputfilename << endl; @@ -159,10 +159,9 @@ int main(int argc, char** argv) try { outfile.close(); } return 0; - + } catch (const std::exception &e) { std::cerr << "Program threw an exception: " << e.what() << "\n"; throw; } - diff --git a/examples/implicitcap_test.cpp b/examples/implicitcap_test.cpp index 9ab9e73..e7faf67 100644 --- a/examples/implicitcap_test.cpp +++ b/examples/implicitcap_test.cpp @@ -81,5 +81,3 @@ catch (const std::exception &e) { std::cerr << "Program threw an exception: " << e.what() << "\n"; throw; } - - diff --git a/examples/mimetic_periodic_test.cpp b/examples/mimetic_periodic_test.cpp index a4de88e..76072be 100644 --- a/examples/mimetic_periodic_test.cpp +++ b/examples/mimetic_periodic_test.cpp @@ -109,4 +109,3 @@ catch (const std::exception &e) { std::cerr << "Program threw an exception: " << e.what() << "\n"; throw; } - diff --git a/examples/mimetic_solver_test.cpp b/examples/mimetic_solver_test.cpp index fe87fea..f523061 100644 --- a/examples/mimetic_solver_test.cpp +++ b/examples/mimetic_solver_test.cpp @@ -135,7 +135,7 @@ void test_flowsolver(const GI& g, const RI& r) vtkwriter.addCellData(cell_velocity_flat, "velocity", dim); vtkwriter.addCellData(cell_pressure, "pressure"); vtkwriter.write("testsolution-" + std::to_string(0), Dune::VTK::ascii); -#else +#else solver.printSystem("system"); typedef typename FlowSolver::SolutionType FlowSolution; FlowSolution soln = solver.getSolution(); diff --git a/examples/sim_steadystate_implicit.cpp b/examples/sim_steadystate_implicit.cpp index 7a5a6cf..bf8c87c 100644 --- a/examples/sim_steadystate_implicit.cpp +++ b/examples/sim_steadystate_implicit.cpp @@ -67,5 +67,3 @@ catch (const std::exception &e) { std::cerr << "Program threw an exception: " << e.what() << "\n"; throw; } - - diff --git a/examples/upscale_avg.cpp b/examples/upscale_avg.cpp index 94d6b4c..24badad 100644 --- a/examples/upscale_avg.cpp +++ b/examples/upscale_avg.cpp @@ -19,22 +19,22 @@ /** @file upscale_avg.C * @brief Upscales using simple averages and reports statistics - * - * Processes permeability and porosity for a given Eclipse-model. - * + * + * Processes permeability and porosity for a given Eclipse-model. + * * Input is Eclipse grid format specifying the corner-point * grid (must be of shoebox-shape, but this condition is slightly relaxed * on top and bottom surfaces). - * + * * The input eclipse file must specify the permeability properties * for each cell. - * + * * The grid processing step from the permeability upscaling code is * used as a simple code library for computing the volume of each * cell, if some effort is invested, this can be circumventing by a * simpler code that does nothing but interpret the coordinates and * calculates the volume, which will be much faster. - * + * * */ #include @@ -85,7 +85,7 @@ namespace { /** @brief Computes simple statistics. */ -int main(int varnum, char** vararg) try { +int main(int varnum, char** vararg) try { Dune::MPIHelper::instance(varnum, vararg); @@ -97,7 +97,7 @@ int main(int varnum, char** vararg) try { cout << "Error: No eclipsefile provided" << endl; usage(); exit(1); - } + } /* @@ -108,14 +108,14 @@ int main(int varnum, char** vararg) try { //Default is not to use actnum for this //When use_actnum is used the cell volumes of the cells with actnum 0 is set to 0 - /* Loop over all command line options in order to look - for options. + /* Loop over all command line options in order to look + for options. argidx loops over all the arguments here, and updates the variable 'argeclindex' *if* it finds any legal options, 'argeclindex' is so that vararg[argeclindex] = the eclipse - filename. If options are illegal, argeclindex will be wrong, - + filename. If options are illegal, argeclindex will be wrong, + */ int argeclindex = 0; for (int argidx = 1; argidx < varnum; argidx += 2) { @@ -132,11 +132,11 @@ int main(int varnum, char** vararg) try { usageandexit(); } } - else { - // if vararg[argidx] does not start in '-', + 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, + break; // out of for-loop, } } @@ -150,11 +150,11 @@ int main(int varnum, char** vararg) try { exit(1); } eclipsefile.close(); - + // Variables for timing/profiling clock_t start, finish; double timeused = 0; - + /*********************************************************************** * Step X * Load geometry and data from Eclipse file @@ -163,19 +163,19 @@ int main(int varnum, char** vararg) try { flush(cout); start = clock(); // eclParser_p is here a pointer to an object of type Opm::EclipseGridParser // (this pointer trick is necessary for the try-catch-clause to work) - + auto deck = RelPermUpscaleHelper::parseEclipseFile(ECLIPSEFILENAME); Opm::EclipseGridInspector eclInspector(deck); finish = clock(); timeused = (double(finish)-double(start))/CLOCKS_PER_SEC; cout << " (" << timeused <<" secs)" << endl; - + // Check that we have the information we need from the eclipse file, we will check PERM-fields later - if (! (deck.hasKeyword("SPECGRID") && deck.hasKeyword("COORD") && deck.hasKeyword("ZCORN"))) { - cerr << "Error: Did not find SPECGRID, COORD and ZCORN in Eclipse file " << ECLIPSEFILENAME << endl; - usage(); - exit(1); + if (! (deck.hasKeyword("SPECGRID") && deck.hasKeyword("COORD") && deck.hasKeyword("ZCORN"))) { + cerr << "Error: Did not find SPECGRID, COORD and ZCORN in Eclipse file " << ECLIPSEFILENAME << endl; + usage(); + exit(1); } SinglePhaseUpscaler upscaler; @@ -195,12 +195,12 @@ int main(int varnum, char** vararg) try { else { cout << "Error: option use_actnum set to 1 but Gridfile " << ECLIPSEFILENAME << " does not include any field ACTNUM." << endl; usage(); - exit(1); + exit(1); } } // Print header for the output - + cout << "Statistics for filename: " << ECLIPSEFILENAME << endl; cout << "-----------------------------------------------------" << endl; bool doporosity = false; @@ -211,7 +211,7 @@ int main(int varnum, char** vararg) try { if (deck.hasKeyword("NTG")) { // Ntg only used together with PORO if (deck.hasKeyword("PORO")) dontg = true; - } + } bool doperm = false; if (deck.hasKeyword("PERMX")) { @@ -237,9 +237,9 @@ int main(int varnum, char** vararg) try { griddims[1] << " x " << griddims[2] << ")" << endl; int pillars = (griddims[0]+1) * (griddims[1]+1); - cout << " Pillars: " << pillars << " (" << griddims[0]+1 << + cout << " Pillars: " << pillars << " (" << griddims[0]+1 << " x " << griddims[1]+1 << ")" << endl; - + // Find max and min in x-, y- and z-directions std::array gridlimits = eclInspector.getGridLimits(); cout << " x-limits: " << gridlimits[0] << " -- " << gridlimits[1] << endl; @@ -267,13 +267,13 @@ int main(int varnum, char** vararg) try { permys = deck["PERMY"].back().getRawDoubleData(); permzs = deck["PERMZ"].back().getRawDoubleData(); } - + } const std::vector& ecl_idx = upscaler.grid().globalCell(); Dune::CpGrid::Codim<0>::LeafIterator c = upscaler.grid().leafbegin<0>(); for (; c != upscaler.grid().leafend<0>(); ++c) { size_t cell_idx = ecl_idx[c->index()]; - + //for (size_t cell_idx = 0; cell_idx < num_eclipse_cells; ++cell_idx) { if (!use_actnum){ cellVolumes[cell_idx] = c->geometry().volume(); @@ -292,19 +292,19 @@ int main(int varnum, char** vararg) try { } } } - + cout << " Active cells: " << active_cell_count << " (" << (double)active_cell_count/(double)num_eclipse_cells*100.0 << "%)" << endl; double volume = std::accumulate(cellVolumes.begin(), cellVolumes.end(), 0.0); cout << " Total volume: " << volume << endl; if (doporosity) { - double poreVolume = std::accumulate(cellPoreVolumes.begin(), + double poreVolume = std::accumulate(cellPoreVolumes.begin(), cellPoreVolumes.end(), 0.0); cout << " Total porevolume: " << poreVolume << endl; cout << " Upscaled porosity: " << poreVolume/volume << endl; - + int zeroporocells = 0; int negativeporocells = 0; for (size_t cell_idx = 0; cell_idx < (size_t)num_eclipse_cells; ++cell_idx) { @@ -327,20 +327,20 @@ int main(int varnum, char** vararg) try { netCellVolumes.end(), 0.0); cout << " Total net volume: " << netVolume << endl; - cout << " Upscaled NTG: " << netVolume/volume << endl; - double netPoreVolume = std::accumulate(netCellPoreVolumes.begin(), + cout << " Upscaled NTG: " << netVolume/volume << endl; + double netPoreVolume = std::accumulate(netCellPoreVolumes.begin(), netCellPoreVolumes.end(), 0.0); cout << " Total net porevolume: " << netPoreVolume << endl; - cout << " Upscaled net porosity: " << netPoreVolume/netVolume << endl; + cout << " Upscaled net porosity: " << netPoreVolume/netVolume << endl; } - + double permxsum = 0.0, permysum = 0.0, permzsum = 0.0; double invpermxsum = 0.0, invpermysum = 0.0, invpermzsum = 0.0; double volpermxsum = 0.0, volpermysum = 0.0, volpermzsum = 0.0; double invvolpermxsum = 0.0, invvolpermysum = 0.0, invvolpermzsum = 0.0; double logpermxsum = 0.0, logpermysum = 0.0, logpermzsum = 0.0; - double logvolpermxsum = 0.0, logvolpermysum = 0.0, logvolpermzsum=0.0; + double logvolpermxsum = 0.0, logvolpermysum = 0.0, logvolpermzsum=0.0; if (doperm) { int zeropermcells = 0; int negativepermcells = 0; @@ -365,7 +365,7 @@ int main(int varnum, char** vararg) try { invvolpermysum += cellVolumes[cell_idx] / permys[cell_idx]; logpermysum += log(permys[cell_idx]); logvolpermysum += log(permys[cell_idx]) * cellVolumes[cell_idx]; - + permzsum += permzs[cell_idx]; volpermzsum += permzs[cell_idx] * cellVolumes[cell_idx]; invpermzsum += 1.0/permzs[cell_idx]; @@ -375,14 +375,14 @@ int main(int varnum, char** vararg) try { } } } - + cout << "Total arithmetic permeability average: " << volpermxsum/volume << endl; cout << " Total harmonic permeability average: " << volume/invvolpermxsum << endl; cout << " Total geometric permeability average: " << exp(logvolpermxsum/volume) << endl; cout << "Total arithmetic permeability average: " << permxsum/((double)active_cell_count) << " (not volume-weighted)" << endl; cout << " Total harmonic permeability average: " << ((double)active_cell_count)/invpermxsum << " (not volume-weighted)" << endl; cout << " Total geometric permeability average: " << exp(logpermxsum/(double)active_cell_count) << " (not volume-weighted)" << endl; - + if (anisotropic_input) { cout << endl; cout << "Total arithmetic permeability (y) average: " << volpermysum/volume << endl; @@ -398,7 +398,7 @@ int main(int varnum, char** vararg) try { cout << "Total arithmetic permeability (z) average: " << permzsum/((double)active_cell_count) << " (not volume-weighted)" << endl; cout << " Total harmonic permeability (z) average: " << ((double)active_cell_count)/invpermzsum << " (not volume-weighted)" << endl; cout << " Total geometric permeability (z) average: " << exp(logpermzsum/(double)active_cell_count) << " (not volume-weighted)" << endl; - + } if (zeropermcells > 0) { @@ -409,13 +409,13 @@ int main(int varnum, char** vararg) try { } } - // Then do statistics on rocktype by rocktype basis + // Then do statistics on rocktype by rocktype basis bool dosatnums = false; vector satnums; if (deck.hasKeyword("SATNUM")) { dosatnums = true; satnums = deck["SATNUM"].back().getIntData(); - } // If SATNUM was not present, maybe ROCKTYPE is there, + } // If SATNUM was not present, maybe ROCKTYPE is there, // if so, we will use it as SATNUM. else if (deck.hasKeyword("ROCKTYPE")) { dosatnums = true; @@ -431,7 +431,7 @@ int main(int varnum, char** vararg) try { if (satnums[i] > maxsatnumvalue) { maxsatnumvalue = satnums[i]; } - if (satnums[i] < 0 || satnums[i] > 1000) { + if (satnums[i] < 0 || satnums[i] > 1000) { cerr << "satnums[" << i << "] = " << satnums[i] << ", not sane, quitting." << endl; exit(1); } @@ -439,14 +439,14 @@ int main(int varnum, char** vararg) try { vector permxsum_rocktype; permxsum_rocktype.resize(maxsatnumvalue+1, 0.0); - + vector invpermxsum_rocktype; invpermxsum_rocktype.resize(maxsatnumvalue+1, 0.0); // for harmonic average vector volpermxsum_rocktype; // volume weighted volpermxsum_rocktype.resize(maxsatnumvalue+1, 0.0); - - vector invvolpermxsum_rocktype; // volume weighted + + vector invvolpermxsum_rocktype; // volume weighted invvolpermxsum_rocktype.resize(maxsatnumvalue+1, 0.0); // for harmonic average vector porevolumesum_rocktype; @@ -476,7 +476,7 @@ int main(int varnum, char** vararg) try { volpermxsum_rocktype[satnums[cell_idx]] += permxs[cell_idx] * cellVolumes[cell_idx]; invvolpermxsum_rocktype[satnums[cell_idx]] += cellVolumes[cell_idx] / permxs[cell_idx]; - + } if (doporosity) { porevolumesum_rocktype[satnums[cell_idx]] += cellVolumes[cell_idx] * poros[cell_idx]; @@ -484,7 +484,7 @@ int main(int varnum, char** vararg) try { } } - // Compute the sample variance in porosity per rock type + // Compute the sample variance in porosity per rock type /* phi_avg = 1/V*sum(v_i*phi_i) phi_var = 1/V*sum(v_i*(phi_i-phi_avg)^2) @@ -492,23 +492,23 @@ int main(int varnum, char** vararg) try { if (doporosity) { for (size_t cell_idx = 0; cell_idx < (size_t)num_eclipse_cells; ++cell_idx) { if (cellVolumes[cell_idx] > emptycellvolumecutoff) { - porosityvariancesum_rocktype[satnums[cell_idx]] += cellVolumes[cell_idx] + porosityvariancesum_rocktype[satnums[cell_idx]] += cellVolumes[cell_idx] * pow((poros[cell_idx]-porevolumesum_rocktype[satnums[cell_idx]]/volumesum_rocktype[satnums[cell_idx]]),2); - + } - } - } - + } + } + // Now loop over rocktypes in order to print statistics // Does SATNUM/ROCKTYPE start with 0 or 1? for (int rocktype_idx = 0; rocktype_idx <= maxsatnumvalue; ++rocktype_idx) { if (activecellcount_rocktype[rocktype_idx] > 0) { cout << endl << "Statistics for rocktype " << rocktype_idx << endl; - + cout << " Volume: " << volumesum_rocktype[rocktype_idx] << " (" << volumesum_rocktype[rocktype_idx]/volume*100 << "%)" << endl; cout << " Total cells: " << totalcellcount_rocktype[rocktype_idx] << " (" << (double)totalcellcount_rocktype[rocktype_idx]/(double)num_eclipse_cells*100 << "%)" << endl; - cout << " Active cells: " << activecellcount_rocktype[rocktype_idx] << " (" << + cout << " Active cells: " << activecellcount_rocktype[rocktype_idx] << " (" << (double)activecellcount_rocktype[rocktype_idx]/(double)totalcellcount_rocktype[rocktype_idx]*100.0 << "% within rocktype)" << endl; if (doporosity) { cout << " Porevolume: " << porevolumesum_rocktype[rocktype_idx] << endl; @@ -520,17 +520,16 @@ int main(int varnum, char** vararg) try { cout << " Perm harm. avg: " << volumesum_rocktype[rocktype_idx]/invvolpermxsum_rocktype[rocktype_idx] << endl; cout << "Perm arith. avg: " << permxsum_rocktype[rocktype_idx]/activecellcount_rocktype[rocktype_idx] << " (not volume-weighted)" << endl; cout << " Perm harm. avg: " << activecellcount_rocktype[rocktype_idx]/invpermxsum_rocktype[rocktype_idx] << " (not volume-weighted)" << endl; - } + } } } } - - + + return 0; - + } catch (const std::exception &e) { std::cerr << "Program threw an exception: " << e.what() << "\n"; throw; } - diff --git a/examples/upscale_cap.cpp b/examples/upscale_cap.cpp index 2e66bc2..d726755 100644 --- a/examples/upscale_cap.cpp +++ b/examples/upscale_cap.cpp @@ -22,8 +22,8 @@ @file upscale_cap.C @brief Upscales capillary pressure function, water saturation, and water saturation pr. rocktype. - Description: - + Description: + Reads in a lithofacies geometry in Eclipse format, reads in J(S_w) permeability and porosity, and upscale the capillary pressure function and water saturation, both for the whole geometry, and on a per rocktype basis @@ -39,17 +39,17 @@ 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 + 2: Read Eclipse file 3: Read and J-function for each stone-type. - 4: Find minimum and maximum capillary pressure from the + 4: Find minimum and maximum capillary pressure from the J-functions in each cell. 5: Upscale water saturation as a function of capillary pressure 6: Print output to screen and optionally to file. - + @author: Håvard Berland statoil.com> */ #include @@ -125,7 +125,7 @@ namespace { int main(int varnum, char** vararg) try -{ +{ /****************************************************************************** * Step 1: @@ -145,9 +145,9 @@ try map options; options.insert(make_pair("points", "50" )); // Number of saturation points (uniformly distributed within saturation endpoints) options.insert(make_pair("jFunctionCurve", "4")); // Which column in the rock type file is the J-function curve - options.insert(make_pair("output", "")); // If this is set, output goes to screen and to this file. + options.insert(make_pair("output", "")); // If this is set, output goes to screen and to this file. options.insert(make_pair("outputprecision", "8")); // number of decimals to print - options.insert(make_pair("surfaceTension", "11")); // Surface tension given in dynes/cm + options.insert(make_pair("surfaceTension", "11")); // Surface tension given in dynes/cm options.insert(make_pair("maxPermContrast", "1e7")); // maximum allowed contrast in each single-phase computation options.insert(make_pair("minPerm", "1e-12")); // absoluted minimum allowed minimal cell permeability options.insert(make_pair("minPoro", "0.0001")); // this limit is necessary for pcmin/max computation @@ -160,14 +160,14 @@ try // Reference: http://www.spe.org/spe-site/spe/spe/papers/authors/Metric_Standard.pdf - /* Loop over all command line options in order to look - for options. + /* Loop over all command line options in order to look + for options. argidx loops over all the arguments here, and updates the variable 'argeclindex' *if* it finds any legal options, 'argeclindex' is so that vararg[argeclindex] = the eclipse - filename. If options are illegal, argeclindex will be wrong, - + filename. If options are illegal, argeclindex will be wrong, + */ int argeclindex = 0; for (int argidx = 1; argidx < varnum; argidx += 2) { @@ -184,59 +184,59 @@ try usageandexit(); } } - else { - // if vararg[argidx] does not start in '-', + 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, + break; // out of for-loop, } } - + // 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 JFindex = argeclindex; - + /* Check if at least one J-function is supplied on command line */ if (varnum <= JFindex) { cerr << "Error: No J-functions found on command line." << endl; usageandexit(); } - - + + /*********************************************************************** * Step 2: * Load geometry and data from Eclipse file */ - - // Read data from the Eclipse file and + + // Read data from the Eclipse file and // populate our vectors with data from the file const double emptycellvolumecutoff = 1e-10; - + // Test if filename exists and is readable ifstream eclipsefile(ECLIPSEFILENAME, ios::in); if (eclipsefile.fail()) { cerr << "Error: Filename " << ECLIPSEFILENAME << " not found or not readable." << endl; usageandexit(); } - eclipsefile.close(); + eclipsefile.close(); cout << "Parsing Eclipse file <" << ECLIPSEFILENAME << "> ... " << endl; auto deck = RelPermUpscaleHelper::parseEclipseFile(ECLIPSEFILENAME); - - // Check that we have the information we need from the eclipse file: - if (! (deck.hasKeyword("SPECGRID") && deck.hasKeyword("COORD") && deck.hasKeyword("ZCORN") - && deck.hasKeyword("PORO") && deck.hasKeyword("PERMX") && deck.hasKeyword("SATNUM"))) { - cerr << "Error: Did not find SPECGRID, COORD and ZCORN in Eclipse file " << ECLIPSEFILENAME << endl; - usageandexit(); - } - - vector satnums = deck["SATNUM"].back().getIntData(); + + // Check that we have the information we need from the eclipse file: + if (! (deck.hasKeyword("SPECGRID") && deck.hasKeyword("COORD") && deck.hasKeyword("ZCORN") + && deck.hasKeyword("PORO") && deck.hasKeyword("PERMX") && deck.hasKeyword("SATNUM"))) { + cerr << "Error: Did not find SPECGRID, COORD and ZCORN in Eclipse file " << ECLIPSEFILENAME << endl; + usageandexit(); + } + + vector satnums = deck["SATNUM"].back().getIntData(); vector poros = deck["PORO"].back().getRawDoubleData(); vector permxs = deck["PERMX"].back().getRawDoubleData(); vector griddims(3); @@ -254,12 +254,12 @@ try - Check that SATNUM are set sensibly, that is => 0 and < 1000, error if not. - Check that porosity is between 0 and 1, error if not. Set to minPoro if zero or less than minPoro (due to pcmin/max computation) - - Check that permeability is zero or positive. Error if negative. + - Check that permeability is zero or positive. Error if negative. Set to minPerm if zero or less than minPerm. - Check maximum number of SATNUM values (can be number of rock types present) */ for (unsigned int i = 0; i < satnums.size(); ++i) { - if (satnums[i] < 0 || satnums[i] > 1000) { + if (satnums[i] < 0 || satnums[i] > 1000) { cerr << "satnums[" << i << "] = " << satnums[i] << ", not sane, quitting." << endl; usageandexit(); } @@ -283,40 +283,40 @@ try // Explicitly handle "no rock" cells, set them to minimum perm and zero porosity. if (satnums[i] == 0) { permxs[i] = minPerm; - poros[i] = 0; // zero poro is fine for these cells, as they are not + poros[i] = 0; // zero poro is fine for these cells, as they are not // used in pcmin/max computation. } - } + } /*************************************************************************** * 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 + * We read columns from text-files, syntax allowed is determined + * by MonotCubicInterpolator which actually opens and parses the * text files. */ // 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. int stone_types = int(*(max_element(satnums.begin(), satnums.end()))); std::vector InvJfunctions; // Holds the inverse of the loaded J-functions. - + std::vector JfunctionNames; // Placeholder for the names of the loaded J-functions. // Input for surfaceTension is dynes/cm // SI units are Joules/square metre - const double surfaceTension = atof(options["surfaceTension"].c_str()) * 1e-3; // multiply with 10^-3 to obtain SI units + const double surfaceTension = atof(options["surfaceTension"].c_str()) * 1e-3; // multiply with 10^-3 to obtain SI units const int jFunctionCurve = atoi(options["jFunctionCurve"].c_str()); const int interpolationPoints = atoi(options["points"].c_str()); 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 == JFindex + stone_types) { for (int i=0 ; i < stone_types; ++i) { const char* ROCKFILENAME = vararg[JFindex+i]; @@ -326,10 +326,10 @@ try cerr << "Error: Filename " << ROCKFILENAME << " not found or not readable." << endl; usageandexit(); } - rockfile.close(); + rockfile.close(); MonotCubicInterpolator Jtmp; try { - Jtmp = MonotCubicInterpolator(ROCKFILENAME, 1, jFunctionCurve); + Jtmp = MonotCubicInterpolator(ROCKFILENAME, 1, jFunctionCurve); } catch (const char * errormessage) { cerr << "Error: " << errormessage << endl; @@ -345,9 +345,9 @@ try cerr << "Error: Jfunction " << i+1 << " in rock file " << ROCKFILENAME << " was not invertible." << endl; usageandexit(); } - } + } } - + else if (varnum == JFindex + 1) { for (int i=0; i < stone_types; ++i) { const char* ROCKFILENAME = vararg[JFindex]; @@ -357,7 +357,7 @@ try cerr << "Error: Filename " << ROCKFILENAME << " not found or not readable." << endl; usageandexit(); } - rockfile.close(); + rockfile.close(); MonotCubicInterpolator Jtmp; try { Jtmp = MonotCubicInterpolator(ROCKFILENAME, 1, jFunctionCurve); @@ -382,7 +382,7 @@ try cerr << "Error: Wrong number of stone-functions provided. " << endl; usageandexit(); } - + /****************************************************************************** * Step 5: @@ -398,7 +398,7 @@ try cout << "Illegal contrast value" << endl; usageandexit(); } - + // Construct an object for single-phase upscaling, since we need to get some // information from the grid. @@ -407,7 +407,7 @@ try Opm::unit::convert::from(minPerm, Opm::prefix::milli*Opm::unit::darcy), 1e-8, 0, 1, false); // options on this line are noops for upscale_cap - vector cellVolumes, cellPoreVolumes; + vector cellVolumes, cellPoreVolumes; cellVolumes.resize(satnums.size(), 0.0); cellPoreVolumes.resize(satnums.size(), 0.0); @@ -422,12 +422,12 @@ try /* Find minimium and maximum capillary pressure values in each cell, and use the global min/max as the two initial pressure points for computations. - - Also find max single-phase permeability, used to obey the + + Also find max single-phase permeability, used to obey the maxPermContrast option. Also find properly upscaled saturation endpoints, these are - printed out to stdout for reference during computations, but will + printed out to stdout for reference during computations, but will automatically appear as the lowest and highest saturation points in finished output. */ @@ -446,34 +446,34 @@ try cellVolumes[cell_idx] = c->geometry().volume(); cellPoreVolumes[cell_idx] = cellVolumes[cell_idx] * poros[cell_idx]; - + double Pcmincandidate = InvJfunctions[int(satnums[cell_idx])-1].getMinimumX().first / sqrt(permxs[cell_idx] * milliDarcyToSqMetre/poros[cell_idx]) * surfaceTension; Pcmin = min(Pcmincandidate, Pcmin); - + double Pcmaxcandidate = InvJfunctions[int(satnums[cell_idx])-1].getMaximumX().first / sqrt(permxs[cell_idx] * milliDarcyToSqMetre/poros[cell_idx]) * surfaceTension; Pcmax = max(Pcmaxcandidate, Pcmax); - + maxSinglePhasePerm = max( maxSinglePhasePerm, permxs[cell_idx]); - + cellporevolume_rocktype[satnums[cell_idx]] += cellPoreVolumes[cell_idx]; double minSw = InvJfunctions[int(satnums[cell_idx])-1].getMinimumF().second; double maxSw = InvJfunctions[int(satnums[cell_idx])-1].getMaximumF().second; - + // cout << "minSwc: " << minSw << endl; // cout << "maxSwc: " << maxSw << endl; - + // Add irreducible water saturation volume Swirvolume += minSw * cellPoreVolumes[cell_idx]; Sworvolume += maxSw * cellPoreVolumes[cell_idx]; - + } ++tesselatedCells; // keep count. } //double minSinglePhasePerm = max(maxSinglePhasePerm/maxPermContrast, minPerm); - + cout << "Pcmin: " << Pcmin << endl; cout << "Pcmax: " << Pcmax << endl; @@ -483,7 +483,7 @@ try } // Total porevolume and total volume -> upscaled porosity: - double poreVolume = std::accumulate(cellPoreVolumes.begin(), + double poreVolume = std::accumulate(cellPoreVolumes.begin(), cellPoreVolumes.end(), 0.0); double volume = std::accumulate(cellVolumes.begin(), @@ -499,9 +499,9 @@ try cout << "Upscaled Swir: " << Swir << endl; cout << "Upscaled Swmax: " << Swor << endl; //Swor=1-Swmax - // Sometimes, if Swmax=1 or Swir=0 in the input tables, the upscaled + // Sometimes, if Swmax=1 or Swir=0 in the input tables, the upscaled // values can be a little bit larger (within machine precision) and - // the check below fails. Hence, check if these values are within the + // the check below fails. Hence, check if these values are within the // the [0 1] interval within some precision double linsolver_tolerance = atof(options["linsolver_tolerance"].c_str()); if (Swor > 1.0 && Swor - linsolver_tolerance < 1.0) { @@ -513,7 +513,7 @@ try if (Swir < 0.0 || Swir > 1.0 || Swor < 0.0 || Swor > 1.0) { cerr << "ERROR: Swir/Swor unsensible. Check your input. Exiting"; usageandexit(); - } + } /*************************************************************************** * Step 6: @@ -529,9 +529,9 @@ try */ MonotCubicInterpolator WaterSaturationVsCapPressure; - + double largestSaturationInterval = Swor-Swir; - + double Ptestvalue; vector watersaturation_rocktype; @@ -539,10 +539,10 @@ try MonotCubicInterpolator tmp; watersaturation_rocktype.push_back(tmp); } - + while (largestSaturationInterval > (Swor-Swir)/200.0) { if (Pcmax == Pcmin) { - // This is a dummy situation, we go through once and then + // This is a dummy situation, we go through once and then // we are finished (this will be triggered by zero permeability) Ptestvalue = Pcmin; largestSaturationInterval = 0; @@ -565,7 +565,7 @@ try Ptestvalue = SatDiff.first; largestSaturationInterval = SatDiff.second; } - + // Check for saneness of Ptestvalue: if (std::isnan(Ptestvalue) || std::isinf(Ptestvalue)) { cerr << "ERROR: Ptestvalue was inf or nan" << endl; @@ -577,7 +577,7 @@ try for (unsigned int satidx = 0; satidx <= maxSatnum; ++satidx) { watervolume_rocktype[satidx] = 0.0; } - + double waterVolume = 0.0; for (unsigned int cell_idx = 0; cell_idx < satnums.size(); ++cell_idx) { if (cellVolumes[cell_idx] > emptycellvolumecutoff) { @@ -585,15 +585,15 @@ try if (satnums[cell_idx] > 0) { // handle "no rock" cells with satnum zero double PtestvalueCell; PtestvalueCell = Ptestvalue; - double Jvalue = sqrt(permxs[cell_idx] * milliDarcyToSqMetre/poros[cell_idx]) + double Jvalue = sqrt(permxs[cell_idx] * milliDarcyToSqMetre/poros[cell_idx]) * PtestvalueCell / surfaceTension; //cout << "JvalueCell: " << Jvalue << endl; - waterSaturationCell + waterSaturationCell = InvJfunctions[int(satnums[cell_idx])-1].evaluate(Jvalue); } waterVolume += waterSaturationCell * cellPoreVolumes[cell_idx]; watervolume_rocktype[satnums[cell_idx]] += waterSaturationCell * cellPoreVolumes[cell_idx]; - + } } WaterSaturationVsCapPressure.addPair(Ptestvalue, waterVolume/poreVolume); @@ -614,7 +614,7 @@ try // Check if the saturation vs cap pressure curve is monotone // If not, it would have been a problem for upscale_relperm, but // it is not as critical here, so we only issue a warning - // (upscale_relperm solves this by issung chopFlatEndpoints and possibly shrinkFlatAreas, + // (upscale_relperm solves this by issung chopFlatEndpoints and possibly shrinkFlatAreas, // but this is trickier to implement in this code due to watersaturation_rocktype[satidx]) if (!WaterSaturationVsCapPressure.isStrictlyMonotone()) { { @@ -622,10 +622,10 @@ try cerr << " Unphysical input data?." << endl; } } - MonotCubicInterpolator CapPressureVsWaterSaturation(WaterSaturationVsCapPressure.get_fVector(), + MonotCubicInterpolator CapPressureVsWaterSaturation(WaterSaturationVsCapPressure.get_fVector(), WaterSaturationVsCapPressure.get_xVector()); - + /********************************************************************************* * Step 9 * @@ -633,9 +633,9 @@ try * file if the '-outputWater'-option and/or '-outputOil' has been set, as this option is an * empty string by default. */ - vector Pvalues = WaterSaturationVsCapPressure.get_xVector(); - vector Satvalues = WaterSaturationVsCapPressure.get_fVector(); - + vector Pvalues = WaterSaturationVsCapPressure.get_xVector(); + vector Satvalues = WaterSaturationVsCapPressure.get_fVector(); + vector > watersaturation_rocktype_values; vector tmp; watersaturation_rocktype_values.push_back(tmp); // dummy zero index element @@ -643,14 +643,14 @@ try watersaturation_rocktype_values.push_back(watersaturation_rocktype[satidx].get_fVector()); } stringstream outputtmp; - + // Print a table of all computed values: outputtmp << "######################################################################" << endl; outputtmp << "# Results from upscaling capillary pressure and water saturations."<< endl; outputtmp << "#" << endl; time_t now = std::time(NULL); outputtmp << "# Finished: " << asctime(localtime(&now)); - + utsname hostname; uname(&hostname); outputtmp << "# Hostname: " << hostname.nodename << endl; @@ -669,11 +669,11 @@ try outputtmp << "# jFunctionCurve: " << options["jFunctionCurve"] << endl; outputtmp << "# points: " << options["points"] << endl; outputtmp << "# maxPermContrast: " << options["maxPermContrast"] << endl; - outputtmp << "# surfaceTension: " << options["surfaceTension"] << endl; + outputtmp << "# surfaceTension: " << options["surfaceTension"] << endl; outputtmp << "######################################################################" << endl; - outputtmp << "# Pc (Pa) Sw Sw1 Sw2 Sw3 etc.." << endl; - - + outputtmp << "# Pc (Pa) Sw Sw1 Sw2 Sw3 etc.." << endl; + + // If user wants interpolated output, do monotone cubic interpolation // by modifying the data vectors that are to be printed if (interpolationPoints > 0) { @@ -717,24 +717,24 @@ try const int fieldwidth = outputprecision + 8; for (unsigned int i=0; i < Satvalues.size(); ++i) { - outputtmp << showpoint << setw(fieldwidth) << setprecision(outputprecision) << Pvalues[i]; - outputtmp << showpoint << setw(fieldwidth) << setprecision(outputprecision) << Satvalues[i]; - for (unsigned int satidx = 1; satidx <= maxSatnum; ++satidx) { - outputtmp << showpoint << setw(fieldwidth) << setprecision(outputprecision) - << watersaturation_rocktype_values[satidx][i]; - } - outputtmp << endl; - + outputtmp << showpoint << setw(fieldwidth) << setprecision(outputprecision) << Pvalues[i]; + outputtmp << showpoint << setw(fieldwidth) << setprecision(outputprecision) << Satvalues[i]; + for (unsigned int satidx = 1; satidx <= maxSatnum; ++satidx) { + outputtmp << showpoint << setw(fieldwidth) << setprecision(outputprecision) + << watersaturation_rocktype_values[satidx][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(); + outfile.close(); } @@ -744,4 +744,3 @@ catch (const std::exception &e) { std::cerr << "Program threw an exception: " << e.what() << "\n"; throw; } - diff --git a/examples/upscale_cond.cpp b/examples/upscale_cond.cpp index 476ad12..4220e49 100644 --- a/examples/upscale_cond.cpp +++ b/examples/upscale_cond.cpp @@ -20,25 +20,25 @@ /** @file upscale_cond.C @brief Upscales resistivity as a function of water saturation - - Description: - + + Description: + Reads in a lithofacies geometry in Eclipse format, reads in J(S_w) - for each stone type, and calculates upscaled + for each stone type, and calculates upscaled resistivity values for values of S_w. - + Steps in the code: - + 1: Process command line options. - 2: Read and parse Eclipse file + 2: Read and parse Eclipse file 3: Read a 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 conductivity ( = 1/resistivity ) for each saturation point 8: Print output to screen and optionally to file. - - The relative permeability computation is based on + + The relative permeability computation is based on - Capillary equilibrium, p_c is spatially invariant. Units handling: @@ -76,7 +76,7 @@ #include #include - + using namespace Opm; using namespace std; @@ -140,7 +140,7 @@ namespace { int main(int varnum, char** vararg) try -{ +{ // Variables used for timing/profiling: clock_t start, finish; double timeused = 0, timeused_tesselation = 0; @@ -157,7 +157,7 @@ try } - + /* Populate options-map with default values */ @@ -171,10 +171,10 @@ try options.insert(make_pair("output", "" )); // Output to file as well as screen if provided options.insert(make_pair("mudresistivity", "1.4" )); // Constant resistivity for mud options.insert(make_pair("mud1rocktype", "0" )); // First mud rock type. 0 if there is no mud types - options.insert(make_pair("mud2rocktype", "0" )); // Second mud rock type. + options.insert(make_pair("mud2rocktype", "0" )); // Second mud rock type. options.insert(make_pair("jFunctionCurve", "4" )); // column number of J-function values in input files options.insert(make_pair("surfaceTension", "11" )); // Surface tension given in dynes/cm - options.insert(make_pair("interpolate", "0" )); // default is not to interpolate + options.insert(make_pair("interpolate", "0" )); // default is not to interpolate options.insert(make_pair("minPerm", "1e-12" )); // minimum modelled permeability (for saturation distr) options.insert(make_pair("minPoro", "0.0001")); // this limit is necessary for pcmin/max computation @@ -208,20 +208,20 @@ try Opm::unit::convert::to(1.0*Opm::prefix::milli*Opm::unit::darcy, Opm::unit::square(Opm::unit::meter)); // Reference: http://www.spe.org/spe-site/spe/spe/papers/authors/Metric_Standard.pdf - + /* - Look for strings in args matching the entries in the options map, + Look for strings in args matching the entries in the options map, if found, replace default values with command line values. */ - /* Loop over all command line options in order to look - for options. + /* Loop over all command line options in order to look + for options. argidx loops over all the arguments here, and updates the variable 'argeclindex' *if* it finds any legal options, 'argeclindex' is so that vararg[argeclindex] = the eclipse - filename. If options are illegal, argeclindex will be wrong, - + filename. If options are illegal, argeclindex will be wrong, + */ int argeclindex = 0; for (int argidx = 1; argidx < varnum; argidx += 2) { @@ -238,33 +238,33 @@ try usageandexit(); } } - else { - // if vararg[argidx] does not start in '-', + 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, + break; // out of for-loop, } } - - + + // argeclindex should now point at the eclipse file static char* ECLIPSEFILENAME(vararg[argeclindex]); argeclindex += 1; // areglcindex jumps to next input argument, now it points to the first J-file - + // argcindex now points to the first J-function. This index is not // to be touched now. static int JFindex = argeclindex; - + /* Check if at least one J-function is supplied on command line */ if (varnum <= JFindex) { if (isMaster) cerr << "Error: No J-functions found on command line." << endl; usageandexit(); } - // Check that boundary conditions are valid , and make booleans + // Check that boundary conditions are valid , and make booleans // for boundary conditions. This allows more readable code later. bool isFixed = false, isLinear = false, isPeriodic = false; - SinglePhaseUpscaler::BoundaryConditionType boundaryCondition = SinglePhaseUpscaler::Fixed ; + SinglePhaseUpscaler::BoundaryConditionType boundaryCondition = SinglePhaseUpscaler::Fixed ; int tensorElementCount = 0; // Number of independent elements in resulting tensor if (options["bc"].substr(0,1) == "f") { @@ -287,8 +287,8 @@ try if (isMaster) cout << "Invalid boundary condition: " << options["bc"] << endl; usageandexit(); } - - // Read data from the Eclipse file and + + // Read data from the Eclipse file and // populate our vectors with data from the file // Test if filename exists and is readable ifstream eclipsefile(ECLIPSEFILENAME, ios::in); @@ -296,7 +296,7 @@ try if (isMaster) cerr << "Error: Filename " << ECLIPSEFILENAME << " not found or not readable." << endl; usageandexit(); } - eclipsefile.close(); + eclipsefile.close(); if (isMaster) cout << "Parsing Eclipse file <" << ECLIPSEFILENAME << "> ... "; flush(cout); start = clock(); @@ -318,13 +318,13 @@ try vector satnums = deck["SATNUM"].back().getIntData(); vector poros = deck["PORO"].back().getRawDoubleData(); vector permxs = deck["PERMX"].back().getRawDoubleData(); - const double minPerm = atof(options["minPerm"].c_str()); + const double minPerm = atof(options["minPerm"].c_str()); const double minPoro = atof(options["minPoro"].c_str()); // Read in J-functions for all stone-types. // 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. @@ -337,7 +337,7 @@ try const int jFunctionCurve = atoi(options["jFunctionCurve"].c_str()); // default 4 // Input for surfaceTension is dynes/cm // SI units are Joules/square metre - const double surfaceTension = atof(options["surfaceTension"].c_str()) * 1e-3; // multiply + const double surfaceTension = atof(options["surfaceTension"].c_str()) * 1e-3; // multiply if (varnum == JFindex + stone_types) { for (int i=0 ; i < stone_types; ++i) { const char* ROCKFILENAME = vararg[JFindex+i]; @@ -347,11 +347,11 @@ try if (isMaster) cerr << "Error: Filename " << ROCKFILENAME << " not found or not readable." << endl; usageandexit(); } - rockfile.close(); - + rockfile.close(); + MonotCubicInterpolator Jtmp; try { - Jtmp = MonotCubicInterpolator(vararg[JFindex + i], 1, jFunctionCurve); + Jtmp = MonotCubicInterpolator(vararg[JFindex + i], 1, jFunctionCurve); } catch (const char * errormessage) { if (isMaster) { @@ -360,7 +360,7 @@ try } usageandexit(); } - + // Invert J-function, now we get saturation as a function of pressure: if (Jtmp.isStrictlyMonotone()) { InvJfunctions.push_back(MonotCubicInterpolator(Jtmp.get_fVector(), Jtmp.get_xVector())); @@ -372,7 +372,7 @@ try if (isMaster) cerr << "Error: Jfunction " << i+1 << " in rock file " << ROCKFILENAME << " was not invertible." << endl; usageandexit(); } - } + } } else if (varnum == JFindex + 1) { for (int i=0; i < stone_types; ++i) { @@ -383,7 +383,7 @@ try if (isMaster) cerr << "Error: Filename " << ROCKFILENAME << " not found or not readable." << endl; usageandexit(); } - rockfile.close(); + rockfile.close(); MonotCubicInterpolator Jtmp; try { @@ -409,8 +409,8 @@ try } } else { - if (isMaster) cerr << "Error: Wrong number of J-functions provided. " << endl << - "Note that all input arguments after the eclipsefile " << endl << + if (isMaster) cerr << "Error: Wrong number of J-functions provided. " << endl << + "Note that all input arguments after the eclipsefile " << endl << "are interpreted as J-functions." << endl; usageandexit(); } @@ -418,10 +418,10 @@ try /* Step 4 - + Generate tesselated grid. Tesselation depends on boundary conditions. For periodic boundary conditions, the grid needs to be massaged slightly - (crop top and bottom). These modifications ruin the computations for + (crop top and bottom). These modifications ruin the computations for fixed and linear boundary conditions. */ const auto& specgridRecord = deck["SPECGRID"].back().getRecord(0); @@ -439,9 +439,9 @@ try linsolver_tolerance, linsolver_verbosity, linsolver_type, twodim_hack); finish = clock(); timeused_tesselation = (double(finish)-double(start))/CLOCKS_PER_SEC; - if (isMaster) cout << " (" << timeused_tesselation <<" secs)" << endl; - - + if (isMaster) cout << " (" << timeused_tesselation <<" secs)" << endl; + + int maxSatnum = 0; int tesselatedCells = 0; @@ -454,7 +454,7 @@ try - Check maximum number of SATNUM values (can be number of rock types present) */ for (size_t i = 0; i < satnums.size(); ++i) { - if (satnums[i] < 0 || satnums[i] > 1000) { + if (satnums[i] < 0 || satnums[i] > 1000) { if (isMaster) cerr << "satnums[" << i << "] = " << satnums[i] << ", not sane, quitting." << endl; usageandexit(); } @@ -493,7 +493,7 @@ try */ - vector cellVolumes, cellPoreVolumes; + vector cellVolumes, cellPoreVolumes; cellVolumes.resize(satnums.size(), 0.0); cellPoreVolumes.resize(satnums.size(), 0.0); @@ -506,19 +506,19 @@ try if (satnums[cell_idx] > 0) { // Satnum zero is "no rock" cellVolumes[cell_idx] = c->geometry().volume(); cellPoreVolumes[cell_idx] = cellVolumes[cell_idx] * poros[cell_idx]; - + double Pcmincandidate, Pcmaxcandidate, minSw, maxSw; - + Pcmincandidate = InvJfunctions[int(satnums[cell_idx])-1].getMinimumX().first / sqrt(permxs[cell_idx] * milliDarcyToSqMetre / poros[cell_idx]) * surfaceTension; Pcmaxcandidate = InvJfunctions[int(satnums[cell_idx])-1].getMaximumX().first / sqrt(permxs[cell_idx] * milliDarcyToSqMetre/poros[cell_idx]) * surfaceTension; minSw = InvJfunctions[int(satnums[cell_idx])-1].getMinimumF().second; maxSw = InvJfunctions[int(satnums[cell_idx])-1].getMaximumF().second; - + Pcmin = min(Pcmincandidate, Pcmin); Pcmax = max(Pcmaxcandidate, Pcmax); - + // Add irreducible water saturation volume Swirvolume += minSw * cellPoreVolumes[cell_idx]; Sworvolume += maxSw * cellPoreVolumes[cell_idx]; @@ -537,7 +537,7 @@ try } // Total porevolume and total volume -> upscaled porosity: - double poreVolume = std::accumulate(cellPoreVolumes.begin(), + double poreVolume = std::accumulate(cellPoreVolumes.begin(), cellPoreVolumes.end(), 0.0); double volume = std::accumulate(cellVolumes.begin(), @@ -555,9 +555,9 @@ try cout << "Saturation points to be computed: " << points << endl; } - // Sometimes, if Swmax=1 or Swir=0 in the input tables, the upscaled + // Sometimes, if Swmax=1 or Swir=0 in the input tables, the upscaled // values can be a little bit larger (within machine precision) and - // the check below fails. Hence, check if these values are within the + // the check below fails. Hence, check if these values are within the // the [0 1] interval within some precision (use linsolver_precision) if (Swor > 1.0 && Swor - linsolver_tolerance < 1.0) { Swor = 1.0; @@ -568,8 +568,8 @@ try if (Swir < 0.0 || Swir > 1.0 || Swor < 0.0 || Swor > 1.0) { if (isMaster) cerr << "ERROR: Swir/Swor unsensible. Check your input. Exiting"; usageandexit(); - } - + } + /************************************************************ * Handle petrophysical constants/parameters @@ -586,7 +586,7 @@ try vector cementationexponents, saturationexponents; cementationexponents.resize(maxSatnum + 1); // index 0 unused!! saturationexponents.resize(maxSatnum + 1); // index 0 unused!! - + /* The saturation and cementation exponents may be set globally for all rock types using the options cementationExponent and saturationExponent. However, these can also be set on a rock by @@ -597,7 +597,7 @@ try stringstream rocktypestring; rocktypestring << i; // integer to string conversion. - + //cout << "rock" + rocktypestring.str() + "cemexp" << endl; double cemoptionvalue = atof(options["rock" + rocktypestring.str() + "cemexp"].c_str()); double satoptionvalue = atof(options["rock" + rocktypestring.str() + "satexp"].c_str()); @@ -616,19 +616,19 @@ try //cout << "rocktype " << i << endl; //cout << " cem: " << cementationexponents[i] << endl; //cout << " sat: " << saturationexponents[i] << endl; - + } if (waterresistivity <= 0) { if (isMaster) cout << "Error: Water resistivity must be positive." << endl; usageandexit(); } - - // 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()); - bool doInterpolate = false; + + // 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()); + bool doInterpolate = false; if (interpolationPoints > 1) { doInterpolate = true; } @@ -656,15 +656,15 @@ try */ MonotCubicInterpolator WaterSaturationVsCapPressure; - + double largestSaturationInterval = Swor-Swir; - + double Ptestvalue = Pcmax; - + while (largestSaturationInterval > (Swor-Swir)/500.0) { // cout << Ptestvalue << endl; if (Pcmax == Pcmin) { - // This is a dummy situation, we go through once and then + // This is a dummy situation, we go through once and then // we are finished (this will be triggered by zero permeability) Ptestvalue = Pcmin; largestSaturationInterval = 0; @@ -687,14 +687,14 @@ try Ptestvalue = SatDiff.first; largestSaturationInterval = SatDiff.second; } - + // Check for saneness of Ptestvalue: if (std::isnan(Ptestvalue) || std::isinf(Ptestvalue)) { if (isMaster) cerr << "ERROR: Ptestvalue was inf or nan" << endl; break; // Jump out of while-loop, just print out the results // up to now and exit the program } - + double waterVolume = 0.0; for (auto c = upscaler.grid().leafbegin<0>(); c != upscaler.grid().leafend<0>(); ++c) { size_t cell_idx = ecl_idx[c->index()]; @@ -706,17 +706,17 @@ try PtestvalueCell = Ptestvalue; double Jvalue = sqrt(permxs[cell_idx] * milliDarcyToSqMetre /poros[cell_idx]) * PtestvalueCell / surfaceTension; //cout << "JvalueCell: " << Jvalue << endl; - waterSaturationCell + waterSaturationCell = InvJfunctions[int(satnums[cell_idx])-1].evaluate(Jvalue); } waterVolume += waterSaturationCell * cellPoreVolumes[cell_idx]; } - + WaterSaturationVsCapPressure.addPair(Ptestvalue, waterVolume/poreVolume); } - + //cout << WaterSaturationVsCapPressure.toString(); - + // Now, it may happen that we have a large number of cells, and // some cells with near zero poro and perm. This may cause that // Pcmax has been estimated so high that it does not affect Sw @@ -733,18 +733,18 @@ try } usageandexit(); } - MonotCubicInterpolator CapPressureVsWaterSaturation(WaterSaturationVsCapPressure.get_fVector(), + MonotCubicInterpolator CapPressureVsWaterSaturation(WaterSaturationVsCapPressure.get_fVector(), WaterSaturationVsCapPressure.get_xVector()); /***************************************************************** * Step 7: - * + * * Loop through a given number of uniformly distributed saturation points * and upscale conductivity 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 + * 1) Loop over all cells to find the saturation value given the * capillary pressure found in (a). Given the saturation value, find the * conductivity in the cell given its physical properties * 2) Upscale conductivity for the geometry. @@ -756,7 +756,7 @@ try Matrix zeroMatrix(3,3,(double*)0); zero(zeroMatrix); - + vector WaterSaturation; // This will hold re-upscaled water saturation for the computed pressure points. vector > UpscaledConductivity; // 'tensorElementCount' phaseperm values per pressurepoint. @@ -766,14 +766,14 @@ try WaterSaturation.push_back(0.0); vector tmp; UpscaledConductivity.push_back(tmp); - for (int voigtIdx=0; voigtIdx < tensorElementCount; ++voigtIdx) { - UpscaledConductivity[idx].push_back(0.0); - } + for (int voigtIdx=0; voigtIdx < tensorElementCount; ++voigtIdx) { + UpscaledConductivity[idx].push_back(0.0); + } } // Make vector of capillary pressure points corresponding to uniformly distributed // saturation points between Swor and Swir. - + vector pressurePoints; for (int pointidx = 1; pointidx <= points; ++pointidx) { // pointidx=1 corresponds to Swir, pointidx=points to Swor. @@ -789,7 +789,7 @@ try for (int idx=0; idx < points; ++idx) { node_vs_pressurepoint.push_back(0); } - + #ifdef HAVE_MPI // Distribute work load over mpi nodes. for (int idx=0; idx < points; ++idx) { @@ -799,27 +799,27 @@ try /*if (isMaster) { cout << "Pressure point " << idx << " assigned to node " << node_vs_pressurepoint[idx] << endl; }*/ - } + } #endif clock_t start_upscale_wallclock = clock(); - + double waterVolumeLF = 0.0; // Now loop through the vector of capillary pressure points that // this node should compute. for (int pointidx = 0; pointidx < points; ++pointidx) { - + double accRes = 0.0; // Should "I" (mpi-wise) compute this pressure point? if (node_vs_pressurepoint[pointidx] == mpi_rank) { - + Ptestvalue = pressurePoints[pointidx]; - + cout << "Upscaling resistivity for Pc = " << Ptestvalue; flush(cout); - + //start = clock(); - + // Loop over each cell again to find saturations given this particular // capillary pressure: waterVolumeLF = 0.0; @@ -833,26 +833,26 @@ try //cout << endl << "Cell no. " << cell_idx << " satnum: " << satnums[cell_idx] << endl; double Jvalue = sqrt(permxs[cell_idx] * milliDarcyToSqMetre / poros[cell_idx]) * Ptestvalue / surfaceTension; //cout << "JvalueCell: " << Jvalue << endl; - double waterSaturationCell + double waterSaturationCell = InvJfunctions[int(satnums[cell_idx])-1].evaluate(Jvalue); //cout << "WatersaturationCell: " << waterSaturationCell << endl; waterVolumeLF += waterSaturationCell * cellPoreVolumes[cell_idx]; - + // Compute cell resistivity. We use a cutoff-value as we // easily divide by zero here. When water saturation is // zero, we get 'inf', which is circumvented by the cutoff value. - + if ((satnums[cell_idx] == mud1rocktype) || (satnums[cell_idx] == mud2rocktype)) { // Handle mud specially resistivityCell = mudresistivity; } else { - resistivityCell - = min(a_lithologycoeff * waterresistivity / pow(poros[cell_idx], cementationexponents[satnums[cell_idx]]) - / pow(waterSaturationCell, saturationexponents[satnums[cell_idx]]), + resistivityCell + = min(a_lithologycoeff * waterresistivity / pow(poros[cell_idx], cementationexponents[satnums[cell_idx]]) + / pow(waterSaturationCell, saturationexponents[satnums[cell_idx]]), resistivityCutoff); } - } + } // Insert conductivity (reciprocal of resistivity) into // the grid for upscaling: Matrix cellCond = zeroMatrix; @@ -865,18 +865,18 @@ try //cout << "cell " << cell_idx << " resistivity" << resistivityCell << endl; } } - - + + // Output average resistity cout << ", Arith. mean res = " << accRes/float(tesselatedCells) << " ohms, "; - + // Call upscaling code (SINTEF/FRAUNHOFER) - + Matrix condTensor; condTensor = upscaler.upscaleSinglePhase(); - - - + + + // Here we recalculate the upscaled water saturation, // although it is already known when we asked for the // pressure point to compute for. Nonetheless, we @@ -885,7 +885,7 @@ try // points are not perfectly uniformly distributed) WaterSaturation[pointidx] = waterVolumeLF/poreVolume; cout << WaterSaturation[pointidx] << endl; - + invert(condTensor); Matrix resTensor(condTensor); @@ -959,10 +959,10 @@ try #ifdef HAVE_MPI // Sum the upscaling time used by all processes double timeused_total; - MPI_Reduce(&timeused_upscale_wallclock, &timeused_total, 1, MPI_DOUBLE, + MPI_Reduce(&timeused_upscale_wallclock, &timeused_total, 1, MPI_DOUBLE, MPI_SUM, 0, MPI_COMM_WORLD); double avg_upscaling_time_pr_point = timeused_total/(double)points; - + #else double avg_upscaling_time_pr_point = timeused_upscale_wallclock / (double)points; #endif @@ -976,10 +976,10 @@ try if (WaterSaturation.size() == 0) { return(1); // non-zero return value, this means something wrong with input data. } - + if (isMaster) { stringstream outputtmp; - + // Print a table of all computed values: outputtmp << "######################################################################" << endl; outputtmp << "# Results from upscaling resistivity."<< endl; @@ -989,10 +989,10 @@ try #endif time_t now = std::time(NULL); outputtmp << "# Finished: " << asctime(localtime(&now)); - + utsname hostname; uname(&hostname); outputtmp << "# Hostname: " << hostname.nodename << endl; - + outputtmp << "#" << endl; outputtmp << "# Eclipse file: " << ECLIPSEFILENAME << endl; outputtmp << "# cells: " << tesselatedCells << endl; @@ -1032,8 +1032,8 @@ try outputtmp << "# minPerm: " << options["minPerm"] << endl; outputtmp << "# minPoro: " << options["minPoro"] << endl; - if (doInterpolate) { - outputtmp << "# interpolate: " << options["interpolate"] << " points" << endl; + if (doInterpolate) { + outputtmp << "# interpolate: " << options["interpolate"] << " points" << endl; } outputtmp << "#" << endl; outputtmp << "# Archie options (global):" << endl; @@ -1054,37 +1054,37 @@ try } } outputtmp << "# " << endl; - if (doInterpolate) { - outputtmp << "# NB: Data points shown are interpolated." << endl; + if (doInterpolate) { + outputtmp << "# NB: Data points shown are interpolated." << endl; } outputtmp << "######################################################################" << endl; if (isFixed) { outputtmp << "# Pc (Pa) Sw Rxx Ryy Rzz" << endl; } else if (isLinear) { - outputtmp << "# Pc (Pa) Sw Rxx Ryy Rzz Ryz Rxz Rxy Rzy Rzx Ryx " << endl; + outputtmp << "# Pc (Pa) Sw Rxx Ryy Rzz Ryz Rxz Rxy Rzy Rzx Ryx " << endl; } else if (isPeriodic) { - outputtmp << "# Pc (Pa) Sw Rxx Ryy Rzz Ryz Rxz Rxy Rzy Rzx Ryx " << endl; + outputtmp << "# Pc (Pa) Sw Rxx Ryy Rzz Ryz Rxz Rxy Rzy Rzx Ryx " << endl; } - + vector Pvalues = pressurePoints; // WaterSaturation.get_xVector(); vector Satvalues = WaterSaturation; //.get_fVector(); - + /* Rearrange the UpscaledConductivity array so that voigtIdx is the first index - (to facilitate interpolation in the then last variable) + (to facilitate interpolation in the then last variable) (results from this is only to be trusted on the master node) */ - + vector > ResDirValues; // voigtIdx is first index. for (int voigtIdx=0; voigtIdx < tensorElementCount; ++voigtIdx) { vector tmp; ResDirValues.push_back(tmp); } - // Loop over all pressure points + // Loop over all pressure points for (int idx=0; idx < points; ++idx) { Matrix condTensor(zeroMatrix); for (int voigtIdx = 0; voigtIdx < tensorElementCount; ++voigtIdx) { @@ -1099,7 +1099,7 @@ try /* for (int voigtIdx = 0; voigtIdx < tensorElementCount; ++voigtIdx) { ResDirValues.push_back(Res[voigtIdx].get_fVector()); }*/ - + /* If user wants interpolated output, do monotone cubic interpolation by modifying the data vectors that are to be printed */ if (doInterpolate) { @@ -1121,7 +1121,7 @@ try } // Now capillary pressure and computed conductivity-values must be viewed as functions // of saturation, and then interpolated on the uniform saturation grid. - + // Now overwrite existing Pvalues and ResDirValues-data with interpolated data: MonotCubicInterpolator PvaluesVsSaturation(Satvalues, Pvalues); Pvalues.clear(); @@ -1138,31 +1138,31 @@ try // Now also overwrite Satvalues Satvalues.clear(); Satvalues = SatvaluesInterp; - + } - + /* Output computed resistivity data */ for (size_t i=0; i < Satvalues.size(); ++i) { // Note: The Interpolator-object's values contain the log10 of the real values. outputtmp << std::showpoint << std::setw(14) << Pvalues[i]; - outputtmp << std::showpoint << std::setw(14) << Satvalues[i]; - + outputtmp << std::showpoint << std::setw(14) << Satvalues[i]; + for (int voigtIdx = 0; voigtIdx < tensorElementCount; ++voigtIdx) { outputtmp << showpoint << setw(14) << ResDirValues[voigtIdx][i]; } outputtmp << endl; } - + cout << outputtmp.str(); - + /* Possibly write to output file */ 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(); + outfile.close(); } } @@ -1172,4 +1172,3 @@ catch (const std::exception &e) { std::cerr << "Program threw an exception: " << e.what() << "\n"; throw; } - diff --git a/examples/upscale_elasticity.cpp b/examples/upscale_elasticity.cpp index a66acb3..0fef77c 100644 --- a/examples/upscale_elasticity.cpp +++ b/examples/upscale_elasticity.cpp @@ -10,7 +10,7 @@ //! //============================================================================== #ifdef HAVE_CONFIG_H -# include "config.h" +# include "config.h" #endif #include @@ -470,7 +470,7 @@ int main(int argc, char** argv) try { try { - if (argc < 2 || strcmp(argv[1],"-h") == 0 + if (argc < 2 || strcmp(argv[1],"-h") == 0 || strcmp(argv[1],"--help") == 0 || strcmp(argv[1],"-?") == 0) { syntax(argv); diff --git a/examples/upscale_perm.cpp b/examples/upscale_perm.cpp index 6cf8989..818804c 100644 --- a/examples/upscale_perm.cpp +++ b/examples/upscale_perm.cpp @@ -22,18 +22,18 @@ /** @file upscale_perm.C * @brief Upscales permeability - * - * Upscales permeability for a given Eclipse-model. + * + * Upscales permeability for a given Eclipse-model. * Upscaling is performed for periodic, linear and * fixed boundary conditions and output is ASCII to standard out. - * + * * Input is Eclipse grid format specifying the corner-point * grid (must be of shoebox-shape, but this condition is slightly relaxed * on top and bottom surfaces). - * + * * The input eclipse file must specify the permeability properties * for each cell. - * + * * If only PERMX is supplied, isotropic permeability is assumed in each cell * * If PERMX, PERMY and PERMZ are supplied, diagonal anisotropic permeability @@ -78,7 +78,7 @@ void usage() { " the terminal (standard out)." << endl << "-bc -- which boundary conditions to compute for. " << endl << " may contain a combination of the" << endl << - " letters lfp to compute linear (l), fixed (f) " << endl << + " letters lfp to compute linear (l), fixed (f) " << endl << " and periodic (p) boundary conditions." << endl << " Default: f (fixed boundary conditions)" << endl << "-minPerm -- Minimum floating point value allowed for" << endl << @@ -88,7 +88,7 @@ void usage() { /** @brief Upscales permeability - + @param varnum Number of input arguments @param vararg Input arguments @return int @@ -99,12 +99,12 @@ int upscale(int varnum, char** vararg) { cout << "Error: No eclipsefile provided" << endl; usage(); exit(1); - } + } map options; - options.insert(make_pair("output", "")); // If this is set, output goes to screen and to this file - options.insert(make_pair("bc", "f")); // Fixed boundary conditions are default + options.insert(make_pair("output", "")); // If this is set, output goes to screen and to this file + options.insert(make_pair("bc", "f")); // Fixed boundary conditions are default options.insert(make_pair("minPerm", "1e-9")); // Minimum allowable permeability value (for diagonal tensor entries) - + options.insert(make_pair("linsolver_tolerance", "1e-8")); // residual tolerance for linear solver options.insert(make_pair("linsolver_verbosity", "0")); // verbosity level for linear solver options.insert(make_pair("linsolver_max_iterations", "0")); // Maximum number of iterations allow, specify 0 for default @@ -128,11 +128,11 @@ int upscale(int varnum, char** vararg) { cout << "Parsed command line option: " << searchfor << " := " << vararg[argidx+1] << endl; eclipseindex += 2; } - } + } const char* ECLIPSEFILENAME(vararg[eclipseindex]); - + // Test if filename exists and is readable ifstream eclipsefile(ECLIPSEFILENAME, ios::in); if (eclipsefile.fail()) { @@ -141,15 +141,15 @@ int upscale(int varnum, char** vararg) { exit(1); } eclipsefile.close(); - - // Check validity of boundary conditions chosen, and make booleans + + // Check validity of boundary conditions chosen, and make booleans // for boundary conditions, this allows more readable code later bool isFixed, isLinear, isPeriodic; isFixed = isLinear = isPeriodic = false; - + // Read in default or user-specified boundary conditions: string boundcond(options["bc"]); - + // Length of string must be between 1 and 3: if (boundcond.length()>= 1 && boundcond.length() <= 3) { if (boundcond.find(string("p")) < 3) { @@ -161,12 +161,12 @@ int upscale(int varnum, char** vararg) { if (boundcond.find(string("l")) < 3) { isLinear = true; } - + // If no boundary conditions are set now, issue error: if (!isFixed && !isLinear && !isPeriodic) { cerr << "Error: No boundary conditions specified: " << boundcond << endl; usage(); - exit(1); + exit(1); } } else { @@ -175,20 +175,20 @@ int upscale(int varnum, char** vararg) { exit(1); } - + // Variables for timing/profiling clock_t start, finish; double timeused = 0; // reusable variable double timeused_periodic_tesselation = 0, timeused_nonperiodic_tesselation = 0; double timeused_periodic = 0, timeused_fixed = 0, timeused_linear = 0; - + cout << endl; - + // Storage for upscaled results: using Opm::SinglePhaseUpscaler; typedef SinglePhaseUpscaler::permtensor_t Matrix; Matrix Kfixed, Klinear, Kperiodic; - + /*********************************************************************** * Load geometry and data from Eclipse file @@ -200,18 +200,18 @@ int upscale(int varnum, char** vararg) { finish = clock(); timeused = (double(finish)-double(start))/CLOCKS_PER_SEC; cout << " (" << timeused <<" secs)" << endl; - + // Check that we have the information we need from the eclipse file, we will check PERM-fields later - if (! (deck.hasKeyword("SPECGRID") && deck.hasKeyword("COORD") && deck.hasKeyword("ZCORN"))) { - cerr << "Error: Did not find SPECGRID, COORD and ZCORN in Eclipse file " << ECLIPSEFILENAME << endl; - usage(); - exit(1); + if (! (deck.hasKeyword("SPECGRID") && deck.hasKeyword("COORD") && deck.hasKeyword("ZCORN"))) { + cerr << "Error: Did not find SPECGRID, COORD and ZCORN in Eclipse file " << ECLIPSEFILENAME << endl; + usage(); + exit(1); } /***************************************************************** - * Tesselate grid - * + * Tesselate grid + * * Possibly twice because, the grid must be massaged slightly * (crop top and bottom) for periodic boundary conditions. These * modifications ruin the computations for linear and fixed @@ -237,9 +237,9 @@ int upscale(int varnum, char** vararg) { if (isFixed || isLinear) { cout << "Tesselating non-periodic grid ..."; start = clock(); - upscaler_nonperiodic.init(deck, + upscaler_nonperiodic.init(deck, isFixed ? SinglePhaseUpscaler::Fixed : SinglePhaseUpscaler::Linear, - minPerm, linsolver_tolerance, linsolver_verbosity, linsolver_type, + minPerm, linsolver_tolerance, linsolver_verbosity, linsolver_type, twodim_hack, linsolver_maxit, linsolver_prolongate_factor, smooth_steps); finish = clock(); timeused_nonperiodic_tesselation = (double(finish)-double(start))/CLOCKS_PER_SEC; @@ -257,8 +257,8 @@ int upscale(int varnum, char** vararg) { } - - + + /********************************************************************* * Do porosity upscaling * @@ -273,11 +273,11 @@ int upscale(int varnum, char** vararg) { upscaledPorosity = upscaler_nonperiodic.upscalePorosity(); } } - + /********************************************************************* - * Do single-phase permeability upscaling + * Do single-phase permeability upscaling */ - + if (isFixed) { cout << "Compute for fixed boundary conditions: ... "; start = clock(); @@ -290,7 +290,7 @@ int upscale(int varnum, char** vararg) { cout << Kfixed << endl; cout << endl; } - + if (isLinear) { cout << "Compute for linear boundary conditions: ... " << endl; start = clock(); @@ -303,26 +303,26 @@ int upscale(int varnum, char** vararg) { cout << " ( " << timeused_linear << " secs)" << endl; cout << endl << endl; } - - if (isPeriodic) { + + if (isPeriodic) { cout << "Compute for periodic boundary conditions: ... "; start = clock(); upscaler_periodic.setBoundaryConditionType(SinglePhaseUpscaler::Periodic); Kperiodic = upscaler_periodic.upscaleSinglePhase(); Kperiodic *= 1.0/(Opm::prefix::milli*Opm::unit::darcy); finish = clock(); - timeused_periodic = (double(finish)-double(start))/CLOCKS_PER_SEC; - cout << " (" << timeused_periodic << " secs)" << endl; + timeused_periodic = (double(finish)-double(start))/CLOCKS_PER_SEC; + cout << " (" << timeused_periodic << " secs)" << endl; cout << Kperiodic << endl; cout << endl; - } - + } + /*********************************************************************** * Output results to stdout or optionally to file */ - + stringstream outputtmp; - + // Print a table of all computed values: outputtmp << "###############################################################################" << endl; outputtmp << "# Results from upscaling permeability."<< endl; @@ -338,7 +338,7 @@ int upscale(int varnum, char** vararg) { outputtmp << "#" << endl; outputtmp << "# Options used:" << endl; outputtmp << "# Boundary conditions: "; - if (isFixed) outputtmp << "Fixed (no-flow) "; + if (isFixed) outputtmp << "Fixed (no-flow) "; if (isPeriodic) outputtmp << "Periodic "; if (isLinear) outputtmp << "Linear "; outputtmp << endl; @@ -346,7 +346,7 @@ int upscale(int varnum, char** vararg) { outputtmp << "#" << endl; outputtmp << "# If both linear and fixed boundary conditions are calculated, " << endl << "# the nonperiodic tesselation is done only once" << endl << "# " << endl << "#" << endl; - + if (isFixed) { outputtmp << "# Upscaled permeability for fixed boundary conditions:" << endl; outputtmp << "# Tesselation time: " << timeused_nonperiodic_tesselation << " s" << endl; @@ -367,9 +367,9 @@ int upscale(int varnum, char** vararg) { outputtmp << "# Tesselation time: " << timeused_periodic_tesselation << " s" << endl; outputtmp << "# Computation time: " << timeused_periodic << " s" << endl; outputtmp << Kperiodic; - } + } cout << endl << outputtmp.str(); - + if (options["output"] != "") { cout << "Writing results to " << options["output"] << endl; ofstream outfile; @@ -377,14 +377,14 @@ int upscale(int varnum, char** vararg) { outfile << outputtmp.str(); outfile.close(); } - return 0; + return 0; } } // namespace anonymous /** @brief Upscales permeability - + @param varnum Number of input arguments @param vararg Input arguments @return int @@ -396,4 +396,3 @@ catch (const std::exception &e) { std::cerr << "Program threw an exception: " << e.what() << "\n"; throw; } - diff --git a/examples/upscale_relpermvisc.cpp b/examples/upscale_relpermvisc.cpp index bec02af..2d4ffe1 100644 --- a/examples/upscale_relpermvisc.cpp +++ b/examples/upscale_relpermvisc.cpp @@ -20,12 +20,12 @@ /** @file upscale_relpermvisc.C @brief Upscales relative permeability as a funtion of watersaturation in the viscous limit. - - Description: - + + Description: + Upscaling of relative permeability (two-phase) using steady-state in the viscous limit, on Eclipse corner-point geometries in shoe-box format - + Reads in a lithofacies geometry in Eclipse format, reads in viscosities for oil and water from command line, relpermcurve(S_w) for each stone type from file, and calculates upscaled (three @@ -34,35 +34,35 @@ Fixed, linear and periodic boundary conditions are supported, yielding 3, 9 or 9 values respectively (but only 6 significant for periodic) for relative permeability at each saturation point. - - The relative permeability computation is based on + + The relative permeability computation is based on - Steady-state, viscous limit. v_w/v_o = constant, === lambda_w/lamda_o = constant - No gravitational effects. - No capillary effects. - + Steps in the code: - + 1: Process command line options. - 2: Read Eclipse file + 2: Read Eclipse file 3: Read relperm-function for each stone-type. 4: Tesselate the grid (Sintef code) 5: Generate simple statistics by looping over the cells. 6: Upscale fractional flow ratio vs. water saturation - 6: Find upper and lower bounds for fractional flow ratio v_w/v_o + 6: Find upper and lower bounds for fractional flow ratio v_w/v_o 7: Upscale single phase permeability (in order to make relative perm later) - 8: For uniformly distributed water saturation values between upscaled Swir and Swor, + 8: For uniformly distributed water saturation values between upscaled Swir and Swor, a: Find fractional flow ratio corresponding to wanted upscaled water saturation b: Model water saturation in each cell given fractional flow ratio c: Model phase permeability in each cell given cell water saturation and inputted relative permeability curves. d: Upscale phase permeability - 9: Repeat step 8 for the oil phase. - 10: Print output to screen and optionally to files. - - The relperm-functions must be defined in text files, with water saturation in - the first column, and then relative permeability for water in column 2, and relative + 9: Repeat step 8 for the oil phase. + 10: Print output to screen and optionally to files. + + The relperm-functions must be defined in text files, with water saturation in + the first column, and then relative permeability for water in column 2, and relative permeability for oil in column 3. Lines starting with -- or # are ignored. - + Units for (dynamic) viscosity are assumed to be in Pascal seconds. (Pa s) 1000 Pa s = 1 centiPoise. */ @@ -99,18 +99,18 @@ using namespace std; namespace { /** - The usage() function displays a message with syntax and available - options for the program, and is typically called whenever the code - encounters any errors in the options or in the command line syntax - (missing filenames etc.) + The usage() function displays a message with syntax and available + options for the program, and is typically called whenever the code + encounters any errors in the options or in the command line syntax + (missing filenames etc.) */ void usage() { cout << "Usage: upscale_relpermvisc rock1.txt rock2.txt... (isotropic case)" << endl << " upscale_relpermvisc rock1_water.txt rock1_oil.txt rock2_water.txt rock2_oil.txt ... (anisotropic case)" << endl << " where the options are:" << endl << - "-bc -- which boundary conditions to use." << endl << - " Possible values are f (fixed), " << endl << + "-bc -- which boundary conditions to use." << endl << + " Possible values are f (fixed), " << endl << " l (linear) and p (periodic). Default f." << endl << "-points -- Number of saturation points to upscale for." << endl << " Uniformly distributed within saturation endpoints." << endl << @@ -121,14 +121,14 @@ void usage() "-waterCurve -- the column number in the stone files that represent" << endl << " relative permeability for water. Default 2" << endl << "-oilCurve -- the column number in the stone files that represents " << endl << - " relative permeability for oil. Default 3." << endl << + " relative permeability for oil. Default 3." << endl << "-outputWater -- filename for where to write upscaled values for" << endl << " water relperm. If not supplied, output will only" << endl << " go to the terminal (standard out)." << endl << "-outputOil -- ditto" << endl << - "-interpolate -- If supplied and > 1, the output data points will be" << endl << - " interpolated using monotone cubic interpolation" << endl << - " on a uniform grid with the specified number of" << endl << + "-interpolate -- If supplied and > 1, 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 << "" << endl << "-minPerm -- Minimum floating point value allowed for" << endl << " phase permeability in computations. If set to zero," << endl << @@ -142,7 +142,7 @@ void usage() // -maxPermContrast // -minPoro // -saturationThreshold - + } void usageandexit() { @@ -154,17 +154,17 @@ void usageandexit() { int main(int varnum, char** vararg) try -{ +{ // Variables used for timing/profiling: clock_t start, finish; - double timeused = 0, timeused_tesselation = 0; //, timeused_upscale_acc_water = 0, timeused_upscale_acc_oil = 0; + double timeused = 0, timeused_tesselation = 0; //, timeused_upscale_acc_water = 0, timeused_upscale_acc_oil = 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 @@ -175,7 +175,7 @@ try usage(); exit(1); } - + /* Populate options-map with default values */ @@ -207,15 +207,15 @@ try if (isMaster) cerr << "Error: No Eclipsefile or stonefiles found on command line." << endl; usageandexit(); } - - /* Loop over all command line options in order to look - for options. - + + /* Loop over all command line options in order to look + for options. + argidx loops over all the arguments here, and updates the variable 'argeclindex' *if* it finds any legal options, 'argeclindex' is so that vararg[argeclindex] = the eclipse - filename. If options are illegal, argeclindex will be wrong, - + filename. If options are illegal, argeclindex will be wrong, + */ int argeclindex = 0; for (int argidx = 1; argidx < varnum; argidx += 2) { @@ -232,135 +232,135 @@ try usageandexit(); } } - else { - // if vararg[argidx] does not start in '-', + 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, + break; // out of for-loop, } } - + // 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. - + // Boolean set to true of input permeability in eclipse-file has diagonal anisotropy. // (full-tensor anisotropy will be ignored) bool anisotropic_input = false; - + // argcindex 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) { if (isMaster) cerr << "Error: No J-functions found on command line." << endl; usageandexit(); } - - /* Check validity of boundary conditions chosen, and make booleans + + /* Check validity of boundary conditions chosen, and make booleans for boundary conditions, this allows more readable code later. */ - bool isFixed, isLinear, isPeriodic; - SinglePhaseUpscaler::BoundaryConditionType boundaryCondition; - int tensorElementCount; // Number of independent elements in resulting tensor. - if (options["bc"].substr(0,1) == "f") { - isFixed = true; isLinear = false; isPeriodic = false; - boundaryCondition = SinglePhaseUpscaler::Fixed; - tensorElementCount = 3; // Diagonal - } - else if (options["bc"].substr(0,1) == "l") { - isLinear = true; isFixed = false; isPeriodic = false; - boundaryCondition = SinglePhaseUpscaler::Linear; - tensorElementCount = 9; // Full-tensor - } - else if (options["bc"].substr(0,1) == "p") { - isPeriodic = true; isLinear = false; isFixed = false; - boundaryCondition = SinglePhaseUpscaler::Periodic; - tensorElementCount = 9; // Symmetric. - } - else { - if (isMaster) cout << "Invalid boundary condition. Only one of the letters f, l or p are allowed." << endl; + bool isFixed, isLinear, isPeriodic; + SinglePhaseUpscaler::BoundaryConditionType boundaryCondition; + int tensorElementCount; // Number of independent elements in resulting tensor. + if (options["bc"].substr(0,1) == "f") { + isFixed = true; isLinear = false; isPeriodic = false; + boundaryCondition = SinglePhaseUpscaler::Fixed; + tensorElementCount = 3; // Diagonal + } + else if (options["bc"].substr(0,1) == "l") { + isLinear = true; isFixed = false; isPeriodic = false; + boundaryCondition = SinglePhaseUpscaler::Linear; + tensorElementCount = 9; // Full-tensor + } + else if (options["bc"].substr(0,1) == "p") { + isPeriodic = true; isLinear = false; isFixed = false; + boundaryCondition = SinglePhaseUpscaler::Periodic; + tensorElementCount = 9; // Symmetric. + } + else { + if (isMaster) cout << "Invalid boundary condition. Only one of the letters f, l or p are allowed." << endl; usageandexit(); - } - + } + // Index values for the phases. We are going to use the same code // twice to avoid source code duplication in this file, first for // water, and then for oil. To accomplish this, we have some data // in a vector where the first element is the waterphase and the // second is the oilphase, defined by these constants. - const int waterPhaseIndex = 0; - const int oilPhaseIndex = 1; - + const int waterPhaseIndex = 0; + const int oilPhaseIndex = 1; + const double points = atof(options["points"].c_str()); - + const int waterCurveColumn = atoi(options["waterCurve"].c_str()); const int oilCurveColumn = atoi(options["oilCurve"].c_str()); //const int outputprecision = atoi(options["outputprecision"].c_str()); - - - // Vector of column index for input rel perm curves - vector relPermCurves; - relPermCurves.push_back(waterCurveColumn); - relPermCurves.push_back(oilCurveColumn); - // Now this vector can be accessed with - // relPermCurves[waterPhaseIndex] etc. - - // Vector of viscosities - vector viscosities; - viscosities.push_back(atof(options["waterViscosity"].c_str())); - viscosities.push_back(atof(options["oilViscosity"].c_str())); - // Now this vector can be accessed with - // viscosities[waterPhaseIndex] to give the viscosity for water. - - - // 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()); - bool doInterpolate = false; + + + // Vector of column index for input rel perm curves + vector relPermCurves; + relPermCurves.push_back(waterCurveColumn); + relPermCurves.push_back(oilCurveColumn); + // Now this vector can be accessed with + // relPermCurves[waterPhaseIndex] etc. + + // Vector of viscosities + vector viscosities; + viscosities.push_back(atof(options["waterViscosity"].c_str())); + viscosities.push_back(atof(options["oilViscosity"].c_str())); + // Now this vector can be accessed with + // viscosities[waterPhaseIndex] to give the viscosity for water. + + + // 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()); + bool doInterpolate = false; if (interpolationPoints > 1) { doInterpolate = true; } - + /*********************************************************************** * Step 2: * Load geometry and data from Eclipse file */ - // Read data from the Eclipse file and + // Read data from the Eclipse file and // populate our vectors with data from the file - + // Test if filename exists and is readable ifstream eclipsefile(ECLIPSEFILENAME, ios::in); if (eclipsefile.fail()) { if (isMaster) cerr << "Error: Filename " << ECLIPSEFILENAME << " not found or not readable." << endl; usageandexit(); } - eclipsefile.close(); - + eclipsefile.close(); + if (isMaster) cout << "Parsing Eclipse file <" << ECLIPSEFILENAME << "> ... "; flush(cout); start = clock(); - + auto deck = RelPermUpscaleHelper::parseEclipseFile(ECLIPSEFILENAME); finish = clock(); timeused = (double(finish)-double(start))/CLOCKS_PER_SEC; if (isMaster) cout << " (" << timeused <<" secs)" << endl; - - // Check that we have the information we need from the eclipse file: - if (! (deck.hasKeyword("SPECGRID") && deck.hasKeyword("COORD") && deck.hasKeyword("ZCORN") - && deck.hasKeyword("PORO") && deck.hasKeyword("PERMX") && deck.hasKeyword("SATNUM"))) { - cerr << "Error: Did not find SPECGRID, COORD and ZCORN in Eclipse file " << ECLIPSEFILENAME << endl; - usage(); - exit(1); - } - + + // Check that we have the information we need from the eclipse file: + if (! (deck.hasKeyword("SPECGRID") && deck.hasKeyword("COORD") && deck.hasKeyword("ZCORN") + && deck.hasKeyword("PORO") && deck.hasKeyword("PERMX") && deck.hasKeyword("SATNUM"))) { + cerr << "Error: Did not find SPECGRID, COORD and ZCORN in Eclipse file " << ECLIPSEFILENAME << endl; + usage(); + exit(1); + } + vector poros = deck["PORO"].back().getRawDoubleData(); vector permxs = deck["PERMX"].back().getRawDoubleData(); - + // Load anisotropic (only diagonal supported) input if present in grid vector permys, permzs; - + if (deck.hasKeyword("PERMY") && deck.hasKeyword("PERMZ")) { anisotropic_input = true; permys = deck["PERMY"].back().getRawDoubleData(); @@ -369,35 +369,35 @@ try cout << "Info: PERMY and PERMZ present, going into anisotropic input mode, no J-functions\n" << " Options -relPermCurve and -jFunctionCurve is meaningless.\n"; } - } - - - /* Initialize a default satnums-vector with only "ones" (meaning only one rocktype) */ - vector satnums(poros.size(), 1); - + } + + + /* Initialize a default satnums-vector with only "ones" (meaning only one rocktype) */ + vector satnums(poros.size(), 1); + if (deck.hasKeyword("SATNUM")) { satnums = deck["SATNUM"].back().getIntData(); - } + } else if (deck.hasKeyword("ROCKTYPE")) { satnums = deck["ROCKTYPE"].back().getIntData(); - } - else { - if (isMaster) cout << "Warning: SATNUM or ROCKTYPE not found in input file, assuming only one rocktype" << endl; - } - + } + else { + if (isMaster) cout << "Warning: SATNUM or ROCKTYPE not found in input file, assuming only one rocktype" << endl; + } + /* Sanity check/fix on input for each cell: - Check that SATNUM are set sensibly, that is => 0 and < 1000, error if not. - Check that porosity is between 0 and 1, error if not. Set to minPoro if zero or less than minPoro (due to pcmin/max computation) - - Check that permeability is zero or positive. Error if negative. + - Check that permeability is zero or positive. Error if negative. Set to minPerm if zero or less than minPerm. - Check maximum number of SATNUM values (can be number of rock types present) */ const double maxPermContrast = atof(options["maxPermContrast"].c_str()); const double minPerm = atof(options["minPerm"].c_str()); const double maxPerm = atof(options["maxPerm"].c_str()); - const double minPoro = atof(options["minPoro"].c_str()); + const double minPoro = atof(options["minPoro"].c_str()); const double saturationThreshold = atof(options["saturationThreshold"].c_str()); double maxPermInInputFile = 0.0; int cells_truncated_from_below_poro = 0; @@ -410,10 +410,10 @@ try if (isMaster) cout << "Illegal maxPermContrast value" << endl; usageandexit(); } - + int maxSatnum = 0; for (unsigned int i = 0; i < satnums.size(); ++i) { - if (satnums[i] < 0 || satnums[i] > 1000) { + if (satnums[i] < 0 || satnums[i] > 1000) { if (isMaster) cerr << "satnums[" << i << "] = " << satnums[i] << ", not sane, quitting." << endl; usageandexit(); } @@ -459,14 +459,14 @@ try // Explicitly handle "no rock" cells, set them to minimum perm and zero porosity. if (satnums[i] == 0) { permxs[i] = minPerm; - poros[i] = 0; // zero poro is fine for these cells, as they are not + poros[i] = 0; // zero poro is fine for these cells, as they are not // used in pcmin/max computation. if (anisotropic_input) { permys[i] = minPerm; permzs[i] = minPerm; } } - } + } if (cells_truncated_from_below_poro > 0) { cout << "Cells with truncated porosity: " << cells_truncated_from_below_poro << endl; } @@ -480,15 +480,15 @@ try /*************************************************************************** * Step 3: * Load relperm-curves for the stone types. - * We read columns from text-files, syntax allowed is determined - * by MonotCubicInterpolator which actually opens and parses the + * We read columns from text-files, syntax allowed is determined + * by MonotCubicInterpolator which actually opens and parses the * text files. - * + * * How do we determine which mode of operation? * - If PERMY and PERMZ are present in grdecl-file, we are in the anisotropic mode - * + * */ - + int stone_types = int(*(max_element(satnums.begin(), satnums.end()))); std::vector Krw; // Holds relperm-water-curves for each stone type std::vector Kro; // Holds relperm-oil-curves for each stone type @@ -496,12 +496,12 @@ try std::vector Krwx, Krwy, Krwz, Krox, Kroy, Kroz; std::vector rockTypeNames; // Placeholder for the names of the loaded rock perm files, // to be used in final output. - + // Handle two command line input formats, either stone.txt for all stone types // *or* one each. If there is only one stone type, both code blocks below are equivalent. if (! anisotropic_input) { if (varnum == rockfileindex + stone_types) { // one .txt for each stone type in isotropic case - for (int i=0 ; i < stone_types; ++i) { + for (int i=0 ; i < stone_types; ++i) { const char* ROCKFILENAME = vararg[rockfileindex+i]; // Check if rock file exists and is readable: ifstream rockfile(ROCKFILENAME, ios::in); @@ -509,10 +509,10 @@ try if (isMaster) cerr << "Error: Filename " << ROCKFILENAME << " not found or not readable." << endl; usageandexit(); } - rockfile.close(); + rockfile.close(); MonotCubicInterpolator Krwtmp; try { - Krwtmp = MonotCubicInterpolator(ROCKFILENAME, 1, relPermCurves[waterPhaseIndex]); + Krwtmp = MonotCubicInterpolator(ROCKFILENAME, 1, relPermCurves[waterPhaseIndex]); } catch (const char * errormessage) { if (isMaster) { @@ -530,7 +530,7 @@ try } MonotCubicInterpolator Krotmp; try { - Krotmp = MonotCubicInterpolator(ROCKFILENAME, 1, relPermCurves[oilPhaseIndex]); + Krotmp = MonotCubicInterpolator(ROCKFILENAME, 1, relPermCurves[oilPhaseIndex]); } catch (const char * errormessage) { if (isMaster) { @@ -539,7 +539,7 @@ try } usageandexit(); } - + if (!Krotmp.isStrictlyMonotone()) { if (isMaster) { cerr << "Error: Data in column " << relPermCurves[oilPhaseIndex] << " of file " << ROCKFILENAME << endl << @@ -547,26 +547,26 @@ try } usageandexit(); } - + if ((Krwtmp.isStrictlyIncreasing() && Krotmp.isStrictlyIncreasing()) || (Krwtmp.isStrictlyDecreasing() && Krotmp.isStrictlyDecreasing())) { if (isMaster) { cerr << "Error: Input relperm curves are both increasing or decreasing in file " << ROCKFILENAME << endl; } usageandexit(); } - - Krw.push_back(Krwtmp); - Kro.push_back(Krotmp); - + + Krw.push_back(Krwtmp); + Kro.push_back(Krotmp); + if (isMaster) { cout << "Loaded rock file: " << ROCKFILENAME << ", for stone type " << i+1 << endl; } - - rockTypeNames.push_back(ROCKFILENAME); + + rockTypeNames.push_back(ROCKFILENAME); } - } - else if (varnum == rockfileindex + 1) { // one .txt for all stone types + } + else if (varnum == rockfileindex + 1) { // one .txt for all stone types const char* ROCKFILENAME = vararg[rockfileindex]; // Check if rock files exists and is readable: ifstream rockfile(ROCKFILENAME, ios::in); @@ -579,7 +579,7 @@ try rockfile.close(); MonotCubicInterpolator Krwtmp; try { - Krwtmp = MonotCubicInterpolator(ROCKFILENAME, 1, relPermCurves[waterPhaseIndex]); + Krwtmp = MonotCubicInterpolator(ROCKFILENAME, 1, relPermCurves[waterPhaseIndex]); } catch (const char * errormessage) { if (isMaster) { @@ -597,7 +597,7 @@ try } MonotCubicInterpolator Krotmp; try { - Krotmp = MonotCubicInterpolator(ROCKFILENAME, 1, relPermCurves[oilPhaseIndex]); + Krotmp = MonotCubicInterpolator(ROCKFILENAME, 1, relPermCurves[oilPhaseIndex]); } catch (const char * errormessage) { if (isMaster) { @@ -620,22 +620,22 @@ try usageandexit(); } for (int i=0 ; i < stone_types; ++i) { //Insert the same input curves for all rock types - Krw.push_back(Krwtmp); - Kro.push_back(Krotmp); - rockTypeNames.push_back(ROCKFILENAME); + Krw.push_back(Krwtmp); + Kro.push_back(Krotmp); + rockTypeNames.push_back(ROCKFILENAME); } - + if (isMaster) { cout << "Loaded rock file: " << ROCKFILENAME << ", for all stone types" << endl; } } - else { - cerr << "Error: Wrong number of stone-functions provided. " << endl - << "Note that all input arguments after eclipse file are " << endl - << "interpreted as stone functions." << endl; - return 1; - } + else { + cerr << "Error: Wrong number of stone-functions provided. " << endl + << "Note that all input arguments after eclipse file are " << endl + << "interpreted as stone functions." << endl; + return 1; + } } else { // Anisotropic case here! (double set of input curves needed cout << "rock types: " << stone_types << endl; @@ -643,7 +643,7 @@ try cout << "rockfileindex: " << rockfileindex << endl; if (varnum == rockfileindex + 2*stone_types) { // two .txt for each stone type in anisotropic case int rockidx=0; - for (int i=0 ; i < 2*stone_types; i+=2) { + for (int i=0 ; i < 2*stone_types; i+=2) { rockidx++; const char* WATERFILENAME = vararg[rockfileindex+i]; const char* OILFILENAME = vararg[rockfileindex+i+1]; @@ -656,20 +656,20 @@ try } usageandexit(); } - waterfile.close(); + waterfile.close(); if (oilfile.fail()) { if (isMaster) { cerr << "Error: Filename " << OILFILENAME << " not found or not readable." << endl; } usageandexit(); } - oilfile.close(); - + oilfile.close(); + MonotCubicInterpolator Krwxtmp, Krwytmp, Krwztmp, Kroxtmp, Kroytmp, Kroztmp; try { - Krwxtmp = MonotCubicInterpolator(WATERFILENAME, 2, 3); - Krwytmp = MonotCubicInterpolator(WATERFILENAME, 2, 4); - Krwztmp = MonotCubicInterpolator(WATERFILENAME, 2, 5); + Krwxtmp = MonotCubicInterpolator(WATERFILENAME, 2, 3); + Krwytmp = MonotCubicInterpolator(WATERFILENAME, 2, 4); + Krwztmp = MonotCubicInterpolator(WATERFILENAME, 2, 5); } catch (const char * errormessage) { if (isMaster) { @@ -700,9 +700,9 @@ try usageandexit(); } try { - Kroxtmp = MonotCubicInterpolator(OILFILENAME, 2, 3); - Kroytmp = MonotCubicInterpolator(OILFILENAME, 2, 4); - Kroztmp = MonotCubicInterpolator(OILFILENAME, 2, 5); + Kroxtmp = MonotCubicInterpolator(OILFILENAME, 2, 3); + Kroytmp = MonotCubicInterpolator(OILFILENAME, 2, 4); + Kroztmp = MonotCubicInterpolator(OILFILENAME, 2, 5); } catch (const char * errormessage) { if (isMaster) { @@ -732,17 +732,17 @@ try } usageandexit(); } - Krwx.push_back(Krwxtmp); Krwy.push_back(Krwytmp); Krwz.push_back(Krwztmp); - Krox.push_back(Kroxtmp); Kroy.push_back(Kroytmp); Kroz.push_back(Kroztmp); + Krwx.push_back(Krwxtmp); Krwy.push_back(Krwytmp); Krwz.push_back(Krwztmp); + Krox.push_back(Kroxtmp); Kroy.push_back(Kroytmp); Kroz.push_back(Kroztmp); if (isMaster) { cout << "Loaded rock files: " << WATERFILENAME << " and " << OILFILENAME << ", for stone type " << rockidx << endl; } - rockTypeNames.push_back(WATERFILENAME); - rockTypeNames.push_back(OILFILENAME); + rockTypeNames.push_back(WATERFILENAME); + rockTypeNames.push_back(OILFILENAME); } } - else if (varnum == rockfileindex + 2) { // one waterfile and one oilfile for all stone types + else if (varnum == rockfileindex + 2) { // one waterfile and one oilfile for all stone types const char* WATERFILENAME = vararg[rockfileindex]; const char* OILFILENAME = vararg[rockfileindex+1]; // Check if rock files exists and is readable: @@ -754,19 +754,19 @@ try } usageandexit(); } - waterfile.close(); + waterfile.close(); if (oilfile.fail()) { if (isMaster) { cerr << "Error: Filename " << OILFILENAME << " not found or not readable." << endl; } usageandexit(); } - oilfile.close(); + oilfile.close(); MonotCubicInterpolator Krwxtmp, Krwytmp, Krwztmp; try { - Krwxtmp = MonotCubicInterpolator(WATERFILENAME, 2, 3); - Krwytmp = MonotCubicInterpolator(WATERFILENAME, 2, 4); - Krwztmp = MonotCubicInterpolator(WATERFILENAME, 2, 5); + Krwxtmp = MonotCubicInterpolator(WATERFILENAME, 2, 3); + Krwytmp = MonotCubicInterpolator(WATERFILENAME, 2, 4); + Krwztmp = MonotCubicInterpolator(WATERFILENAME, 2, 5); } catch (const char * errormessage) { if (isMaster) { @@ -798,9 +798,9 @@ try } MonotCubicInterpolator Kroxtmp, Kroytmp, Kroztmp; try { - Kroxtmp = MonotCubicInterpolator(OILFILENAME, 2, 3); - Kroytmp = MonotCubicInterpolator(OILFILENAME, 2, 4); - Kroztmp = MonotCubicInterpolator(OILFILENAME, 2, 5); + Kroxtmp = MonotCubicInterpolator(OILFILENAME, 2, 3); + Kroytmp = MonotCubicInterpolator(OILFILENAME, 2, 4); + Kroztmp = MonotCubicInterpolator(OILFILENAME, 2, 5); } catch (const char * errormessage) { if (isMaster) { @@ -830,30 +830,30 @@ try } usageandexit(); } - for (int i=0 ; i < 2*stone_types; i+=2) { - Krwx.push_back(Krwxtmp); Krwy.push_back(Krwytmp); Krwz.push_back(Krwztmp); - Krox.push_back(Kroxtmp); Kroy.push_back(Kroytmp); Kroz.push_back(Kroztmp); - rockTypeNames.push_back(WATERFILENAME); - rockTypeNames.push_back(OILFILENAME); + for (int i=0 ; i < 2*stone_types; i+=2) { + Krwx.push_back(Krwxtmp); Krwy.push_back(Krwytmp); Krwz.push_back(Krwztmp); + Krox.push_back(Kroxtmp); Kroy.push_back(Kroytmp); Kroz.push_back(Kroztmp); + rockTypeNames.push_back(WATERFILENAME); + rockTypeNames.push_back(OILFILENAME); } if (isMaster) { cout << "Loaded rock files: " << WATERFILENAME << " and " << OILFILENAME << ", for all stone types" << endl; } } - else { - cerr << "Error: Wrong number of stone-functions provided. " << endl - << "Note that all input arguments after eclipse file are " << endl - << "interpreted as input functions." << endl; - return 1; - } + else { + cerr << "Error: Wrong number of stone-functions provided. " << endl + << "Note that all input arguments after eclipse file are " << endl + << "interpreted as input functions." << endl; + return 1; + } } /***************************************************************************** * Step 4: * Generate tesselated grid: - * This is a step needed for the later discretization code to figure out which + * 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. */ if (isMaster) cout << "Tesselating grid... "; @@ -866,18 +866,18 @@ try upscaler.init(deck, boundaryCondition, Opm::unit::convert::from(minPerm, Opm::prefix::milli*Opm::unit::darcy), linsolver_tolerance, linsolver_verbosity, linsolver_type, twodim_hack); - + finish = clock(); timeused_tesselation = (double(finish)-double(start))/CLOCKS_PER_SEC; if (isMaster) cout << " (" << timeused_tesselation <<" secs)" << endl; - + /****************************************************************************** * Step 5: * Loop over cells to calculate * - total volume for each rock type * - Upscaled Swir * - Upscaled Swor - * - maximum single phase perm, to be used for permeability contrast control - * + * - maximum single phase perm, to be used for permeability contrast control + * * Find total volume for each stone type * This is information to be used later on as it allows for faster computation * (because in the viscous limit, only the rock type matters for each cells @@ -890,25 +890,25 @@ try vector rocktypeVolume; // rocktypeVolume[rockIdx] == aggregated pore-volume for that type rocktypeVolume.resize(stone_types); - + int tesselatedCells = 0; - + double maxSinglePhasePerm = 0; double Swirvolume = 0; double Sworvolume = 0; - - // Loop through all cells. Add the cells porevolume to the corresponding rock type volume - // Also determine bounds for fractional flow ratio + + // Loop through all cells. Add the cells porevolume to the corresponding rock type volume + // Also determine bounds for fractional flow ratio const std::vector& ecl_idx = upscaler.grid().globalCell(); Dune::CpGrid::Codim<0>::LeafIterator c = upscaler.grid().leafbegin<0>(); for (; c != upscaler.grid().leafend<0>(); ++c) { unsigned int cell_idx = ecl_idx[c->index()]; if (satnums[cell_idx] > 0) { // Satnum zero is "no rock" cellVolumes[cell_idx] = c->geometry().volume(); - cellPoreVolumes[cell_idx] = cellVolumes[cell_idx] * poros[cell_idx]; - - rocktypeVolume[int(satnums[cell_idx])-1] += cellPoreVolumes[cell_idx]; - + cellPoreVolumes[cell_idx] = cellVolumes[cell_idx] * poros[cell_idx]; + + rocktypeVolume[int(satnums[cell_idx])-1] += cellPoreVolumes[cell_idx]; + // Also find max single-phase perm in input file: maxSinglePhasePerm = max( maxSinglePhasePerm, permxs[cell_idx]); double minSw, maxSw; @@ -920,28 +920,28 @@ try } else { minSw = Krwx[int(satnums[cell_idx])-1].getMinimumX().first; - maxSw = Krwx[int(satnums[cell_idx])-1].getMaximumX().first; + maxSw = Krwx[int(satnums[cell_idx])-1].getMaximumX().first; } // Add irreducible water saturation volume Swirvolume += minSw * cellPoreVolumes[cell_idx]; Sworvolume += maxSw * cellPoreVolumes[cell_idx]; } ++tesselatedCells; // keep count (also counts non-rock-cells) - } - + } + double minSinglePhasePerm = max(maxSinglePhasePerm/maxPermContrast, minPerm); - - // Total porevolume and total volume -> upscaled porosity: - double poreVolume = accumulate(cellPoreVolumes.begin(), - cellPoreVolumes.end(), - 0.0); - double volume = accumulate(cellVolumes.begin(), - cellVolumes.end(), - 0.0); - + + // Total porevolume and total volume -> upscaled porosity: + double poreVolume = accumulate(cellPoreVolumes.begin(), + cellPoreVolumes.end(), + 0.0); + double volume = accumulate(cellVolumes.begin(), + cellVolumes.end(), + 0.0); + double Swir = Swirvolume/poreVolume; double Swor = Sworvolume/poreVolume; - + if (isMaster) { cout << "LF Pore volume: " << poreVolume << endl; cout << "LF Volume: " << volume << endl; @@ -950,11 +950,11 @@ try cout << "Upscaled Swmax: " << Swor << endl; //Swor=1-Swmax cout << "Saturation points to be computed: " << points << endl; } - - // Sometimes, if Swmax=1 or Swir=0 in the input tables, the upscaled + + // Sometimes, if Swmax=1 or Swir=0 in the input tables, the upscaled // values can be a little bit larger (within machine precision) and - // the check below fails. Hence, check if these values are within the + // the check below fails. Hence, check if these values are within the // the [0 1] interval within some precision (use linsolver_precision) if (Swor > 1.0 && Swor - linsolver_tolerance < 1.0) { Swor = 1.0; @@ -967,10 +967,10 @@ try cerr << "ERROR: Swir/Swor unsensible. Check your input. Exiting"; } usageandexit(); - } + } - /*********************************************************************** - * Step 6 + /*********************************************************************** + * Step 6 * Upscale fractional flow ratio vs water saturation * * This is upscaled in advance in order to be able to have uniformly distributed @@ -981,43 +981,43 @@ try * is 1/500 of the saturation interval. Monotone cubic * interpolation will be used afterwards for accessing the * tabulated values. - */ - + */ + MonotCubicInterpolator WaterSaturationVsFractionalFlow; double largestSaturationInterval = Swor - Swir; - - // Make fractional flow rato vs watersaturation curve for each stone type - const int interpPoints = 1000; - const double lowerRelPermLimit = 1e-12; - vector FracFlowRatioInv; // We only need the inverse of this + + // Make fractional flow rato vs watersaturation curve for each stone type + const int interpPoints = 1000; + const double lowerRelPermLimit = 1e-12; + vector FracFlowRatioInv; // We only need the inverse of this vector FracFlowRatioInvY, FracFlowRatioInvZ; - for (int stone_idx = 0; stone_idx < stone_types; ++stone_idx) { - MonotCubicInterpolator waterCurve, oilCurve, + for (int stone_idx = 0; stone_idx < stone_types; ++stone_idx) { + MonotCubicInterpolator waterCurve, oilCurve, waterCurveY, oilCurveY, waterCurveYInv, oilCurveYInv, - waterCurveZ, oilCurveZ, waterCurveZInv, oilCurveZInv; + waterCurveZ, oilCurveZ, waterCurveZInv, oilCurveZInv; // NOTE: This is a huge hack for anisotropic_input, we only look at // x-directions, and hope this is good enough for establishing - // which fracflows to choose. The only side effect is probably slightly + // which fracflows to choose. The only side effect is probably slightly // variations in distance between each saturation point. if (! anisotropic_input) { waterCurve = MonotCubicInterpolator(Krw[stone_idx]); - oilCurve = MonotCubicInterpolator(Kro[stone_idx]); + oilCurve = MonotCubicInterpolator(Kro[stone_idx]); } else { waterCurve = MonotCubicInterpolator(Krwx[stone_idx]); - oilCurve = MonotCubicInterpolator(Krox[stone_idx]); + oilCurve = MonotCubicInterpolator(Krox[stone_idx]); waterCurveY = MonotCubicInterpolator(Krwy[stone_idx]); oilCurveY = MonotCubicInterpolator(Kroy[stone_idx]); waterCurveZ = MonotCubicInterpolator(Krwz[stone_idx]); - oilCurveZ = MonotCubicInterpolator(Kroz[stone_idx]); + oilCurveZ = MonotCubicInterpolator(Kroz[stone_idx]); } - MonotCubicInterpolator waterCurveInv(waterCurve.get_fVector(), waterCurve.get_xVector()); - MonotCubicInterpolator oilCurveInv(oilCurve.get_fVector(), oilCurve.get_xVector()); - // Finding min and max bounds for watersaturation - double waterSwMin = waterCurve.getMinimumX().first; - double waterSwMax = waterCurve.getMaximumX().first; - double oilSwMin = oilCurve.getMinimumX().first; - double oilSwMax = oilCurve.getMaximumX().first; + MonotCubicInterpolator waterCurveInv(waterCurve.get_fVector(), waterCurve.get_xVector()); + MonotCubicInterpolator oilCurveInv(oilCurve.get_fVector(), oilCurve.get_xVector()); + // Finding min and max bounds for watersaturation + double waterSwMin = waterCurve.getMinimumX().first; + double waterSwMax = waterCurve.getMaximumX().first; + double oilSwMin = oilCurve.getMinimumX().first; + double oilSwMax = oilCurve.getMaximumX().first; if (anisotropic_input) { waterSwMin = min(waterSwMin, waterCurveY.getMinimumX().first); waterSwMax = max(waterSwMax, waterCurveY.getMaximumX().first); @@ -1027,77 +1027,77 @@ try oilSwMax = max(oilSwMax, oilCurveY.getMaximumX().first); oilSwMin = min(oilSwMin, oilCurveZ.getMinimumX().first); oilSwMax = max(oilSwMax, oilCurveZ.getMaximumX().first); - } + } // Next two lines limit the Sw-span. CHECK!! - double leftBound = min(oilSwMin, waterSwMin); // Maximum of the two curves left bounds - double rightBound = max(oilSwMax, waterSwMax); // Minimum of the two curves right bounds - // Make sure the rel perm values are not too small (avoid infinity in the fracflow) - if (waterCurve.evaluate(leftBound) < lowerRelPermLimit) { - leftBound = waterCurveInv.evaluate(lowerRelPermLimit); - } - if (oilCurve.evaluate(rightBound) < lowerRelPermLimit) { - rightBound = oilCurveInv.evaluate(lowerRelPermLimit); + double leftBound = min(oilSwMin, waterSwMin); // Maximum of the two curves left bounds + double rightBound = max(oilSwMax, waterSwMax); // Minimum of the two curves right bounds + // Make sure the rel perm values are not too small (avoid infinity in the fracflow) + if (waterCurve.evaluate(leftBound) < lowerRelPermLimit) { + leftBound = waterCurveInv.evaluate(lowerRelPermLimit); + } + if (oilCurve.evaluate(rightBound) < lowerRelPermLimit) { + rightBound = oilCurveInv.evaluate(lowerRelPermLimit); // Assume that also holds for aniso input - } - + } + // Uniformly distributed watersaturation points for each rocktype - vector Sw(interpPoints, leftBound); - vector f(interpPoints, 0); - vector fY(interpPoints, 0); - vector fZ(interpPoints, 0); - double stepSize = (rightBound-leftBound)/(interpPoints-1); - for (int i = 0; i < interpPoints; ++i) { - Sw[i] += i*stepSize; - f[i] = max(waterCurve.evaluate(Sw[i]), lowerRelPermLimit) / + vector Sw(interpPoints, leftBound); + vector f(interpPoints, 0); + vector fY(interpPoints, 0); + vector fZ(interpPoints, 0); + double stepSize = (rightBound-leftBound)/(interpPoints-1); + for (int i = 0; i < interpPoints; ++i) { + Sw[i] += i*stepSize; + f[i] = max(waterCurve.evaluate(Sw[i]), lowerRelPermLimit) / max(oilCurve.evaluate(Sw[i]), lowerRelPermLimit) - * viscosities[oilPhaseIndex]/viscosities[waterPhaseIndex]; + * viscosities[oilPhaseIndex]/viscosities[waterPhaseIndex]; if (anisotropic_input) { - fY[i] = max(waterCurveY.evaluate(Sw[i]), lowerRelPermLimit) / + fY[i] = max(waterCurveY.evaluate(Sw[i]), lowerRelPermLimit) / max(oilCurveY.evaluate(Sw[i]), lowerRelPermLimit) - * viscosities[oilPhaseIndex]/viscosities[waterPhaseIndex]; - fZ[i] = max(waterCurveZ.evaluate(Sw[i]), lowerRelPermLimit) / + * viscosities[oilPhaseIndex]/viscosities[waterPhaseIndex]; + fZ[i] = max(waterCurveZ.evaluate(Sw[i]), lowerRelPermLimit) / max(oilCurveZ.evaluate(Sw[i]), lowerRelPermLimit) - * viscosities[oilPhaseIndex]/viscosities[waterPhaseIndex]; + * viscosities[oilPhaseIndex]/viscosities[waterPhaseIndex]; } - } - FracFlowRatioInv.push_back(MonotCubicInterpolator(f,Sw)); + } + FracFlowRatioInv.push_back(MonotCubicInterpolator(f,Sw)); if (anisotropic_input) { FracFlowRatioInvY.push_back(MonotCubicInterpolator(fY,Sw)); FracFlowRatioInvZ.push_back(MonotCubicInterpolator(fY,Sw)); } - } - + } + // Find max/min fracflowratio over all rock types. - double fracflowratioMin = numeric_limits().max(); - double fracflowratioMax = 0; - for (int rockIdx = 0; rockIdx < stone_types; ++rockIdx) { - double fracflowratioMinRock = FracFlowRatioInv[rockIdx].getMinimumX().first; - double fracflowratioMaxRock = FracFlowRatioInv[rockIdx].getMaximumX().first; - fracflowratioMin = min(fracflowratioMin, fracflowratioMinRock); - fracflowratioMax = max(fracflowratioMax, fracflowratioMaxRock); + double fracflowratioMin = numeric_limits().max(); + double fracflowratioMax = 0; + for (int rockIdx = 0; rockIdx < stone_types; ++rockIdx) { + double fracflowratioMinRock = FracFlowRatioInv[rockIdx].getMinimumX().first; + double fracflowratioMaxRock = FracFlowRatioInv[rockIdx].getMaximumX().first; + fracflowratioMin = min(fracflowratioMin, fracflowratioMinRock); + fracflowratioMax = max(fracflowratioMax, fracflowratioMaxRock); if (anisotropic_input) { - double fracflowratioMinRockY = FracFlowRatioInvY[rockIdx].getMinimumX().first; - double fracflowratioMaxRockY = FracFlowRatioInvY[rockIdx].getMaximumX().first; - fracflowratioMin = min(fracflowratioMin, fracflowratioMinRockY); - fracflowratioMax = max(fracflowratioMax, fracflowratioMaxRockY); - double fracflowratioMinRockZ = FracFlowRatioInvZ[rockIdx].getMinimumX().first; - double fracflowratioMaxRockZ = FracFlowRatioInvZ[rockIdx].getMaximumX().first; - fracflowratioMin = min(fracflowratioMin, fracflowratioMinRockZ); - fracflowratioMax = max(fracflowratioMax, fracflowratioMaxRockZ); - } - } - + double fracflowratioMinRockY = FracFlowRatioInvY[rockIdx].getMinimumX().first; + double fracflowratioMaxRockY = FracFlowRatioInvY[rockIdx].getMaximumX().first; + fracflowratioMin = min(fracflowratioMin, fracflowratioMinRockY); + fracflowratioMax = max(fracflowratioMax, fracflowratioMaxRockY); + double fracflowratioMinRockZ = FracFlowRatioInvZ[rockIdx].getMinimumX().first; + double fracflowratioMaxRockZ = FracFlowRatioInvZ[rockIdx].getMaximumX().first; + fracflowratioMin = min(fracflowratioMin, fracflowratioMinRockZ); + fracflowratioMax = max(fracflowratioMax, fracflowratioMaxRockZ); + } + } + if (isMaster) { cout << endl << "Lower fracflowratio: " << fracflowratioMin << ", Upper fracflowratio: " << fracflowratioMax << endl; } - // Now upscale fractional flow vs water saturation + // Now upscale fractional flow vs water saturation // (i.e., populate the vector WaterSaturationVsFractionalFlow) - + double fracFlowRatioTestvalue; - + while (largestSaturationInterval > (Swor-Swir)/500.0) { if (fracflowratioMax == fracflowratioMin) { - // This is a dummy situation, we go through once and then + // This is a dummy situation, we go through once and then // we are finished (this will be triggered by zero permeability) // (not checked for applicatiblity in viscous limit, this is copied from cap limit) fracFlowRatioTestvalue = fracflowratioMin; @@ -1120,7 +1120,7 @@ try fracFlowRatioTestvalue = SatDiff.first; largestSaturationInterval = SatDiff.second; } - + // Check for saneness of fracFlowRatioTestvalue if (std::isnan(fracFlowRatioTestvalue) || std::isinf(fracFlowRatioTestvalue)) { if (isMaster) { @@ -1128,12 +1128,12 @@ try } break; // Jump out out while-loop, just print the results up to now and exit } - + // Do the saturation modelling with the current fracFlowRatioTestvalue double waterVolume = 0.0; vector waterSaturationRockType; waterSaturationRockType.resize(stone_types); - + for (int rockIdx = 0; rockIdx < stone_types; ++rockIdx) { waterSaturationRockType[rockIdx] = FracFlowRatioInv[rockIdx].evaluate(fracFlowRatioTestvalue); if (anisotropic_input) { @@ -1148,7 +1148,7 @@ try // In case we have all flat endpoints in this curve, we chop them off: // (We still preserve Swir and Swor) WaterSaturationVsFractionalFlow.chopFlatEndpoints(saturationThreshold); - + // In case we have a monotone function, but not strictly monotone, // we can remove some data points in order to make it strictly // monotone. This has been seen to happen on one occasion, quite @@ -1156,7 +1156,7 @@ try // increase lowerRelPermLimit from 1e-12 to 1e-8. We might as well // do both. WaterSaturationVsFractionalFlow.shrinkFlatAreas(); - + // Now we can also invert the upscaled water saturation // (it should be monotonic) if (!WaterSaturationVsFractionalFlow.isStrictlyMonotone()) { @@ -1167,10 +1167,10 @@ try } usageandexit(); } - MonotCubicInterpolator FractionalFlowVsWaterSaturation(WaterSaturationVsFractionalFlow.get_fVector(), + MonotCubicInterpolator FractionalFlowVsWaterSaturation(WaterSaturationVsFractionalFlow.get_fVector(), WaterSaturationVsFractionalFlow.get_xVector()); - - + + /***************************************************************************** * Step 7: * Upscale single phase permeability @@ -1184,7 +1184,7 @@ try zero(zeroMatrix); Matrix permTensor = zeroMatrix; Matrix permTensorInv = zeroMatrix; - + if (isMaster) { //cout << "Rank " << mpi_rank << " upscaling single-phase permeability..."; flush(cout); Matrix cellperm = zeroMatrix; @@ -1210,11 +1210,11 @@ try } /***************************************************************** - * Step 8 and 9 - * (step 8 is for water, 9 is for oil, the code is identical, we just - * swap some variables) - * - * For uniformly distributed water saturation values between upscaled Swir and Swor, + * Step 8 and 9 + * (step 8 is for water, 9 is for oil, the code is identical, we just + * swap some variables) + * + * For uniformly distributed water saturation values between upscaled Swir and Swor, * a: Find fractional flow ratio corresponding to wanted upscaled water saturation * b: Model water saturation in each cell given fractional flow ratio * c: Model phase permeability in each cell given cell water saturation and inputted @@ -1223,28 +1223,28 @@ try */ vector WaterSaturation; - + vector > > PhasePerm; // 'phases' * 'tensorElementCount' phaseperm values per fracflowpoint - + vector > tmp1; vector > tmp2; PhasePerm.push_back(tmp1); PhasePerm.push_back(tmp2); - - // Put empty interpolator objects for the upscaled results, number of - // interpolators depending on the boundary condition. + + // Put empty interpolator objects for the upscaled results, number of + // interpolators depending on the boundary condition. for (int idx=0; idx < points; ++idx) { WaterSaturation.push_back(0.0); // Pad with zeroes vector foo1, foo2; PhasePerm[waterPhaseIndex].push_back(foo1); PhasePerm[oilPhaseIndex].push_back(foo2); - for (int voigtIdx=0; voigtIdx < tensorElementCount; ++voigtIdx) { - PhasePerm[waterPhaseIndex][idx].push_back(0.0); // Pad with zeroes. - PhasePerm[oilPhaseIndex][idx].push_back(0.0); - } + for (int voigtIdx=0; voigtIdx < tensorElementCount; ++voigtIdx) { + PhasePerm[waterPhaseIndex][idx].push_back(0.0); // Pad with zeroes. + PhasePerm[oilPhaseIndex][idx].push_back(0.0); + } } - // Make vector of frac flow ratio points corresponding to uniformly distributed + // Make vector of frac flow ratio points corresponding to uniformly distributed // saturation points between Swor and Swir vector fracFlowRatioPoints; for (int pointidx=1; pointidx <= points; ++pointidx) { @@ -1252,7 +1252,7 @@ try double saturation = Swir + (Swor-Swir)/(points-1)*(pointidx-1); fracFlowRatioPoints.push_back(FractionalFlowVsWaterSaturation.evaluate(saturation)); } - + // Construct a vector that tells for each fracflowratio point which mpi-node (rank) should compute for that // particular fracflowratio point vector node_vs_fracflowratiopoint; @@ -1260,7 +1260,7 @@ try for (int idx=0; idx < points; ++idx) { node_vs_fracflowratiopoint.push_back(0); } - + #ifdef HAVE_MPI // Distribute work load over mpi nodes. for (int idx=0; idx < points; ++idx) { @@ -1270,7 +1270,7 @@ try /*if (isMaster) { cout << "Fracflowratio point " << idx << " assigned to node " << node_vs_fracflowratiopoint[idx] << endl; }*/ - } + } #endif clock_t start_upscale_wallclock = clock(); @@ -1279,7 +1279,7 @@ try // Now loop through the vector of fractional flow ratios that // this node should compute. for (int phase = waterPhaseIndex; phase <= oilPhaseIndex; ++phase) { - + string phaseName; if (phase == waterPhaseIndex) { phaseName = string("water"); @@ -1289,23 +1289,23 @@ try } if (isMaster) cout << endl << "Upscaling relative permeability for " << phaseName << "... " << endl; for (int pointidx = 0; pointidx < points; ++pointidx) { - + // Should "I" (mpi-wise) compute this fracflowratio point? if (node_vs_fracflowratiopoint[pointidx] == mpi_rank) { - + fracFlowRatioTestvalue = fracFlowRatioPoints[pointidx]; - + double maxPhasePerm = 0.0; - + vector phasePermValues; vector > phasePermValuesDiag; phasePermValues.resize(satnums.size()); phasePermValuesDiag.resize(satnums.size()); waterVolumeLF = 0.0; - + vector waterSaturationRockType; waterSaturationRockType.resize(stone_types); - + for (int rockIdx = 0; rockIdx < stone_types; ++rockIdx) { waterSaturationRockType[rockIdx] = FracFlowRatioInv[rockIdx].evaluate(fracFlowRatioTestvalue); if (anisotropic_input) { @@ -1314,7 +1314,7 @@ try waterSaturationRockType[rockIdx] /= 3.0; // arithmetic average of three directions. } waterVolumeLF += waterSaturationRockType[rockIdx] * rocktypeVolume[rockIdx]; - } + } for (unsigned int i = 0; i < ecl_idx.size(); ++i) { unsigned int cell_idx = ecl_idx[i]; double cellPhasePerm = minPerm; @@ -1324,12 +1324,12 @@ try cellPhasePermDiag.push_back(minPerm); if (satnums[cell_idx] > 0) { // Satnum zero is "no rock", model those with minPerm. - + // Water saturation is only a function of the rock type - double saturationCell + double saturationCell = waterSaturationRockType[int(satnums[cell_idx])-1]; if (! anisotropic_input) { - double cellRelPerm; + double cellRelPerm; if (phase == waterPhaseIndex) { cellRelPerm = Krw[int(satnums[cell_idx])-1].evaluate(saturationCell); } @@ -1340,19 +1340,19 @@ try } else { //anisotropic if (phase == waterPhaseIndex) { - cellPhasePermDiag[0] = Krwx[int(satnums[cell_idx])-1].evaluate(saturationCell) * + cellPhasePermDiag[0] = Krwx[int(satnums[cell_idx])-1].evaluate(saturationCell) * permxs[cell_idx]; - cellPhasePermDiag[1] = Krwy[int(satnums[cell_idx])-1].evaluate(saturationCell) * + cellPhasePermDiag[1] = Krwy[int(satnums[cell_idx])-1].evaluate(saturationCell) * permys[cell_idx]; - cellPhasePermDiag[2] = Krwz[int(satnums[cell_idx])-1].evaluate(saturationCell) * + cellPhasePermDiag[2] = Krwz[int(satnums[cell_idx])-1].evaluate(saturationCell) * permzs[cell_idx]; } else { // oil - cellPhasePermDiag[0] = Krox[int(satnums[cell_idx])-1].evaluate(saturationCell) * + cellPhasePermDiag[0] = Krox[int(satnums[cell_idx])-1].evaluate(saturationCell) * permxs[cell_idx]; - cellPhasePermDiag[1] = Kroy[int(satnums[cell_idx])-1].evaluate(saturationCell) * + cellPhasePermDiag[1] = Kroy[int(satnums[cell_idx])-1].evaluate(saturationCell) * permys[cell_idx]; - cellPhasePermDiag[2] = Kroz[int(satnums[cell_idx])-1].evaluate(saturationCell) * + cellPhasePermDiag[2] = Kroz[int(satnums[cell_idx])-1].evaluate(saturationCell) * permzs[cell_idx]; } } @@ -1363,12 +1363,12 @@ try maxPhasePerm = max(maxPhasePerm, *max_element(cellPhasePermDiag.begin(), cellPhasePermDiag.end())); } - + // Now we can determine the smallest permitted permeability we can calculate for // We have both a fixed bottom limit, as well as a possible higher limit determined // by a maximum allowable permeability. double minPhasePerm = max(maxPhasePerm/maxPermContrast, minPerm); - + // Now remodel the phase permeabilities obeying minPhasePerm. Matrix cellperm = zeroMatrix; for (unsigned int i = 0; i < ecl_idx.size(); ++i) { @@ -1393,7 +1393,7 @@ try upscaler.setPermeability(i, cellperm); } Matrix phasePermTensor = upscaler.upscaleSinglePhase(); - + // Here we recalculate the upscaled water saturation, // although it is already known when we asked for the // fracflowratio point to compute for. Nonetheless, we @@ -1402,17 +1402,17 @@ try // points are not perfectly uniformly distributed) //cout << waterVolumeLF/poreVolume; WaterSaturation[pointidx] = waterVolumeLF/poreVolume; - + #ifdef HAVE_MPI cout << "Rank " << mpi_rank << ": " << endl;; #endif cout << fracFlowRatioTestvalue << "\t" << WaterSaturation[pointidx]; // Store and print phase-perm-result - for (int voigtIdx=0; voigtIdx < tensorElementCount; ++voigtIdx) { + for (int voigtIdx=0; voigtIdx < tensorElementCount; ++voigtIdx) { PhasePerm[phase][pointidx][voigtIdx] = ::Opm::getVoigtValue(phasePermTensor,voigtIdx); cout << "\t" << ::Opm::getVoigtValue(phasePermTensor,voigtIdx); - } - cout << endl; + } + cout << endl; } } // end loop over saturation points } // end phase loop @@ -1465,19 +1465,19 @@ try #ifdef HAVE_MPI // Sum the upscaling time used by all processes double timeused_total; - MPI_Reduce(&timeused_upscale_wallclock, &timeused_total, 1, MPI_DOUBLE, + MPI_Reduce(&timeused_upscale_wallclock, &timeused_total, 1, MPI_DOUBLE, MPI_SUM, 0, MPI_COMM_WORLD); double avg_upscaling_time_pr_point = timeused_total/(2.0*(double)points); - + #else double avg_upscaling_time_pr_point = timeused_upscale_wallclock / (2.0*(double)points); #endif - /* + /* * Step Xc: Make relperm values from phaseperms * (only master node can do this) */ - + vector > > RelPermValues; // phase is first index voigtIdx is second index. vector > tmp10, tmp20; RelPermValues.push_back(tmp10); @@ -1489,7 +1489,7 @@ try } if (isMaster) { for (int phase = waterPhaseIndex; phase <= oilPhaseIndex; ++phase) { - // Loop over all fracflowratio points + // Loop over all fracflowratio points for (int idx=0; idx < points; ++idx) { Matrix phasePermTensor = zeroMatrix; zero(phasePermTensor); @@ -1516,28 +1516,28 @@ try * file if the '-outputWater'-option and/or '-outputOil' has been set, as this option is an * empty string by default. */ - + if (isMaster) { // If no data computed, we do not have more to do: if (WaterSaturation.size() == 0) { return(1); // non-zero return value, this means something wrong with input data. } - + stringstream outheadtmp; stringstream outwatertmp; stringstream outoiltmp; - + // Print a table of all computed values: outheadtmp << "######################################################################" << endl; outheadtmp << "# Results from upscaling relative permeability."<< endl; outheadtmp << "#" << endl; time_t now = std::time(NULL); outheadtmp << "# Finished: " << asctime(localtime(&now)); - + utsname hostname; uname(&hostname); outheadtmp << "# Hostname: " << hostname.nodename << endl; - + outheadtmp << "#" << endl; outheadtmp << "# Eclipse file: " << ECLIPSEFILENAME << endl; outheadtmp << "# cells: " << tesselatedCells << endl; @@ -1576,11 +1576,11 @@ try outheadtmp << "# points: " << options["points"] << endl; outheadtmp << "# maxPermContrast: " << options["maxPermContrast"] << endl; outheadtmp << "# minPerm: " << options["minPerm"] << endl; - outheadtmp << "# minPoro: " << options["minPoro"] << endl; + outheadtmp << "# minPoro: " << options["minPoro"] << endl; outheadtmp << "# waterViscosity: " << options["waterViscosity"] << " Pa s" << endl; outheadtmp << "# oilViscosity: " << options["oilViscosity"] << " Pa s" << endl; - if (doInterpolate) { - outheadtmp << "# interpolate: " << options["interpolate"] << " points" << endl; + if (doInterpolate) { + outheadtmp << "# interpolate: " << options["interpolate"] << " points" << endl; } outheadtmp << "# " << endl; outheadtmp << "# Single phase permeability" << endl; @@ -1588,10 +1588,10 @@ try outheadtmp << "# |Kyx Kyy Kyz| = " << permTensor(1,0) << " " << permTensor(1,1) << " " << permTensor(1,2) << endl; outheadtmp << "# |Kzx Kzy Kzz| = " << permTensor(2,0) << " " << permTensor(2,1) << " " << permTensor(2,2) << endl; outheadtmp << "# " << endl; - if (doInterpolate) { - outheadtmp << "# NB: Data points shown are interpolated." << endl; + if (doInterpolate) { + outheadtmp << "# NB: Data points shown are interpolated." << endl; } outheadtmp << "######################################################################" << endl; - + if (isFixed) { outwatertmp << "# v_w/v_o Sw Krwxx Krwyy Krwzz" << endl; outoiltmp << "# v_w/v_o Sw Kroxx Kroyy Krozz" << endl; @@ -1604,106 +1604,106 @@ try outwatertmp << "# v_w/v_o Sw Krwxx Krwyy Krwzz Krwyz Krwxz Krwxy Krwzy Krwzx Krwyx" << endl; outoiltmp << "# v_w/v_o Sw Kroxx Kroyy Krozz Kroyz Kroxz Kroxy Krozy Krozx Kroyx" << endl; } - - // If user wants interpolated output, do monotone cubic interpolation - // by modifying the data vectors that are to be printed + + // If user wants interpolated output, do monotone cubic interpolation + // by modifying the data vectors that are to be printed if (doInterpolate) { - // Find min and max for saturation values + // Find min and max for saturation values double satmin = +DBL_MAX; double satmax = -DBL_MAX; - for (unsigned int i = 0; i < WaterSaturation.size(); ++i) { - if (WaterSaturation[i] < satmin) { - satmin = WaterSaturation[i]; - } - if (WaterSaturation[i] > satmax) { - satmax = WaterSaturation[i]; - } - } + for (unsigned int i = 0; i < WaterSaturation.size(); ++i) { + if (WaterSaturation[i] < satmin) { + satmin = WaterSaturation[i]; + } + if (WaterSaturation[i] > satmax) { + satmax = WaterSaturation[i]; + } + } // Make uniform grid in saturation axis (for both water and oil) vector SatvaluesInterp, fracFlowRatioPointsInterp; - for (int i = 0; i < interpolationPoints; ++i) { - SatvaluesInterp.push_back(satmin + ((double)i)/((double)interpolationPoints-1)*(satmax-satmin)); + for (int i = 0; i < interpolationPoints; ++i) { + SatvaluesInterp.push_back(satmin + ((double)i)/((double)interpolationPoints-1)*(satmax-satmin)); fracFlowRatioPointsInterp.push_back(WaterSaturationVsFractionalFlow.evaluate(SatvaluesInterp[i])); - } + } - // Now fracflowratio and computed relperm-values must be viewed as functions - // of saturation, and then interpolated on the uniform saturation grid. - - // Now overwrite existing FlowRatioValues and relperm-data with interpolated data: + // Now fracflowratio and computed relperm-values must be viewed as functions + // of saturation, and then interpolated on the uniform saturation grid. + + // Now overwrite existing FlowRatioValues and relperm-data with interpolated data: MonotCubicInterpolator FracFlowRatioValuesVsSaturationWater, FracFlowRatioValuesVsSaturationOil; - - for (int voigtIdx = 0; voigtIdx < tensorElementCount; ++voigtIdx) { - MonotCubicInterpolator RelPermWaterVsSaturation(WaterSaturation, RelPermValues[waterPhaseIndex][voigtIdx]); - MonotCubicInterpolator RelPermOilVsSaturation( WaterSaturation, RelPermValues[oilPhaseIndex][voigtIdx]); - RelPermValues[waterPhaseIndex][voigtIdx].clear(); + + for (int voigtIdx = 0; voigtIdx < tensorElementCount; ++voigtIdx) { + MonotCubicInterpolator RelPermWaterVsSaturation(WaterSaturation, RelPermValues[waterPhaseIndex][voigtIdx]); + MonotCubicInterpolator RelPermOilVsSaturation( WaterSaturation, RelPermValues[oilPhaseIndex][voigtIdx]); + RelPermValues[waterPhaseIndex][voigtIdx].clear(); RelPermValues[oilPhaseIndex][voigtIdx].clear(); - for (int i=0; i < interpolationPoints; ++i) { - RelPermValues[waterPhaseIndex][voigtIdx].push_back(RelPermWaterVsSaturation.evaluate(SatvaluesInterp[i])); - RelPermValues[oilPhaseIndex][voigtIdx].push_back(RelPermOilVsSaturation.evaluate(SatvaluesInterp[i])); - } - } - + for (int i=0; i < interpolationPoints; ++i) { + RelPermValues[waterPhaseIndex][voigtIdx].push_back(RelPermWaterVsSaturation.evaluate(SatvaluesInterp[i])); + RelPermValues[oilPhaseIndex][voigtIdx].push_back(RelPermOilVsSaturation.evaluate(SatvaluesInterp[i])); + } + } + // Now also overwrite Satvalues and fracFlowRatioPoints - WaterSaturation.clear(); - WaterSaturation = SatvaluesInterp; - + WaterSaturation.clear(); + WaterSaturation = SatvaluesInterp; + fracFlowRatioPoints.clear(); fracFlowRatioPoints = fracFlowRatioPointsInterp; } - - + + for (unsigned int i=0; i < WaterSaturation.size(); ++i) { outwatertmp << showpoint << setw(14) << fracFlowRatioPoints[i]; outwatertmp << showpoint << setw(14) << WaterSaturation[i]; - + for (int voigtIdx=0; voigtIdx < tensorElementCount; ++voigtIdx) { - outwatertmp << showpoint << setw(14) << RelPermValues[waterPhaseIndex][voigtIdx][i]; + outwatertmp << showpoint << setw(14) << RelPermValues[waterPhaseIndex][voigtIdx][i]; } outwatertmp << endl; - + // Ignore further output if we have reached no water saturation - if (WaterSaturation[i] == 0.0) break; /* maybe == on floats will fail at some point */ + if (WaterSaturation[i] == 0.0) break; /* maybe == on floats will fail at some point */ } - + for (unsigned int i=0; i < WaterSaturation.size(); ++i) { outoiltmp << showpoint << setw(14) << fracFlowRatioPoints[i]; - outoiltmp << showpoint << setw(14) << WaterSaturation[i]; - + outoiltmp << showpoint << setw(14) << WaterSaturation[i]; + for (int voigtIdx=0; voigtIdx < tensorElementCount; ++voigtIdx) { - outoiltmp << showpoint << setw(14) << RelPermValues[oilPhaseIndex][voigtIdx][i]; + outoiltmp << showpoint << setw(14) << RelPermValues[oilPhaseIndex][voigtIdx][i]; } outoiltmp << endl; - + // Ignore further output if we have reached no water saturation - if (WaterSaturation[i] == 0.0) break; /* maybe == on floats will fail at some point */ + if (WaterSaturation[i] == 0.0) break; /* maybe == on floats will fail at some point */ } - + cout << outheadtmp.str(); cout << "# Relperm curve for water: " << endl; cout << outwatertmp.str(); cout << "####################################################################" << endl; cout << "# Relperm curve for oil: " << endl; cout << outoiltmp.str(); - - + + if (options["outputWater"] != "") { cout << "Writing (water) results to " << options["outputWater"] << endl; ofstream outfile; outfile.open(options["outputWater"].c_str(), ios::out | ios::trunc); outfile << outheadtmp.str(); outfile << outwatertmp.str(); - outfile.close(); + outfile.close(); } - + if (options["outputOil"] != "") { cout << "Writing (oil) results to " << options["outputOil"] << endl; ofstream outfile; outfile.open(options["outputOil"].c_str(), ios::out | ios::trunc); outfile << outheadtmp.str(); outfile << outoiltmp.str(); - outfile.close(); + outfile.close(); } } @@ -1713,4 +1713,3 @@ catch (const std::exception &e) { std::cerr << "Program threw an exception: " << e.what() << "\n"; throw; } - diff --git a/examples/upscale_steadystate_implicit.cpp b/examples/upscale_steadystate_implicit.cpp index e62f630..1b8cbd2 100644 --- a/examples/upscale_steadystate_implicit.cpp +++ b/examples/upscale_steadystate_implicit.cpp @@ -126,11 +126,11 @@ std::vector > getExtremeSats(std::string rock_list, std::vec if (!rock_stream) { OPM_THROW(std::runtime_error, "Could not open file " + rockfilename); } - + if (! anisorocks) { //Isotropic input rocks (Sw Krw Kro J) MonotCubicInterpolator Jtmp; try { - Jtmp = MonotCubicInterpolator(rockname, 1, 4); + Jtmp = MonotCubicInterpolator(rockname, 1, 4); } catch (const char * errormessage) { std::cerr << "Error: " << errormessage << std::endl; @@ -159,7 +159,7 @@ std::vector > getExtremeSats(std::string rock_list, std::vec } rocksatendp[i][0] = Pctmp.getMinimumX().first; rocksatendp[i][1] = Pctmp.getMaximumX().first; - } + } } return rocksatendp; } @@ -188,12 +188,12 @@ try std::string gridfilename = param.get("gridfilename"); auto deck = Opm::RelPermUpscaleHelper::parseEclipseFile(gridfilename); - // Check that we have the information we need from the eclipse file: + // Check that we have the information we need from the eclipse file: if (! (deck.hasKeyword("SPECGRID") && deck.hasKeyword("COORD") && deck.hasKeyword("ZCORN") && deck.hasKeyword("PORO") && deck.hasKeyword("PERMX"))) { - std::cerr << "Error: Did not find SPECGRID, COORD, ZCORN, PORO and PERMX in Eclipse file " << gridfilename << std::endl; - usageandexit(); - } + std::cerr << "Error: Did not find SPECGRID, COORD, ZCORN, PORO and PERMX in Eclipse file " << gridfilename << std::endl; + usageandexit(); + } // Set default values if not given as input double linsolver_tolerance = param.getDefault("residual_tolerance", 1e-8); @@ -253,25 +253,25 @@ try std::cerr << "Boundary conditions (bc) must be either fixed, periodic or linear. Value is " << bc << std::endl; usageandexit(); } - + // Compute minimum and maximum (large scale) saturations std::string rock_list = param.getDefault("rock_list", "no_list"); std::vector rockfiles; std::vector > rocksatendpoints_ = getExtremeSats(rock_list,rockfiles); - std::vector poros = deck["PORO"].back().getSIDoubleData(); + std::vector poros = deck["PORO"].back().getSIDoubleData(); // Anisotropic relperm not yet implemented in steadystate_implicit //bool anisorocks = param.getDefault("anisotropicrocks", false); - std::vector satnums(poros.size(), 1); + std::vector satnums(poros.size(), 1); if (deck.hasKeyword("SATNUM")) { satnums = deck["SATNUM"].back().getIntData(); - } + } else if (deck.hasKeyword("ROCKTYPE")) { satnums = deck["ROCKTYPE"].back().getIntData(); - } - else { - std::cout << "Warning: SATNUM or ROCKTYPE not found in input file, assuming only one rocktype" << std::endl; - } + } + else { + std::cout << "Warning: SATNUM or ROCKTYPE not found in input file, assuming only one rocktype" << std::endl; + } // check that number of rock types in rock_list matches number of rock types in grid int num_rock_types_grid = int(*(max_element(satnums.begin(), satnums.end()))); int num_rock_types_rocklist = rocksatendpoints_.size(); @@ -284,7 +284,7 @@ try spupscaler.init(deck, Opm::SinglePhaseUpscaler::Fixed, 0.0, linsolver_tolerance, linsolver_verbosity, linsolver_type, false, linsolver_maxit, linsolver_prolongate_factor, linsolver_smooth_steps); - std::vector cellPoreVolumes; + std::vector cellPoreVolumes; cellPoreVolumes.resize(satnums.size(), 0.0); double swirvolume = 0.0; double sworvolume = 0.0; @@ -300,7 +300,7 @@ try } } // Total porevolume and total volume -> upscaled porosity: - double poreVolume = std::accumulate(cellPoreVolumes.begin(), + double poreVolume = std::accumulate(cellPoreVolumes.begin(), cellPoreVolumes.end(), 0.0); double min_sat = swirvolume/poreVolume; @@ -308,7 +308,7 @@ try // Insert computed Swir and Swor as min_sat and max_sat in param object //param.insertParameter("min_sat",toString(min_sat)); //param.insertParameter("max_sat",toString(max_sat)); - + std::vector saturations; Opm::SparseTable all_pdrops; // Linear range of saturations @@ -332,21 +332,21 @@ try typedef Upscaler::permtensor_t permtensor_t; Upscaler upscaler; upscaler.init(param); - + // Compute single phase permeability permtensor_t upscaled_K = upscaler.upscaleSinglePhase(); permtensor_t singlephaseperm = upscaled_K; singlephaseperm *= (1.0/(Opm::prefix::milli*Opm::unit::darcy)); double porosity = upscaler.upscalePorosity(); int num_cells = upscaler.grid().size(0); - + // Holders for upscaled relpermvalues for each saturation-pressure-point std::vector > > RelPermPhase1, RelPermPhase2; for (int satidx=0; satidx >(all_pdrops[satidx].size(),std::vector(tensorElementCount,0.0))); RelPermPhase2.push_back(std::vector >(all_pdrops[satidx].size(),std::vector(tensorElementCount,0.0))); } - + for (int satidx = 0; satidx < num_sats; ++satidx) { std::vector init_sat; if (start_from_cl) { @@ -390,7 +390,7 @@ try outputtmp << "#" << std::endl; time_t now = std::time(NULL); outputtmp << "# Finished: " << asctime(localtime(&now)); - + utsname hostname; uname(&hostname); outputtmp << "# Hostname: " << hostname.nodename << std::endl; outputtmp << "#" << std::endl; @@ -400,7 +400,7 @@ try outputtmp << "# Porosity: " << porosity << std::endl; outputtmp << "#" << std::endl; outputtmp << "# Rock list: " << rock_list << std::endl; - outputtmp << "# with the following rock files: " << std::endl; + outputtmp << "# with the following rock files: " << std::endl; for (std::size_t ridx=0; ridx upscaler_t; //SteadyStateUpscalerManagerImplicit mgr; //mgr.upscale(param); - + } catch (const std::exception &e) { std::cerr << "Program threw an exception: " << e.what() << "\n"; diff --git a/opm/elasticity/asmhandler.hpp b/opm/elasticity/asmhandler.hpp index a1fe4a6..12d3403 100644 --- a/opm/elasticity/asmhandler.hpp +++ b/opm/elasticity/asmhandler.hpp @@ -96,14 +96,14 @@ class ASMHandler { } //! \brief This function needs to be called before starting - //! the element assembly. + //! the element assembly. void initForAssembly(); //! \brief Add an element matrix/vector to the system //! \param[in] K Pointer to the element matrix. Can be NULL //! \param[in] S Pointer to the element load vector. Can be NULL //! \param[in] cell An iterator pointing to the cell we're assembling for - //! \param[in] b Vector to add contributions to. If not given, + //! \param[in] b Vector to add contributions to. If not given, //! we use the internal vector template void addElement(const Dune::FieldMatrix* K, @@ -121,7 +121,7 @@ class ASMHandler { //! \param[in] it An iterator to the cell we want to extract values for //! \param[out] v Vector holding the values requested template - void extractValues(Dune::FieldVector& v, + void extractValues(Dune::FieldVector& v, const Vector& u, const LeafIterator& it); //! \brief Expand a system vector to a solution vector @@ -141,7 +141,7 @@ class ASMHandler { //! \brief Update/add a fixed node //! \param[in] node The node number //! \param[in] entry The fixed values - void updateFixedNode(int node, + void updateFixedNode(int node, const std::pair& entry); //! \brief Check if a node is marked as fixed (in any direction) @@ -185,7 +185,7 @@ class ASMHandler { //! \param[in] row The equation number/row in matrix void nodeAdjacency(const LeafIterator& it, int vertexsize, int row); - //! \brief Internal function. Calculate adjacency pattern + //! \brief Internal function. Calculate adjacency pattern void determineAdjacencyPattern(); //! \brief Internal function. Assemble entries for a single DOF diff --git a/opm/elasticity/asmhandler_impl.hpp b/opm/elasticity/asmhandler_impl.hpp index d23dd4c..11421fd 100644 --- a/opm/elasticity/asmhandler_impl.hpp +++ b/opm/elasticity/asmhandler_impl.hpp @@ -53,7 +53,7 @@ void ASMHandler::addDOF(int row, int erow, const Dune::FieldMatrix* K, const Dune::FieldVector* S, const LeafIndexSet& set, - const LeafIterator& cell, + const LeafIterator& cell, Vector* bptr, double scale) { @@ -192,7 +192,7 @@ void ASMHandler::addMPC(MPC* mpc) int slaveNode = mpc->getSlave().node*dim+mpc->getSlave().dof-1; fixIt it = fixedNodes.find(mpc->getSlave().node); int flag = 1 << (mpc->getSlave().dof-1); - if (it == fixedNodes.end() || + if (it == fixedNodes.end() || !(it->second.first & flag)) { mpcs.insert(std::make_pair(slaveNode,mpc)); return; diff --git a/opm/elasticity/boundarygrid.cpp b/opm/elasticity/boundarygrid.cpp index f26d322..45ae693 100644 --- a/opm/elasticity/boundarygrid.cpp +++ b/opm/elasticity/boundarygrid.cpp @@ -116,7 +116,7 @@ inline bool EQUAL2(const BoundaryGrid::FaceCoord& x, bool BoundaryGrid::find(Vertex& res, const Vertex& coord) const { - // find first quad with coord within bounding box + // find first quad with coord within bounding box std::vector::const_iterator it = std::find_if(grid.begin(),grid.end(), BoundedPredicate(coord.c)); @@ -131,7 +131,7 @@ bool BoundaryGrid::find(Vertex& res, const Vertex& coord) const ok = true; break; } - } + } if (!ok && Q4inv(res.c,*it,coord.c,1.e-8,1.e-8) > 0) { ok = true; } @@ -351,7 +351,7 @@ std::vector BoundaryGrid::Quad::evalBasis(double xi, double eta) const BoundaryGrid::Vertex minXminY(const std::vector& in) { // find the nodes with minimal X - // then find the minimum Y among these + // then find the minimum Y among these std::vector s(in); std::sort(s.begin(),s.end(),BoundaryGrid::VertexLess(0)); std::sort(s.begin(),s.begin()+2,BoundaryGrid::VertexLess(1)); @@ -361,7 +361,7 @@ BoundaryGrid::Vertex minXminY(const std::vector& in) BoundaryGrid::Vertex maxXminY(const std::vector& in) { // find the nodes with maximum X - // then find the minimum Y among these + // then find the minimum Y among these std::vector s(in); std::sort(s.begin(),s.end(),BoundaryGrid::VertexLess(0)); std::sort(s.begin()+2,s.end(),BoundaryGrid::VertexLess(1)); @@ -371,7 +371,7 @@ BoundaryGrid::Vertex maxXminY(const std::vector& in) BoundaryGrid::Vertex maxXmaxY(const std::vector& in) { // find the nodes with maximum X - // then find the maximum Y among these + // then find the maximum Y among these std::vector s(in); std::sort(s.begin(),s.end(),BoundaryGrid::VertexLess(0)); std::sort(s.begin()+2,s.end(),BoundaryGrid::VertexLess(1)); @@ -381,7 +381,7 @@ BoundaryGrid::Vertex maxXmaxY(const std::vector& in) BoundaryGrid::Vertex minXmaxY(const std::vector& in) { // find the nodes with minimum X - // then find the maximum Y among these + // then find the maximum Y among these std::vector s(in); std::sort(s.begin(),s.end(),BoundaryGrid::VertexLess(0)); std::sort(s.begin(),s.begin()+2,BoundaryGrid::VertexLess(1)); diff --git a/opm/elasticity/boundarygrid.hh b/opm/elasticity/boundarygrid.hh index 772615d..c221642 100644 --- a/opm/elasticity/boundarygrid.hh +++ b/opm/elasticity/boundarygrid.hh @@ -54,7 +54,7 @@ class BoundaryGrid { //! \brief Default (empty) destructor virtual ~BoundaryGrid() {} - //! \brief Holds the indices and relevant coordinates of the vertices + //! \brief Holds the indices and relevant coordinates of the vertices // on a boundary class Quad; @@ -221,7 +221,7 @@ class BoundaryGrid { bool operator()(const Quad& q) { double eps = 1.e-8; - return (coord[0] >= q.bb[0][0]-eps && + return (coord[0] >= q.bb[0][0]-eps && coord[0] <= q.bb[1][0]+eps && coord[1] >= q.bb[0][1]-eps && coord[1] <= q.bb[1][1]+eps); @@ -463,7 +463,7 @@ class HexGeometry<2, cdim, GridImp> //! \brief Returns the inverse, transposed Jacobian //! \param[in] local The local coordinates - const Dune::FieldMatrix + const Dune::FieldMatrix jacobianInverseTransposed(const LocalCoordinate& local) const { Dune::FieldMatrix Jti = jacobianTransposed(local); diff --git a/opm/elasticity/elasticity.hpp b/opm/elasticity/elasticity.hpp index a2c3751..ec1d828 100644 --- a/opm/elasticity/elasticity.hpp +++ b/opm/elasticity/elasticity.hpp @@ -37,7 +37,7 @@ class Elasticity { template void getBVector(Dune::FieldVector& BVector, const Dune::FieldVector& point); - + //! \brief Returns the B matrix in a quadrature point //! \param[out] B The B matrix //! \param[in] point (Reference) coordinates of quadrature point @@ -66,7 +66,7 @@ class Elasticity { //! \param[out] sigma The stress vector in the given quadrature point template void getStressVector(Dune::FieldVector& sigma, - const Dune::FieldVector& v, + const Dune::FieldVector& v, const Dune::FieldVector& eps0, const Dune::FieldMatrix& B, const Dune::FieldMatrix& C); diff --git a/opm/elasticity/elasticity_impl.hpp b/opm/elasticity/elasticity_impl.hpp index 2cc00ae..e0bc013 100644 --- a/opm/elasticity/elasticity_impl.hpp +++ b/opm/elasticity/elasticity_impl.hpp @@ -23,7 +23,7 @@ void Elasticity::getBmatrix(Dune::FieldMatrix& point, const Dune::FieldMatrix& Jinv) { - P1ShapeFunctionSet basis + P1ShapeFunctionSet basis = P1ShapeFunctionSet::instance(); int funcs = funcdim/dim; @@ -77,7 +77,7 @@ void Elasticity::getBmatrix(Dune::FieldMatrix::getBVector(Dune::FieldVector& Bvector, const Dune::FieldVector& point) { - P1ShapeFunctionSet basis + P1ShapeFunctionSet basis = P1ShapeFunctionSet::instance(); Dune::FieldMatrix N; @@ -85,7 +85,7 @@ void Elasticity::getBVector(Dune::FieldVector& Bvector Bvector[i] = basis[i].evaluateFunction(point); } } - + template template void Elasticity::getStiffnessMatrix( @@ -120,7 +120,7 @@ Dune::FieldVector waveSpeeds(const Dune::FieldMatrix& C, d x[0] = r*cos(theta)*cos(phi); x[1] = r*sin(theta)*cos(phi); x[2] = r*sin(phi); - + Dune::FieldMatrix D; D[0][0] = x[0]; D[0][4] = x[2]; diff --git a/opm/elasticity/elasticity_preconditioners.hpp b/opm/elasticity/elasticity_preconditioners.hpp index f38f514..3e27986 100644 --- a/opm/elasticity/elasticity_preconditioners.hpp +++ b/opm/elasticity/elasticity_preconditioners.hpp @@ -134,7 +134,7 @@ struct AMG1 { } }; -//! \brief A FastAMG +//! \brief A FastAMG struct FastAMG { typedef Dune::Amg::FastAMG type; @@ -167,7 +167,7 @@ struct AMG2Level { typename AMG1::Criterion> CoarsePolicy; typedef typename Dune::Amg::SmootherTraits::Arguments SmootherArgs; - + typedef Dune::Amg::TwoLevelMethod type; //! \brief Setup preconditioner diff --git a/opm/elasticity/elasticity_upscale.hpp b/opm/elasticity/elasticity_upscale.hpp index 730a957..316c90f 100644 --- a/opm/elasticity/elasticity_upscale.hpp +++ b/opm/elasticity/elasticity_upscale.hpp @@ -51,7 +51,7 @@ namespace Elasticity { //! \brief An enumeration of available linear solver classes enum Solver { DIRECT, - ITERATIVE + ITERATIVE }; enum Preconditioner { @@ -225,7 +225,7 @@ class ElasticityUpscale //! \brief Vector holding the volume fractions for materials (grouped by SATNUM) std::vector volumeFractions; //! \brief Are volume fractions grouped by SATNUM? - bool bySat; + bool bySat; //! \brief Upscaled density double upscaledRho; @@ -237,7 +237,7 @@ class ElasticityUpscale //! \param[in] file The eclipse grid file //! \param[in] rocklist If not blank, file is a rocklist //! \param[in] verbose_ If true, give verbose output - ElasticityUpscale(const GridType& gv_, ctype tol_, ctype Escale_, + ElasticityUpscale(const GridType& gv_, ctype tol_, ctype Escale_, const std::string& file, const std::string& rocklist, bool verbose_) : A(gv_), gv(gv_), tol(tol_), Escale(Escale_), E(gv_), verbose(verbose_), @@ -324,8 +324,8 @@ class ElasticityUpscale //! \param[in] value The constant coordinate describing the plane bool isOnPlane(Direction plane, GlobalCoordinate coord, ctype value); - //! \brief Check if the given coordinate falls on a given line - //! \param[in] dir The line direction + //! \brief Check if the given coordinate falls on a given line + //! \param[in] dir The line direction //! \param[in] coord The coordinates to check //! \param[in] x The first coordinate of the line //! \param[in] y The second coordinate of the line diff --git a/opm/elasticity/elasticity_upscale_impl.hpp b/opm/elasticity/elasticity_upscale_impl.hpp index 7e0e3c8..586dded 100644 --- a/opm/elasticity/elasticity_upscale_impl.hpp +++ b/opm/elasticity/elasticity_upscale_impl.hpp @@ -27,7 +27,7 @@ namespace Elasticity { #define IMPL_FUNC(A,B) template \ A ElasticityUpscale::B -IMPL_FUNC(std::vector, +IMPL_FUNC(std::vector, extractFace(Direction dir, ctype coord)) { std::vector result; @@ -67,13 +67,13 @@ IMPL_FUNC(BoundaryGrid, extractMasterFace(Direction dir, int c = 0; int i = log2(float(dir)); BoundaryGrid result; - // we first group nodes into this map through the coordinate of lower left + // we first group nodes into this map through the coordinate of lower left // vertex. we then split this up into pillars for easy processing later std::map > nodeMap; - for (LeafIterator cell = gv.leafGridView().template begin<0>(); + for (LeafIterator cell = gv.leafGridView().template begin<0>(); cell != gv.leafGridView().template end<0>(); ++cell, ++c) { std::vector verts; - int idx=0; + int idx=0; if (side == LEFT) idx = set.subIndex(*cell,V1[i][0],dim); else if (side == RIGHT) @@ -283,7 +283,7 @@ IMPL_FUNC(void, assembleBBlockMortar(const BoundaryGrid& b1, double alpha)) { // get a set of P1 shape functions for the displacements - P1ShapeFunctionSet ubasis = + P1ShapeFunctionSet ubasis = P1ShapeFunctionSet::instance(); // get a set of PN shape functions for the multipliers @@ -298,7 +298,7 @@ IMPL_FUNC(void, assembleBBlockMortar(const BoundaryGrid& b1, // get a quadrature rule int quadorder = std::max((1.0+n1+0.5)/2.0,(1.0+n2+0.5)/2.0); quadorder = std::max(quadorder, 2); - const Dune::QuadratureRule& rule = + const Dune::QuadratureRule& rule = Dune::QuadratureRules::rule(gt,quadorder); // do the assembly loop @@ -369,7 +369,7 @@ IMPL_FUNC(void, fixPoint(Direction dir, typedef typename GridType::LeafGridView::template Codim::Iterator VertexLeafIterator; const VertexLeafIterator itend = gv.leafGridView().template end(); - // make a mapper for codim 0 entities in the leaf grid + // make a mapper for codim 0 entities in the leaf grid using LeafGridView = Dune::GridView>; Dune::MultipleCodimMultipleGeomTypeMapper mapper(gv.leafGridView(), Dune::mcmgVertexLayout()); @@ -682,8 +682,8 @@ IMPL_FUNC(void, loadMaterialsFromGrid(const std::string& file)) std::cout << "Number of materials: " << cache.size() << std::endl; double totalvolume=0; - for (std::map::iterator it = volume.begin(); - it != volume.end(); ++it) + for (std::map::iterator it = volume.begin(); + it != volume.end(); ++it) totalvolume += it->second; // statistics @@ -813,7 +813,7 @@ IMPL_FUNC(void, periodicBCs(const double* min, const double* max)) periodicPlane(Z,XYZ,slave[2],master[2]); } -IMPL_FUNC(void, periodicBCsMortar(const double* min, +IMPL_FUNC(void, periodicBCsMortar(const double* min, const double* max, int n1, int n2, int p1, int p2)) @@ -832,7 +832,7 @@ IMPL_FUNC(void, periodicBCsMortar(const double* min, // step 1 fixCorners(min,max); - + // step 2 std::cout << "\textracting nodes on top face..." << std::endl; slave.push_back(extractFace(Z,max[2])); @@ -911,7 +911,7 @@ IMPL_FUNC(void, setupSolvers(const LinSolParams& params)) #ifdef HAVE_OPENMP numsolvers = omp_get_max_threads(); #endif - + if (params.type == ITERATIVE) { op.reset(new Operator(A.getOperator())); bool copy; @@ -987,7 +987,7 @@ IMPL_FUNC(void, setupSolvers(const LinSolParams& params)) for (Matrix::ConstRowIterator it = A.getOperator().begin(); it != A.getOperator().end(); ++it, ++row) { double alpha=0; - for (Matrix::ConstColIterator it2 = it->begin(); + for (Matrix::ConstColIterator it2 = it->begin(); it2 != it->end(); ++it2) { if (it2.index() != row) alpha += fabs(*it2); @@ -1051,7 +1051,7 @@ IMPL_FUNC(void, setupSolvers(const LinSolParams& params)) } } } else { - if (B.N()) + if (B.N()) A.getOperator() = MatrixOps::augment(A.getOperator(), B, 0, A.getOperator().M(), true); #if HAVE_SUITESPARSE_UMFPACK || HAVE_SUPERLU diff --git a/opm/elasticity/logutils.hpp b/opm/elasticity/logutils.hpp index 427eeef..b354a2b 100644 --- a/opm/elasticity/logutils.hpp +++ b/opm/elasticity/logutils.hpp @@ -42,7 +42,7 @@ class LoggerHelper { } //! \brief Returns the start and end offsets of a chunk group - //! \param[in] group The requested chunk + //! \param[in] group The requested chunk std::pair group(int group) { return std::make_pair(groups[group],groups[group+1]); diff --git a/opm/elasticity/material.cpp b/opm/elasticity/material.cpp index f24f6f3..b67725a 100644 --- a/opm/elasticity/material.cpp +++ b/opm/elasticity/material.cpp @@ -72,13 +72,13 @@ Material* Material::create(int ID, const std::string& file) } std::string str,str2; f >> str; - if (str == "ti") { // transverse isotropic material with unit axis in z-direction. + if (str == "ti") { // transverse isotropic material with unit axis in z-direction. Dune::DynamicVectormat(21,double(0)); double a = 0.0, b = 0.0, c = 0.0, d = 0.0, e = 0.0, rho = 0.0; f >> a >> b >> c >> d >> e; f >> str2; if (str2 != "density") { - std::cerr << "Rock file " << file << " has wrong format. when isotropycase is \"ti\", " << std::endl + std::cerr << "Rock file " << file << " has wrong format. when isotropycase is \"ti\", " << std::endl << "the next line need 5 a,b,c,d and e for C matrix with C11=C22=a, C33=c, " << std::endl << "C44=C55=d, C66=e, C12=a-2e and C13=b. Then keyword density on next " << std::endl << "line followed by material density (double) on last line." << std::endl; @@ -95,15 +95,15 @@ Material* Material::create(int ID, const std::string& file) double rho; for (int cuidx=0; cuidx<21; cuidx++) { f >> mat[cuidx]; - } + } f >> str2; if (str2 != "density") { - std::cerr << "Rock file " << file << " has wrong format. when isotropycase is \"anisotropic\", " << std::endl + std::cerr << "Rock file " << file << " has wrong format. when isotropycase is \"anisotropic\", " << std::endl << "the next line need 21 entries for upper C matrix. Then keyword density on next " << std::endl << "line followed by material density (double) on last line." << std::endl; exit(1); } - + f >> rho; return new OrthotropicSym(ID,mat); } @@ -141,7 +141,7 @@ Material* Material::create(int ID, const std::string& file) std::cerr << "Material in " << file << " is not isotropic (nu=" << nu << "), bailing" << std::endl; exit(1); } - + return new Isotropic(ID,E,nu,rho); } } diff --git a/opm/elasticity/material.hh b/opm/elasticity/material.hh index b77ef1b..9acf452 100644 --- a/opm/elasticity/material.hh +++ b/opm/elasticity/material.hh @@ -39,7 +39,7 @@ protected: //! \brief Prints the material properties to a stream. virtual std::ostream& write(std::ostream& os) const - { + { return os; } diff --git a/opm/elasticity/matrixops.cpp b/opm/elasticity/matrixops.cpp index 0c39282..d3fefa9 100644 --- a/opm/elasticity/matrixops.cpp +++ b/opm/elasticity/matrixops.cpp @@ -76,7 +76,7 @@ void MatrixOps::print(const Matrix& A) it != A.end(); ++it) { for (Matrix::ConstColIterator it2 = it->begin(); it2 != it->end();++it2) { - double val = *it2; + double val = *it2; if (fabs(val) < 1.e-14) continue; std::cout << it.index() << " " << it2.index() << " : " << val << std::endl; @@ -140,14 +140,14 @@ Matrix MatrixOps::augment(const Matrix& A, const Matrix& B, adj.resize(nrow); for (Matrix::ConstIterator it = A.begin(); it != A.end();++it) { - for (Matrix::ConstColIterator it2 = it->begin(); + for (Matrix::ConstColIterator it2 = it->begin(); it2 != it->end();++it2) { adj[it.index()].insert(it2.index()); } } for (Matrix::ConstIterator it = B.begin(); it != B.end();++it) { - for (Matrix::ConstColIterator it2 = it->begin(); + for (Matrix::ConstColIterator it2 = it->begin(); it2 != it->end();++it2) { adj[it.index()+r0].insert(it2.index()+c0); if (symmetric) @@ -163,21 +163,21 @@ Matrix MatrixOps::augment(const Matrix& A, const Matrix& B, fromAdjacency(result,adj,nrow,ncol); for (Matrix::ConstIterator it = A.begin(); it != A.end();++it) { - for (Matrix::ConstColIterator it2 = it->begin(); + for (Matrix::ConstColIterator it2 = it->begin(); it2 != it->end();++it2) { result[it.index()][it2.index()] = *it2; } } for (Matrix::ConstIterator it = B.begin(); it != B.end();++it) { - for (Matrix::ConstColIterator it2 = it->begin(); + for (Matrix::ConstColIterator it2 = it->begin(); it2 != it->end();++it2) { result[it.index()+r0][it2.index()+c0] = *it2; if (symmetric) result[it2.index()+c0][it.index()+r0] = *it2; } } - + return result; } diff --git a/opm/elasticity/matrixops.hpp b/opm/elasticity/matrixops.hpp index f1dc8a2..bb8e4c2 100644 --- a/opm/elasticity/matrixops.hpp +++ b/opm/elasticity/matrixops.hpp @@ -14,7 +14,7 @@ #include -#include +#include #include #include diff --git a/opm/elasticity/mortar_schur_precond.hpp b/opm/elasticity/mortar_schur_precond.hpp index 1c7b0ed..982011f 100644 --- a/opm/elasticity/mortar_schur_precond.hpp +++ b/opm/elasticity/mortar_schur_precond.hpp @@ -25,7 +25,7 @@ namespace Elasticity { * mortar-elasticity system * [A B] * [B' ] - * + * * The preconditioner is * [Apre B] * [ P] @@ -118,7 +118,7 @@ class MortarSchurPre : public Dune::Preconditioner { LUSolver Lpre; //! \brief Whether or not to use a symmetric preconditioner - bool symmetric; + bool symmetric; }; } diff --git a/opm/elasticity/mpc.hh b/opm/elasticity/mpc.hh index 30a839e..4f64009 100644 --- a/opm/elasticity/mpc.hh +++ b/opm/elasticity/mpc.hh @@ -22,7 +22,7 @@ namespace Elasticity { //! \brief An enum for specification of global coordinate directions enum Direction { NONE = 0, X = 1, Y = 2, Z = 4, - XY = 1+2, XZ = 1+4, YZ = 2+4, + XY = 1+2, XZ = 1+4, YZ = 2+4, XYZ = 1+2+4 }; /*! diff --git a/opm/porsol/common/GridInterfaceEuler.hpp b/opm/porsol/common/GridInterfaceEuler.hpp index a44efa6..5c83b4a 100644 --- a/opm/porsol/common/GridInterfaceEuler.hpp +++ b/opm/porsol/common/GridInterfaceEuler.hpp @@ -109,7 +109,7 @@ namespace Opm : pgrid_(&grid), iter_(it), local_index_(loc_ind) { } - + /// @brief /// @todo Doc me! /// @return diff --git a/opm/porsol/common/PeriodicHelpers.hpp b/opm/porsol/common/PeriodicHelpers.hpp index d9e921c..bfc5bc4 100644 --- a/opm/porsol/common/PeriodicHelpers.hpp +++ b/opm/porsol/common/PeriodicHelpers.hpp @@ -16,7 +16,7 @@ /* Copyright 2009, 2010 SINTEF ICT, Applied Mathematics. Copyright 2009, 2010 Statoil ASA. - + This file is part of The Open Reservoir Simulator Project (OpenRS). OpenRS is free software: you can redistribute it and/or modify @@ -258,7 +258,7 @@ namespace Opm double cmin = low[pddir]; double cmax = hi[pddir]; double cdelta = cmax - cmin; - + for (int i = 0; i < num_bdy; ++i) { Vector fcent = bface_iters[i]->centroid(); int canon_pos = -1; diff --git a/opm/porsol/common/ReservoirPropertyCommon.hpp b/opm/porsol/common/ReservoirPropertyCommon.hpp index 55d6af5..16748db 100644 --- a/opm/porsol/common/ReservoirPropertyCommon.hpp +++ b/opm/porsol/common/ReservoirPropertyCommon.hpp @@ -65,7 +65,7 @@ namespace Opm /// @brief Tensor type for read and write access to permeability. typedef SharedCMatrix SharedPermTensor; - /// @brief The number of phases + /// @brief The number of phases enum { NumberOfPhases = 2 }; /// @brief Default constructor. diff --git a/opm/porsol/euler/EulerUpstreamResidual_impl.hpp b/opm/porsol/euler/EulerUpstreamResidual_impl.hpp index 90485d3..3e23ec8 100644 --- a/opm/porsol/euler/EulerUpstreamResidual_impl.hpp +++ b/opm/porsol/euler/EulerUpstreamResidual_impl.hpp @@ -192,7 +192,7 @@ namespace Opm // Computing G. Note that we do not multiply with the mobility, // so this G is wrong in case of anisotropic relperm. const double G = s.method_gravity_ ? - loc_area*inner(loc_normal, grav_influence) + loc_area*inner(loc_normal, grav_influence) : 0.0; const int triv_phase = G >= 0.0 ? 0 : 1; const int ups_cell = loc_flux >= 0.0 ? 0 : 1; diff --git a/opm/porsol/euler/EulerUpstream_impl.hpp b/opm/porsol/euler/EulerUpstream_impl.hpp index e0bfd48..e3e24e4 100644 --- a/opm/porsol/euler/EulerUpstream_impl.hpp +++ b/opm/porsol/euler/EulerUpstream_impl.hpp @@ -354,7 +354,7 @@ namespace Opm - + template template inline void EulerUpstream::smallTimeStep(std::vector& saturation, diff --git a/opm/porsol/euler/ImplicitCapillarity.hpp b/opm/porsol/euler/ImplicitCapillarity.hpp index 97e0803..f1c5a9a 100644 --- a/opm/porsol/euler/ImplicitCapillarity.hpp +++ b/opm/porsol/euler/ImplicitCapillarity.hpp @@ -53,7 +53,7 @@ namespace Opm { { public: typedef IncompFlowSolverHybrid PressureSolver; diff --git a/opm/porsol/mimetic/IncompFlowSolverHybrid.hpp b/opm/porsol/mimetic/IncompFlowSolverHybrid.hpp index 5ebf964..aa02d9c 100644 --- a/opm/porsol/mimetic/IncompFlowSolverHybrid.hpp +++ b/opm/porsol/mimetic/IncompFlowSolverHybrid.hpp @@ -435,7 +435,7 @@ namespace Opm { Opm::SparseTable< int > cellFaces_; std::vector pressure_; Opm::SparseTable outflux_; - + void clear() { std::vector().swap(cellno_); cellFaces_.clear(); @@ -676,16 +676,16 @@ namespace Opm { solveLinearSystem(residual_tolerance, linsolver_verbosity, linsolver_maxit); break; case 1: // AMG preconditioned CG - solveLinearSystemAMG(residual_tolerance, linsolver_verbosity, + solveLinearSystemAMG(residual_tolerance, linsolver_verbosity, linsolver_maxit, prolongate_factor, same_matrix, smooth_steps); break; - + case 2: // KAMG preconditioned CG - solveLinearSystemKAMG(residual_tolerance, linsolver_verbosity, + solveLinearSystemKAMG(residual_tolerance, linsolver_verbosity, linsolver_maxit, prolongate_factor, same_matrix,smooth_steps); break; case 3: // CG preconditioned with AMG that uses less memory badwidth - solveLinearSystemFastAMG(residual_tolerance, linsolver_verbosity, + solveLinearSystemFastAMG(residual_tolerance, linsolver_verbosity, linsolver_maxit, prolongate_factor, same_matrix,smooth_steps); break; default: @@ -1485,7 +1485,7 @@ namespace Opm { S_[0][0] *= 2; } opS_.reset(new Operator(S_)); - + // Construct preconditioner. double relax = 1; typename Precond::SmootherArgs smootherArgs; @@ -1511,7 +1511,7 @@ namespace Opm { Dune::InverseOperatorResult result; soln_ = 0.0; - // Adapt initial guess such Dirichlet boundary conditions are + // Adapt initial guess such Dirichlet boundary conditions are // represented, i.e. soln_i=A_{ii}^-1 rhs_i typedef typename Dune::BCRSMatrix ::ConstRowIterator RowIter; typedef typename Dune::BCRSMatrix ::ConstColIterator ColIter; @@ -1553,7 +1553,7 @@ namespace Opm { S_[0][0] *= 2; } opS_.reset(new Operator(S_)); - + // Construct preconditioner. typedef Dune::Amg::AggregationCriterion > CritBase; @@ -1577,7 +1577,7 @@ namespace Opm { Dune::InverseOperatorResult result; soln_ = 0.0; - // Adapt initial guess such Dirichlet boundary conditions are + // Adapt initial guess such Dirichlet boundary conditions are // represented, i.e. soln_i=A_{ii}^-1 rhs_i typedef typename Dune::BCRSMatrix ::ConstRowIterator RowIter; typedef typename Dune::BCRSMatrix ::ConstColIterator ColIter; @@ -1606,8 +1606,8 @@ namespace Opm { void solveLinearSystemKAMG(double residual_tolerance, int verbosity_level, int maxit, double prolong_factor, bool same_matrix, int smooth_steps) // ---------------------------------------------------------------- - { - + { + typedef Dune::Amg::KAMG Precond; // Adapted from upscaling.cc by Arne Rekdal, 2009 Scalar residTol = residual_tolerance; @@ -1617,7 +1617,7 @@ namespace Opm { S_[0][0] *= 2; } opS_.reset(new Operator(S_)); - + // Construct preconditioner. double relax = 1; typename Precond::SmootherArgs smootherArgs; @@ -1643,7 +1643,7 @@ namespace Opm { Dune::InverseOperatorResult result; soln_ = 0.0; - // Adapt initial guess such Dirichlet boundary conditions are + // Adapt initial guess such Dirichlet boundary conditions are // represented, i.e. soln_i=A_{ii}^-1 rhs_i typedef typename Dune::BCRSMatrix ::ConstRowIterator RowIter; typedef typename Dune::BCRSMatrix ::ConstColIterator ColIter; diff --git a/opm/upscaling/SteadyStateUpscalerImplicit_impl.hpp b/opm/upscaling/SteadyStateUpscalerImplicit_impl.hpp index 473aeae..95ca405 100644 --- a/opm/upscaling/SteadyStateUpscalerImplicit_impl.hpp +++ b/opm/upscaling/SteadyStateUpscalerImplicit_impl.hpp @@ -190,7 +190,7 @@ namespace Opm if (use_gravity_) { gravity[2] = Opm::unit::gravity; } - + if (gravity.two_norm() > 0.0) { OPM_MESSAGE("Warning: Gravity is experimental for flow solver."); } diff --git a/opm/upscaling/SteadyStateUpscalerManager.hpp b/opm/upscaling/SteadyStateUpscalerManager.hpp index 8bed443..ced6de5 100644 --- a/opm/upscaling/SteadyStateUpscalerManager.hpp +++ b/opm/upscaling/SteadyStateUpscalerManager.hpp @@ -101,8 +101,8 @@ namespace Opm { /* const int num_rows = K.numRows(); */ /* const int num_cols = K.numCols(); */ - - /* We write tensor output in Voigt notation, + + /* We write tensor output in Voigt notation, (but we also output the remainding three terms) http://en.wikipedia.org/wiki/Voigt_notation */ @@ -113,7 +113,7 @@ namespace Opm If linear or periodic bc's, output all 9 elements (even though 6 is strictly necessary for periodic bc's) Use the Tensor class to order output into Voigt notation, so that it works for more than 3x3 matrices */ - + os << pdrop << '\t'; os << sat << '\t'; os << K(0,0) << '\t'; /* xx */ @@ -125,7 +125,7 @@ namespace Opm os << K(2,1) << '\t'; /* zy */ os << K(2,0) << '\t'; /* zx */ os << K(1,2); /* yz */ - + os << std::endl; @@ -231,10 +231,10 @@ namespace Opm // Changing initial saturations for next pressure drop to equal the steady state of the last init_sat = upscaler.lastSaturationState(); - + writeRelPerm(krw_out, lambda.first , usat, pdrop); writeRelPerm(kro_out, lambda.second, usat, pdrop); - + } } } diff --git a/opm/upscaling/SteadyStateUpscalerManagerImplicit.hpp b/opm/upscaling/SteadyStateUpscalerManagerImplicit.hpp index 34a02ad..4f37509 100644 --- a/opm/upscaling/SteadyStateUpscalerManagerImplicit.hpp +++ b/opm/upscaling/SteadyStateUpscalerManagerImplicit.hpp @@ -189,7 +189,7 @@ namespace Opm std::string kro_filename = param.getDefault("outputOil", krotmpname); std::ofstream krw_out(krw_filename.c_str()); std::ofstream kro_out(kro_filename.c_str()); - // std::stringstream krw_out; + // std::stringstream krw_out; // std::stringstream kro_out; krw_out << "# Result from steady state upscaling" << std::endl; krw_out << "# Pressuredrop Sw Krxx Kryy Krzz" << std::endl; diff --git a/opm/upscaling/SteadyStateUpscaler_impl.hpp b/opm/upscaling/SteadyStateUpscaler_impl.hpp index 10d51c1..dcfd03b 100644 --- a/opm/upscaling/SteadyStateUpscaler_impl.hpp +++ b/opm/upscaling/SteadyStateUpscaler_impl.hpp @@ -70,7 +70,7 @@ namespace Opm inline void SteadyStateUpscaler::initImpl(const Opm::parameter::ParameterGroup& param) { Super::initImpl(param); - use_gravity_ = param.getDefault("use_gravity", use_gravity_); + use_gravity_ = param.getDefault("use_gravity", use_gravity_); output_vtk_ = param.getDefault("output_vtk", output_vtk_); print_inoutflows_ = param.getDefault("print_inoutflows", print_inoutflows_); simulation_steps_ = param.getDefault("simulation_steps", simulation_steps_); @@ -165,7 +165,7 @@ namespace Opm // Run pressure solver. this->flow_solver_.solve(this->res_prop_, saturation, this->bcond_, src, - this->residual_tolerance_, this->linsolver_verbosity_, + this->residual_tolerance_, this->linsolver_verbosity_, this->linsolver_type_, false, this->linsolver_maxit_, this->linsolver_prolongate_factor_, this->linsolver_smooth_steps_); @@ -271,7 +271,7 @@ namespace Opm // Compute the (anisotropic) upscaled mobilities. // eff_Kw := lambda_w*K - // => lambda_w = eff_Kw*inv(K); + // => lambda_w = eff_Kw*inv(K); permtensor_t lambda_w(matprod(eff_Kw, inverse3x3(upscaled_perm))); permtensor_t lambda_o(matprod(eff_Ko, inverse3x3(upscaled_perm))); diff --git a/opm/upscaling/UpscalerBase_impl.hpp b/opm/upscaling/UpscalerBase_impl.hpp index 6fb6f0d..3ca10e5 100644 --- a/opm/upscaling/UpscalerBase_impl.hpp +++ b/opm/upscaling/UpscalerBase_impl.hpp @@ -255,7 +255,7 @@ namespace Opm // Run pressure solver. bool same_matrix = (bctype_ != Fixed) && (pdd != 0); flow_solver_.solve(fluid, sat, bcond_, src, residual_tolerance_, - linsolver_verbosity_, + linsolver_verbosity_, linsolver_type_, same_matrix, linsolver_maxit_, linsolver_prolongate_factor_, linsolver_smooth_steps_); @@ -330,7 +330,7 @@ namespace Opm side2_flux += flux*norm_comp; side2_area += area; } - } + } } } } @@ -362,7 +362,7 @@ namespace Opm side2_pos += area*pos_comp; side2_area += area; } - } + } } } } diff --git a/opm/upscaling/initCPGrid.cpp b/opm/upscaling/initCPGrid.cpp index cdc4045..8b9033e 100644 --- a/opm/upscaling/initCPGrid.cpp +++ b/opm/upscaling/initCPGrid.cpp @@ -74,4 +74,3 @@ void Opm::initCPGrid(Dune::CpGrid& grid, const Opm::ParameterGroup& param) { OPM_THROW(std::runtime_error, "Unknown file format string: " + fileformat); } } - diff --git a/tests/common/test_boundaryconditions.cpp b/tests/common/test_boundaryconditions.cpp index a1c2b96..b1d80df 100644 --- a/tests/common/test_boundaryconditions.cpp +++ b/tests/common/test_boundaryconditions.cpp @@ -133,4 +133,3 @@ BOOST_AUTO_TEST_CASE(surfvolboundaryconditions) BOOST_CHECK(bc2.surfvolCond(1).isDirichlet()); BOOST_CHECK(bc2.surfvolCond(1).surfvol() == z); } -