From c202ce450719dab1b16f0ad753b54d7821d44347 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?H=C3=A5vard=20Berland?= Date: Sun, 7 Nov 2010 21:08:34 +0100 Subject: [PATCH] first version of code for experimental variogram --- examples/Makefile.am | 5 +- examples/exp_variogram.cpp | 243 +++++++++++++++++++++++++++++++++++++ 2 files changed, 247 insertions(+), 1 deletion(-) create mode 100644 examples/exp_variogram.cpp diff --git a/examples/Makefile.am b/examples/Makefile.am index 131ba4b..5515aea 100644 --- a/examples/Makefile.am +++ b/examples/Makefile.am @@ -6,7 +6,8 @@ noinst_PROGRAMS = upscale_perm \ upscale_relperm \ upscale_avg \ upscale_cond \ - cpchop + cpchop \ + exp_variogram upscale_perm_SOURCES = upscale_perm.cpp @@ -18,6 +19,8 @@ upscale_avg_SOURCES = upscale_avg.cpp cpchop_SOURCES = cpchop.cpp +exp_variogram_SOURCES = exp_variogram.cpp + AM_CPPFLAGS += $(DUNEMPICPPFLAGS) $(BOOST_CPPFLAGS) $(SUPERLU_CPPFLAGS) AM_LDFLAGS += $(DUNEMPILDFLAGS) $(BOOST_LDFLAGS) $(SUPERLU_LDFLAGS) LDADD = $(DUNEMPILIBS) $(BOOST_UNIT_TEST_FRAMEWORK_LIB) \ diff --git a/examples/exp_variogram.cpp b/examples/exp_variogram.cpp new file mode 100644 index 0000000..aed3ba1 --- /dev/null +++ b/examples/exp_variogram.cpp @@ -0,0 +1,243 @@ +/* + Copyright 2010 SINTEF ICT, Applied Mathematics. + + This file is part of the Open Porous Media project (OPM). + + OPM is free software: you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + OPM is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with OPM. If not, see . +*/ + +/* + This program computes data for an experimental + variogram from a cornerpoint geometry with properties. + + This works by choosing pairs of volumes chosen + randomly, and comparing their distance with the + difference in their properties (poro or perm) + + Direction for pairing is set from the command line + +*/ + +#include +#include +#include +#include + +#include +#include +#include +#include + +#include +#include +#include +#include +#include +#include +#include +#include + + +int main(int argc, char** argv) +{ + Dune::parameter::ParameterGroup param(argc, argv); + std::string gridfilename = param.get("gridfilename"); + Dune::CornerPointChopper ch(gridfilename); + + // The cells with i coordinate in [imin, imax) are included, similar for j. + // The z limits may be changed inside the chopper to match actual min/max z. + const int* dims = ch.dimensions(); + int imin = param.getDefault("imin", 0); + int imax = param.getDefault("imax", dims[0]); + int jmin = param.getDefault("jmin", 0); + int jmax = param.getDefault("jmax", dims[1]); + double zmin = param.getDefault("zmin", ch.zLimits().first); + double zmax = param.getDefault("zmax", ch.zLimits().second); + int pairs = param.getDefault("pairs", 100); + std::string direction = param.get("direction"); + int ilen = param.getDefault("ilen", 0); + int jlen = param.getDefault("jlen", 0); + double zlen = param.getDefault("zlen", 0.0); + boost::mt19937::result_type userseed = param.getDefault("seed", 0); + + int outputprecision = param.getDefault("outputprecision", 8); + std::string resultfile = param.getDefault("resultfile", ""); + + double z_tolerance = param.getDefault("z_tolerance", 0.0); + 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); + + if (ilen <= 0) { + std::cerr << "Error: ilen (" << ilen << ") must be greater than zero\n"; + exit(1); + } + if (jlen <= 0) { + std::cerr << "Error: jlen (" << jlen <<") must be greater than zero\n"; + exit(1); + } + if (zlen <= 0.0) { + std::cerr << "Eror: zlen (" << zlen <<") must be greater than zero\n"; + exit(1); + } + + // Check user supplied variogram direction, either horizontal or vertical + enum variogram_directions { undefined, horizontal, vertical }; + variogram_directions variogram_direction = undefined; + std::string distancemetric; + if ( direction == "vertical" || direction == "vert" ) { + variogram_direction = vertical; + distancemetric = "zcorn"; + } + else { // if (direction == "horizontal" || direction == "horiz") { + variogram_direction = horizontal; + distancemetric = "cells"; + } + + + if (variogram_direction == undefined) { + std::cerr << "Error: variogram direction is undefined, user supplied '" << direction << "'.\n"; + exit(1); + } + + + // 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); + + // Random number generator from boost. + boost::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) { + gen.seed(time(NULL)); + } + else { + gen.seed(userseed); + } + + + // Note that end is included in interval for uniform_int. + boost::uniform_int<> disti(imin, imax - ilen); + boost::uniform_int<> distj(jmin, jmax - jlen); + boost::uniform_real<> distz(zmin, std::max(zmax - zlen, zmin)); + boost::variate_generator > ri(gen, disti); + boost::variate_generator > rj(gen, distj); + boost::variate_generator > rz(gen, distz); + + // 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(); + double zstart_1 = rz(); + ch.chop(istart_1, istart_1 + ilen, jstart_1, jstart_1 + jlen, zstart_1, zstart_1 + zlen, false); + + Dune::EclipseGridParser subparser_1 = ch.subparser(); + + Dune::SinglePhaseUpscaler upscaler_1; + upscaler_1.init(subparser_1, Dune::SinglePhaseUpscaler::Fixed, 0.0, z_tolerance, + residual_tolerance, linsolver_verbosity, linsolver_type, false); + Dune::SinglePhaseUpscaler::permtensor_t upscaled_K_1 = upscaler_1.upscaleSinglePhase(); + upscaled_K_1 *= (1.0/(Dune::prefix::milli*Dune::unit::darcy)); + double porosity_1 = upscaler_1.upscalePorosity(); + + // Pick another location to form a location-pair to be compared + int istart_2, jstart_2, zstart_2; + if (variogram_direction == horizontal) { + istart_2 = ri(); + jstart_2 = rj(); + zstart_2 = zstart_1; + } + else if (variogram_direction == vertical) { + 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); + + Dune::EclipseGridParser subparser_2 = ch.subparser(); + + Dune::SinglePhaseUpscaler upscaler_2; + upscaler_2.init(subparser_2, Dune::SinglePhaseUpscaler::Fixed, 0.0, z_tolerance, + residual_tolerance, linsolver_verbosity, linsolver_type, false); + Dune::SinglePhaseUpscaler::permtensor_t upscaled_K_2 = upscaler_2.upscaleSinglePhase(); + upscaled_K_2 *= (1.0/(Dune::prefix::milli*Dune::unit::darcy)); + double porosity_2 = upscaler_2.upscalePorosity(); + + if (variogram_direction == horizontal) { + distances.push_back(sqrt(pow(istart_2 - istart_1,2) + pow(jstart_2 - jstart_1,2))); + } + else if (variogram_direction == vertical) { + distances.push_back(fabs(zstart_2 - zstart_1)); + } + porodiffs.push_back(fabs(porosity_2 - porosity_1)); + permxdiffs.push_back(fabs(upscaled_K_2(0,0) - upscaled_K_1(0,0))); + 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; + outputtmp << "# Options used:" << std::endl; + outputtmp << "# gridfilename: " << gridfilename << std::endl; + outputtmp << "# i; min,len,max: " << imin << " " << ilen << " " << imax << std::endl; + outputtmp << "# j; min,len,max: " << jmin << " " << jlen << " " << jmax << std::endl; + outputtmp << "# z; min,len,max: " << zmin << " " << zlen << " " << zmax << std::endl; + outputtmp << "# direction: " << direction << std::endl; + 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' << + std::showpoint << std::setw(fieldwidth) << std::setprecision(outputprecision) << porodiffs[pair-1] << '\t' << + std::showpoint << std::setw(fieldwidth) << std::setprecision(outputprecision) << permxdiffs[pair-1] << '\t' << + std::showpoint << std::setw(fieldwidth) << std::setprecision(outputprecision) << permydiffs[pair-1] << '\t' << + 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(); + } + + + + std::cout << outputtmp.str(); +}