| Line | Branch | Exec | Source |
|---|---|---|---|
| 1 | // -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- | ||
| 2 | // vi: set et ts=4 sw=4 sts=4: | ||
| 3 | // | ||
| 4 | // SPDX-FileCopyrightInfo: Copyright © DuMux Project contributors, see AUTHORS.md in root folder | ||
| 5 | // SPDX-License-Identifier: GPL-3.0-or-later | ||
| 6 | // | ||
| 7 | /*! | ||
| 8 | * \file | ||
| 9 | * \ingroup Material | ||
| 10 | * \brief Creating random fields using gstat | ||
| 11 | */ | ||
| 12 | #ifndef DUMUX_MATERIAL_GSTAT_RANDOM_FIELD_HH | ||
| 13 | #define DUMUX_MATERIAL_GSTAT_RANDOM_FIELD_HH | ||
| 14 | |||
| 15 | #include <iostream> | ||
| 16 | #include <algorithm> | ||
| 17 | #include <string> | ||
| 18 | #include <cmath> | ||
| 19 | |||
| 20 | #include <dune/common/exceptions.hh> | ||
| 21 | #include <dune/grid/common/mcmgmapper.hh> | ||
| 22 | #include <dune/grid/io/file/vtk.hh> | ||
| 23 | |||
| 24 | namespace Dumux { | ||
| 25 | |||
| 26 | /*! | ||
| 27 | * \ingroup Material | ||
| 28 | * \brief Creating random fields using gstat | ||
| 29 | * | ||
| 30 | * gstat is an open source software tool which can (among other things) generate | ||
| 31 | * geostatistical random fields (see <a href="http://www.gstat.org">http://www.gstat.org</a> | ||
| 32 | * or \cite Pebesma1998a). | ||
| 33 | * | ||
| 34 | * To use this class, execute the installexternal.py from your DuMuX root | ||
| 35 | * directory or download, unpack and install the tarball from the gstat-website. | ||
| 36 | * Then rerun cmake (in the second case set GSTAT_ROOT in your input file to the | ||
| 37 | * path where gstat is installed). | ||
| 38 | */ | ||
| 39 | template<class GridView, class Scalar> | ||
| 40 | ✗ | class GstatRandomField | |
| 41 | { | ||
| 42 | enum { dim = GridView::dimension }; | ||
| 43 | enum { dimWorld = GridView::dimensionworld }; | ||
| 44 | |||
| 45 | using DataVector = std::vector<Scalar>; | ||
| 46 | using Element = typename GridView::Traits::template Codim<0>::Entity; | ||
| 47 | using ElementMapper = Dune::MultipleCodimMultipleGeomTypeMapper<GridView>; | ||
| 48 | |||
| 49 | public: | ||
| 50 | // Add field types if you want to implement e.g. tensor permeabilities. | ||
| 51 | enum FieldType { scalar, log10 }; | ||
| 52 | |||
| 53 | /*! | ||
| 54 | * \brief Constructor | ||
| 55 | * | ||
| 56 | * \param gridView the used gridView | ||
| 57 | * \param elementMapper Maps elements of the given grid view | ||
| 58 | */ | ||
| 59 | ✗ | GstatRandomField(const GridView& gridView, const ElementMapper& elementMapper) | |
| 60 | : gridView_(gridView) | ||
| 61 | , elementMapper_(elementMapper) | ||
| 62 | ✗ | , data_(gridView.size(0)) | |
| 63 | ✗ | {} | |
| 64 | |||
| 65 | /*! | ||
| 66 | * \brief Creates a new field with random variables, if desired. | ||
| 67 | * Otherwise creates a data field from already available data. | ||
| 68 | * For the random field generation three files are necessary. | ||
| 69 | * | ||
| 70 | * A \a gstatControlFile in which all commands and in/output files for gstat are specified. | ||
| 71 | * A \a gstatInputFile contains all coordinates (cell centers) of the grid, so that | ||
| 72 | * gstat can perform its random realization. The filename must be same as in the gstatControlFile. | ||
| 73 | * A \a gstatOutputFile in which gstat writes the random values to this file. | ||
| 74 | * The filename must be the same as in the gstatControlFile. | ||
| 75 | * \param fieldType | ||
| 76 | * \param gstatControlFile name of control file for gstat | ||
| 77 | * \param gstatInputFile name of input file for gstat | ||
| 78 | * \param gstatOutputFile name of the gstat output file | ||
| 79 | * \param createNew set true to create a new field | ||
| 80 | */ | ||
| 81 | ✗ | void create(const std::string& gstatControlFile, | |
| 82 | const std::string& gstatInputFile = "gstatInput.txt", | ||
| 83 | const std::string& gstatOutputFile = "permeab.dat", | ||
| 84 | FieldType fieldType = FieldType::scalar, | ||
| 85 | bool createNew = true) | ||
| 86 | { | ||
| 87 | ✗ | fieldType_ = fieldType; | |
| 88 | ✗ | if (createNew) | |
| 89 | { | ||
| 90 | #if !DUMUX_HAVE_GSTAT | ||
| 91 | ✗ | DUNE_THROW(Dune::InvalidStateException, "Requested data field generation with gstat" | |
| 92 | << " but gstat was not found on your system. Set GSTAT_ROOT to the path where gstat " | ||
| 93 | << " is installed and pass it to CMake, e.g. through an opts file."); | ||
| 94 | #else | ||
| 95 | std::ofstream gstatInput(gstatInputFile); | ||
| 96 | for (const auto& element : elements(gridView_)) | ||
| 97 | // get global coordinates of cell centers | ||
| 98 | gstatInput << element.geometry().center() << std::endl; | ||
| 99 | |||
| 100 | gstatInput.close(); | ||
| 101 | |||
| 102 | std::string syscom; | ||
| 103 | syscom = GSTAT_EXECUTABLE; | ||
| 104 | syscom += " "; | ||
| 105 | syscom += gstatControlFile; | ||
| 106 | |||
| 107 | if (!gstatInput.good()) | ||
| 108 | { | ||
| 109 | DUNE_THROW(Dune::IOError, "Reading the gstat control file: " | ||
| 110 | << gstatControlFile << " failed." << std::endl); | ||
| 111 | } | ||
| 112 | |||
| 113 | if (system(syscom.c_str())) | ||
| 114 | { | ||
| 115 | DUNE_THROW(Dune::IOError, "Executing gstat failed."); | ||
| 116 | } | ||
| 117 | #endif | ||
| 118 | } | ||
| 119 | |||
| 120 | ✗ | std::ifstream gstatOutput(gstatOutputFile); | |
| 121 | ✗ | if (!gstatOutput.good()) | |
| 122 | { | ||
| 123 | ✗ | DUNE_THROW(Dune::IOError, "Reading from file: " | |
| 124 | << gstatOutputFile << " failed." << std::endl); | ||
| 125 | } | ||
| 126 | |||
| 127 | ✗ | std::string line; | |
| 128 | ✗ | std::getline(gstatOutput, line); | |
| 129 | |||
| 130 | Scalar trash, dataValue; | ||
| 131 | ✗ | for (const auto& element : elements(gridView_)) | |
| 132 | { | ||
| 133 | ✗ | std::getline(gstatOutput, line); | |
| 134 | ✗ | std::istringstream curLine(line); | |
| 135 | if (dim == 1) | ||
| 136 | curLine >> trash >> dataValue; | ||
| 137 | else if (dim == 2) | ||
| 138 | ✗ | curLine >> trash >> trash >> dataValue; | |
| 139 | else if (dim == 3) | ||
| 140 | curLine >> trash >> trash >> trash >> dataValue; | ||
| 141 | else | ||
| 142 | DUNE_THROW(Dune::InvalidStateException, "Invalid dimension " << dim); | ||
| 143 | |||
| 144 | ✗ | data_[elementMapper_.index(element)] = dataValue; | |
| 145 | } | ||
| 146 | ✗ | gstatOutput.close(); | |
| 147 | |||
| 148 | // post processing | ||
| 149 | using std::pow; | ||
| 150 | ✗ | if (fieldType_ == FieldType::log10) | |
| 151 | ✗ | std::for_each(data_.begin(), data_.end(), [](Scalar& s){ s = pow(10.0, s); }); | |
| 152 | ✗ | } | |
| 153 | |||
| 154 | //! \brief Return an entry of the data vector | ||
| 155 | Scalar data(const Element& e) const | ||
| 156 | { | ||
| 157 | return data_[elementMapper_.index(e)]; | ||
| 158 | } | ||
| 159 | |||
| 160 | //! \brief Return the data vector for analysis or external vtk output | ||
| 161 | const DataVector& data() const | ||
| 162 | { | ||
| 163 | return data_; | ||
| 164 | } | ||
| 165 | |||
| 166 | //! \brief Write the data to a vtk file | ||
| 167 | void writeVtk(const std::string& vtkName, | ||
| 168 | const std::string& dataName = "data") const | ||
| 169 | { | ||
| 170 | Dune::VTKWriter<GridView> vtkwriter(gridView_); | ||
| 171 | vtkwriter.addCellData(data_, dataName); | ||
| 172 | |||
| 173 | DataVector logPerm; | ||
| 174 | if (fieldType_ == FieldType::log10) | ||
| 175 | { | ||
| 176 | logPerm = data_; | ||
| 177 | using std::log10; | ||
| 178 | std::for_each(logPerm.begin(), logPerm.end(), [](Scalar& s){ s = log10(s); }); | ||
| 179 | vtkwriter.addCellData(logPerm, "log10 of " + dataName); | ||
| 180 | } | ||
| 181 | vtkwriter.write(vtkName, Dune::VTK::OutputType::ascii); | ||
| 182 | } | ||
| 183 | |||
| 184 | private: | ||
| 185 | const GridView gridView_; | ||
| 186 | const ElementMapper& elementMapper_; | ||
| 187 | DataVector data_; | ||
| 188 | FieldType fieldType_; | ||
| 189 | }; | ||
| 190 | |||
| 191 | } // end namespace Dumux | ||
| 192 | |||
| 193 | #endif | ||
| 194 |