diff --git a/CMakeLists.txt b/CMakeLists.txt index 2d6d97067d523c2b94ef40c96c9ccdb850f33c4c..b7268ed947d80db257899795a5b73ec3a59727e7 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -12,8 +12,19 @@ if(NOT (dune-common_DIR ${PROJECT_BINARY_DIR}) endif() +if(NOT (opm-common_DIR + OR opm-common_ROOT + OR "${CMAKE_PREFIX_PATH}" MATCHES ".*opm-common.*")) + string(REPLACE ${CMAKE_PROJECT_NAME} + opm-common opm-common_DIR + ${PROJECT_BINARY_DIR}) +endif() + #find dune-common and set the module path find_package(dune-common) +#find opm-common to support OPM +find_package(opm-common) + list(APPEND CMAKE_MODULE_PATH ${dune-common_MODULE_PATH} "${PROJECT_SOURCE_DIR}/cmake/modules") #include the dune macros diff --git a/bin/installexternal.sh b/bin/installexternal.sh index 20cce997606cc11595a21deffe2bd54b2558a46f..561da744c26850b8dd5f1ad38a3090988afbf21c 100755 --- a/bin/installexternal.sh +++ b/bin/installexternal.sh @@ -53,17 +53,12 @@ installAluGrid() fi } -installErt() +installEcl() { cd $TOPDIR - checkLocationForDuneModules ert - if test $CORRECT_LOCATION_FOR_DUNE_MODULES == "n"; then - return - fi - - if [ ! -e ert ]; then - git clone -b release/2017.04 https://github.com/Ensembles/ert.git + if [ ! -e libecl ]; then + git clone -b 2018.04 https://github.com/Statoil/libecl.git fi if test "$DOWNLOAD_ONLY" == "y"; then @@ -71,21 +66,20 @@ installErt() fi if test "$CLEANUP" == "y"; then - rm -rf ert + rm -rf libecl return fi - # building ert - echo "Building ert" - cd $TOPDIR/ert + # building libecl + echo "Building libecl" + cd $TOPDIR/libecl mkdir build cd build cmake .. make # show additional information - echo "Ert has been built in directory ert/build." - echo "Do not change this directory otherwise opm will not find ert!" + echo "Ecl has been built in directory libecl/build." cd $TOPDIR } @@ -294,72 +288,38 @@ installOPM() fi if [ ! -e opm-common ]; then - git clone -b release/2017.04 https://github.com/OPM/opm-common - fi - - if [ ! -e opm-core ]; then - git clone -b release/2017.04 https://github.com/OPM/opm-core - fi - - if [ ! -e opm-material ]; then - git clone -b release/2017.04 https://github.com/OPM/opm-material - fi - - if [ ! -e opm-parser ]; then - git clone -b release/2017.04 https://github.com/OPM/opm-parser + git clone -b release/2018.04 https://github.com/OPM/opm-common fi if [ ! -e opm-grid ]; then - git clone -b release/2017.04 https://github.com/OPM/opm-grid + git clone -b release/2018.04 https://github.com/OPM/opm-grid fi - if [ ! -e opm-output ]; then - git clone -b release/2017.04 https://github.com/OPM/opm-output - fi - if test "$DOWNLOAD_ONLY" == "y"; then return fi if test "$CLEANUP" == "y"; then rm -rf opm-common - rm -rf opm-core - rm -rf opm-material - rm -rf opm-parser rm -rf opm-grid - rm -rf opm-output return fi # apply patches echo "Applying patch for opm-common" cd $TOPDIR/opm-common - patch -p1 < $TOPDIR/dumux/patches/opm-common-2017.04.patch - - echo "Applying patch for opm-core" - cd $TOPDIR/opm-core - patch -p1 < $TOPDIR/dumux/patches/opm-core-2017.04.patch - - echo "Applying patch for opm-parser" - cd $TOPDIR/opm-parser - patch -p1 < $TOPDIR/dumux/patches/opm-parser-2017.04.patch + patch -p1 < $TOPDIR/dumux/patches/opm-common-2018.04.patch echo "Applying patch for opm-grid" cd $TOPDIR/opm-grid - patch -p1 < $TOPDIR/dumux/patches/opm-grid-2017.04.patch + patch -p1 < $TOPDIR/dumux/patches/opm-grid-2018.04.patch # show additional information echo "In addition, it might be necessary to set manually some" echo "CMake variables in the CMAKE_FLAGS section of the .opts-file:" echo " -DOPM_COMMON_ROOT=/path/to/opm-common \\" - echo " -Dopm-grid_PREFIX=/path/to/opm-grid \\" - echo " -Dopm-common_PREFIX=/path/to/opm-common \\" - echo " -Dopm-core_PREFIX=/path/to/opm-core \\" - echo " -Dopm-material_PREFIX=/path/to/opm-material \\" - echo " -Dopm-parser_PREFIX=/path/to/opm-parser \\" - echo " -Dopm-output_PREFIX=/path/to/opm-output \\" + echo " -Decl_DIR=/path/to/libecl/build \\" echo " -DUSE_MPI=ON \\" - echo " -DHAVE_OPM_GRID=1 \\" cd $TOPDIR } @@ -476,7 +436,7 @@ usage() echo "Where PACKAGES is one or more of the following" echo " all Install everything and the kitchen sink." echo " alugrid Download dune-alugrid." - echo " ert Download and build ert." + echo " ecl Download and build ecl." echo " foamgrid Download dune-foamgrid." echo " glpk Download and install glpk." echo " gstat Download and install gstat." @@ -544,7 +504,7 @@ for TMP in "$@"; do SOMETHING_DONE="y" createExternalDirectory installAluGrid - installErt + installEcl installFoamGrid installGLPK installGStat @@ -561,9 +521,9 @@ for TMP in "$@"; do SOMETHING_DONE="y" installAluGrid ;; - ert) + ecl) SOMETHING_DONE="y" - installErt + installEcl ;; foamgrid|dune-foamgrid) SOMETHING_DONE="y" diff --git a/dumux/discretization/cellcentered/subcontrolvolume.hh b/dumux/discretization/cellcentered/subcontrolvolume.hh index ae1f51a22549592b3d9b70d34f785885d3441324..7bb3627b2821a764a18d55f6d1f63ac506634166 100644 --- a/dumux/discretization/cellcentered/subcontrolvolume.hh +++ b/dumux/discretization/cellcentered/subcontrolvolume.hh @@ -65,6 +65,18 @@ class CCSubControlVolume using LocalIndexType = typename T::LocalIndexType; using Scalar = typename T::Scalar; + // In the following, the correct parameter type for the geometry passed to + // the constructor below is determined. It depends upon whether the + // `geometry()` method of the `Element` returns a copy or a reference. + // In the first case, the correct type is `Geometry&&`, while it is + // `Geometry` for the second case. Although returning by copy is prescribed + // by the Dune interface, the grid implementation CpGrid uses a const + // reference as of Opm 2018.04. Once this is fixed, the parameter type can + // be hardcoded to `Geometry&&` again. + using Element = typename GV::template Codim<0>::Entity; + using GeometryRT = decltype(std::declval<Element>().geometry()); + static constexpr bool grtIsReference = std::is_lvalue_reference<GeometryRT>::value; + using GeometryParamType = std::conditional_t<grtIsReference, Geometry, Geometry&&>; public: //! export the type used for global coordinates using GlobalPosition = typename T::GlobalPosition; @@ -73,8 +85,8 @@ public: CCSubControlVolume() = default; - // the contructor in the cc case - CCSubControlVolume(Geometry&& geometry, + // See the explanation above for deriving `GeometryParamType`. + CCSubControlVolume(GeometryParamType geometry, GridIndexType elementIndex) : ParentType() , geometry_(std::move(geometry)) diff --git a/dumux/io/cpgridcreator.hh b/dumux/io/cpgridcreator.hh index 9cec6ee56d6c93dfb054aeb549309f43816a9cc0..8f1e27ef69c011d169a167bb59aa73d7d2c59eca 100644 --- a/dumux/io/cpgridcreator.hh +++ b/dumux/io/cpgridcreator.hh @@ -25,32 +25,23 @@ #define DUMUX_CPGRID_CREATOR_HH #if HAVE_OPM_GRID -#include <dune/grid/CpGrid.hpp> +#include <opm/grid/CpGrid.hpp> #include <opm/parser/eclipse/Parser/Parser.hpp> #include <opm/parser/eclipse/Parser/ParseContext.hpp> #include <opm/parser/eclipse/Deck/Deck.hpp> -#include <dumux/common/properties.hh> #include <dumux/common/parameters.hh> namespace Dumux { -namespace Properties -{ -NEW_PROP_TAG(Scalar); -NEW_PROP_TAG(Grid); -} - /*! * \ingroup InputOutput * \brief A grid creator that reads Petrel files and generates a CpGrid. */ -template <class TypeTag> class CpGridCreator { - using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar); - using Grid = typename GET_PROP_TYPE(TypeTag, Grid); + using Grid = Dune::CpGrid; using GridPointer = std::shared_ptr<Grid>; using Deck = Opm::Deck; @@ -62,8 +53,9 @@ public: { auto fileName = getParam<std::string>("Grid.File"); - deck() = Opm::Parser().parseFile(fileName); - Opm::EclipseGrid ecl_grid(deck()); + static auto deckLocal = Opm::Parser().parseFile(fileName); + Opm::EclipseGrid ecl_grid(deckLocal); + deck() = deckLocal; gridPtr() = std::make_shared<Grid>(*(new Grid())); gridPtr()->processEclipseFormat(ecl_grid, false, false); @@ -102,7 +94,8 @@ public: */ static void loadBalance() { - gridPtr()->loadBalance(); + if (gridPtr()->comm().size() > 1) + gridPtr()->loadBalance(); } }; } diff --git a/patches/README b/patches/README index 9ea8dad5fef0f209535078c991d16661cb2fb7ef..6518ad12e211946b6f3582d3c383b80c4a299baa 100644 --- a/patches/README +++ b/patches/README @@ -5,14 +5,14 @@ patch -p1 <../dumux/patches/istl-2.4.1.patch - If opm-grid has to be used for, e.g., employing the CpGridCreator, - and Opm is compiled from source, it might be necessary to patch opm-parser and opm-common: - patch -p1 <../dumux/patches/opm-parser-2017.10.patch - patch -p1 <../dumux/patches/opm-common-2017.10.patch - In addition, it might be necessary to set manually some CMake variables in the - CMAKE_FLAGS section of the .opts-file: - -DOPM_COMMON_ROOT=/path/to/opm-common/ \ - -Dopm-grid_PREFIX=/path/to/opm-grid \ - -Dopm-common_PREFIX=/path/to/opm-common \ - -Dopm-parser_PREFIX=/path/to/opm-parser \ - -DHAVE_OPM_GRID=1 \ - Currently, Dumux is supposed to be compatible with the Opm 2017.04 and 2017.10 releases. + and Opm is compiled from source, it is necessary to patch opm-common and + opm-grid, namely, to execute in the corresponding folders: + patch -p1 <../dumux/patches/opm-common-2018.04.patch + patch -p1 <../dumux/patches/opm-grid-2018.04.patch + Currently, Dumux is supposed to be compatible with the Opm 2018.04 release. + In addition, it might be necessary to set manually some + CMake variables in the CMAKE_FLAGS section of the .opts-file: + -DOPM_COMMON_ROOT=/path/to/opm-common \ + -Decl_DIR=/path/to/libecl/build \ + -DUSE_MPI=ON \ + diff --git a/patches/opm-common-2017.04.patch b/patches/opm-common-2017.04.patch deleted file mode 100644 index 8d684466e11a802265ff56203165b46747afc9b3..0000000000000000000000000000000000000000 --- a/patches/opm-common-2017.04.patch +++ /dev/null @@ -1,15 +0,0 @@ -diff --git a/cmake/Modules/Findopm-parser.cmake b/cmake/Modules/Findopm-parser.cmake -index a0b0b08..eb142dd 100644 ---- a/cmake/Modules/Findopm-parser.cmake -+++ b/cmake/Modules/Findopm-parser.cmake -@@ -38,8 +38,8 @@ string(REGEX REPLACE "${PROJECT_SOURCE_DIR}/?(.*)" "\\1" BUILD_DIR_SUFFIX "${PR - # or in relative directories to this one - if (OPM_PARSER_ROOT) - set (_no_default_path "NO_DEFAULT_PATH") -- set (_opm_parser_source "") -- set (_opm_parser_build "") -+ set (_opm_parser_source "${OPM_PARSER_ROOT}") -+ set (_opm_parser_build "${OPM_PARSER_ROOT}/${BUILD_DIR_SUFFIX}") - else () - set (_no_default_path "") - set (_opm_parser_source diff --git a/patches/opm-common-2017.10.patch b/patches/opm-common-2017.10.patch deleted file mode 100644 index 0fc80ce2fb1630cf0b915e29159260a06e3b895b..0000000000000000000000000000000000000000 --- a/patches/opm-common-2017.10.patch +++ /dev/null @@ -1,15 +0,0 @@ -diff --git a/cmake/Modules/Findopm-parser.cmake b/cmake/Modules/Findopm-parser.cmake -index 42b5ae9..2977ee9 100644 ---- a/cmake/Modules/Findopm-parser.cmake -+++ b/cmake/Modules/Findopm-parser.cmake -@@ -31,8 +31,8 @@ string(REGEX REPLACE "${PROJECT_SOURCE_DIR}/?(.*)" "\\1" BUILD_DIR_SUFFIX "${PR - # or in relative directories to this one - if (OPM_PARSER_ROOT) - set (_no_default_path "NO_DEFAULT_PATH") -- set (_opm_parser_source "") -- set (_opm_parser_build "") -+ set (_opm_parser_source "${OPM_PARSER_ROOT}") -+ set (_opm_parser_build "${OPM_PARSER_ROOT}/${BUILD_DIR_SUFFIX}") - else () - set (_no_default_path "") - set (_opm_parser_source diff --git a/patches/opm-common-2018.04.patch b/patches/opm-common-2018.04.patch new file mode 100644 index 0000000000000000000000000000000000000000..f20218a28dd19cf850a784855a9098d43c210b72 --- /dev/null +++ b/patches/opm-common-2018.04.patch @@ -0,0 +1,32 @@ +diff --git a/cmake/Modules/OpmProject.cmake b/cmake/Modules/OpmProject.cmake +index 51d572df..68feda08 100644 +--- a/cmake/Modules/OpmProject.cmake ++++ b/cmake/Modules/OpmProject.cmake +@@ -76,6 +76,7 @@ function (opm_cmake_config name) + set (template_dir "${OPM_MACROS_ROOT}/cmake/Templates") + + # write configuration file to locate library ++ set(DUNE_PREFIX ${PROJECT_SOURCE_DIR}) + set(OPM_PROJECT_EXTRA_CODE ${OPM_PROJECT_EXTRA_CODE_INTREE}) + set(PREREQ_LOCATION "${PROJECT_SOURCE_DIR}") + configure_cmake_file (${name} "config" "") +@@ -123,6 +124,7 @@ function (opm_cmake_config name) + # of the build directory (using the same input template) + set(OPM_PROJECT_EXTRA_CODE ${OPM_PROJECT_EXTRA_CODE_INSTALLED}) + set(PREREQ_LOCATION "${CMAKE_INSTALL_PREFIX}/share/opm/cmake/Modules") ++ set(DUNE_PREFIX ${CMAKE_INSTALL_PREFIX}) + configure_cmake_file (${name} "install" "") + configure_vars ( + FILE CMAKE "${PROJECT_BINARY_DIR}/${${name}_NAME}-install.cmake" +diff --git a/cmake/Templates/opm-project-config.cmake.in b/cmake/Templates/opm-project-config.cmake.in +index 421708a8..a783043d 100644 +--- a/cmake/Templates/opm-project-config.cmake.in ++++ b/cmake/Templates/opm-project-config.cmake.in +@@ -24,6 +24,7 @@ if(NOT @opm-project_NAME@_FOUND) + set(CMAKE_MODULE_PATH "${CMAKE_MODULE_PATH}" @PREREQ_LOCATION@) + include(@opm-project_NAME@-prereqs) + # propagate these properties from one build system to the other ++ set (@opm-project_NAME@_PREFIX "@DUNE_PREFIX@") + set (@opm-project_NAME@_VERSION "@opm-project_VERSION@") + set (@opm-project_NAME@_DEFINITIONS "@opm-project_DEFINITIONS@") + set (@opm-project_NAME@_INCLUDE_DIRS "@opm-project_INCLUDE_DIRS@") diff --git a/patches/opm-core-2017.04.patch b/patches/opm-core-2017.04.patch deleted file mode 100644 index 103d09da17880cf4ad873c41794d5f4112aa152d..0000000000000000000000000000000000000000 --- a/patches/opm-core-2017.04.patch +++ /dev/null @@ -1,49 +0,0 @@ -diff --git a/CMakeLists_files.cmake b/CMakeLists_files.cmake -index fbbdaf05..34cdc425 100644 ---- a/CMakeLists_files.cmake -+++ b/CMakeLists_files.cmake -@@ -36,7 +36,7 @@ list (APPEND MAIN_SOURCE_FILES - opm/core/flowdiagnostics/TofReorder.cpp - opm/core/linalg/LinearSolverFactory.cpp - opm/core/linalg/LinearSolverInterface.cpp -- opm/core/linalg/LinearSolverIstl.cpp -+# opm/core/linalg/LinearSolverIstl.cpp - opm/core/linalg/LinearSolverPetsc.cpp - opm/core/linalg/LinearSolverUmfpack.cpp - opm/core/linalg/call_umfpack.c -@@ -133,12 +133,12 @@ list (APPEND TEST_SOURCE_FILES - tests/test_nonuniformtablelinear.cpp - tests/test_parallelistlinformation.cpp - tests/test_sparsevector.cpp -- tests/test_velocityinterpolation.cpp -+ tests/test_velocityinterpolation.cpp - tests/test_uniformtablelinear.cpp - tests/test_wells.cpp - tests/test_wachspresscoord.cpp -- tests/test_linearsolver.cpp -- tests/test_parallel_linearsolver.cpp -+# tests/test_linearsolver.cpp -+# tests/test_parallel_linearsolver.cpp - tests/test_param.cpp - tests/test_satfunc.cpp - tests/test_shadow.cpp -@@ -153,7 +153,7 @@ list (APPEND TEST_SOURCE_FILES - tests/test_anisotropiceikonal.cpp - tests/test_stoppedwells.cpp - tests/test_relpermdiagnostics.cpp -- tests/test_norne_pvt.cpp -+ tests/test_norne_pvt.cpp - ) - - # originally generated with the command: -@@ -194,8 +194,8 @@ list (APPEND TEST_DATA_FILES - list (APPEND EXAMPLE_SOURCE_FILES - examples/compute_eikonal_from_files.cpp - examples/compute_initial_state.cpp -- examples/compute_tof.cpp -- examples/compute_tof_from_files.cpp -+# examples/compute_tof.cpp -+# examples/compute_tof_from_files.cpp - examples/diagnose_relperm.cpp - tutorials/tutorial1.cpp - tutorials/tutorial2.cpp diff --git a/patches/opm-grid-2017.04.patch b/patches/opm-grid-2017.04.patch deleted file mode 100644 index f0a83a2bec834fe9c0f4278994dbb94bbd1ab54b..0000000000000000000000000000000000000000 --- a/patches/opm-grid-2017.04.patch +++ /dev/null @@ -1,132 +0,0 @@ -diff --git a/dune/grid/cpgrid/Entity.hpp b/dune/grid/cpgrid/Entity.hpp -index eadcb33..0efcedd 100644 ---- a/dune/grid/cpgrid/Entity.hpp -+++ b/dune/grid/cpgrid/Entity.hpp -@@ -148,7 +148,7 @@ namespace Dune - } - - /// Returns the geometry of the entity (does not depend on its orientation). -- const Geometry& geometry() const; -+ Geometry geometry() const; - - /// We do not support refinement, so level() is always 0. - int level() const -@@ -182,7 +182,7 @@ namespace Dune - /// The count of subentities of codimension cc - unsigned int subEntities ( const unsigned int cc ) const - { -- static_assert(codim == 0, ""); -+// static_assert(codim == 0, ""); - if (cc == 0) { - return 1; - } else if (cc == 3) { -@@ -412,7 +412,7 @@ namespace Dune { - namespace cpgrid { - - template <int codim> --const typename Entity<codim>::Geometry& Entity<codim>::geometry() const -+typename Entity<codim>::Geometry Entity<codim>::geometry() const - { - return pgrid_->geomVector<codim>()[*this]; - } -diff --git a/dune/grid/cpgrid/Geometry.hpp b/dune/grid/cpgrid/Geometry.hpp -index 86d82a5..10c7e14 100644 ---- a/dune/grid/cpgrid/Geometry.hpp -+++ b/dune/grid/cpgrid/Geometry.hpp -@@ -41,8 +41,13 @@ - - #include <dune/common/version.hh> - #include <dune/geometry/referenceelements.hh> -+ -+#if DUNE_VERSION_NEWER(DUNE_COMMON, 2, 5 ) -+#include <dune/geometry/type.hh> -+#else - #include <dune/geometry/genericgeometry/geometrytraits.hh> - #include <dune/geometry/genericgeometry/matrixhelper.hh> -+#endif - - #include <opm/common/utility/platform_dependent/reenable_warnings.h> - -@@ -99,6 +104,12 @@ namespace Dune - /// Type of the inverse of the transposed Jacobian matrix - typedef FieldMatrix< ctype, coorddimension, mydimension > JacobianInverseTransposed; - -+#if DUNE_VERSION_NEWER(DUNE_GRID,2,5) -+ typedef Dune::Impl::FieldMatrixHelper< double > MatrixHelperType; -+#else -+ typedef Dune::GenericGeometry::MatrixHelper< Dune::GenericGeometry::DuneCoordTraits<double> > MatrixHelperType; -+#endif -+ - /// @brief Construct from centroid, volume (1- and 0-moments) and - /// corners. - /// @param pos the centroid of the entity -@@ -189,12 +200,11 @@ namespace Dune - LocalCoordinate x = refElement.position(0,0); - LocalCoordinate dx; - do { -- using namespace GenericGeometry; - // DF^n dx^n = F^n, x^{n+1} -= dx^n - JacobianTransposed JT = jacobianTransposed(x); - GlobalCoordinate z = global(x); - z -= y; -- MatrixHelper<DuneCoordTraits<double> >::template xTRightInvA<3, 3>(JT, z, dx ); -+ MatrixHelperType::template xTRightInvA<3, 3>(JT, z, dx ); - x -= dx; - } while (dx.two_norm2() > epsilon*epsilon); - return x; -@@ -206,9 +216,8 @@ namespace Dune - /// and {u_j} are the reference coordinates. - double integrationElement(const LocalCoordinate& local_coord) const - { -- FieldMatrix<ctype, coorddimension, mydimension> Jt = jacobianTransposed(local_coord); -- using namespace GenericGeometry; -- return MatrixHelper<DuneCoordTraits<double> >::template sqrtDetAAT<3, 3>(Jt); -+ JacobianTransposed Jt = jacobianTransposed(local_coord); -+ return MatrixHelperType::template sqrtDetAAT<3, 3>(Jt); - } - - /// Using the cube type for all entities now (cells and vertices), -diff --git a/dune/grid/cpgrid/Intersection.hpp b/dune/grid/cpgrid/Intersection.hpp -index f7152ef..cf01473 100644 ---- a/dune/grid/cpgrid/Intersection.hpp -+++ b/dune/grid/cpgrid/Intersection.hpp -@@ -188,7 +188,7 @@ namespace Dune - /// @brief - /// @todo Doc me! - /// @return -- const Geometry& geometry() const -+ Geometry geometry() const - { - return global_geom_; - } -diff --git a/dune/grid/polyhedralgrid/geometry.hh b/dune/grid/polyhedralgrid/geometry.hh -index 64a15fc..9c3c0b4 100644 ---- a/dune/grid/polyhedralgrid/geometry.hh -+++ b/dune/grid/polyhedralgrid/geometry.hh -@@ -42,6 +42,12 @@ namespace Dune - typedef typename Grid::Traits::ExtraData ExtraData; - typedef typename Grid::Traits::template Codim<codimension>::EntitySeed EntitySeed; - -+#if DUNE_VERSION_NEWER(DUNE_GRID,2,5) -+ typedef Dune::Impl::FieldMatrixHelper< double > MatrixHelperType; -+#else -+ typedef Dune::GenericGeometry::MatrixHelper< Dune::GenericGeometry::DuneCoordTraits<double> > MatrixHelperType; -+#endif -+ - explicit PolyhedralGridBasicGeometry ( ExtraData data ) - : data_( data ), - seed_( ) -@@ -136,12 +142,11 @@ namespace Dune - LocalCoordinate x = refElement.position(0,0); - LocalCoordinate dx; - do { -- using namespace GenericGeometry; - // DF^n dx^n = F^n, x^{n+1} -= dx^n - JacobianTransposed JT = jacobianTransposed(x); - GlobalCoordinate z = global(x); - z -= y; -- MatrixHelper<DuneCoordTraits<double> >::template xTRightInvA<mydimension,coorddimension>(JT, z, dx ); -+ MatrixHelperType::template xTRightInvA<mydimension,coorddimension>(JT, z, dx ); - x -= dx; - } while (dx.two_norm2() > epsilon*epsilon); - return x; diff --git a/patches/opm-grid-2018.04.patch b/patches/opm-grid-2018.04.patch new file mode 100644 index 0000000000000000000000000000000000000000..8b6c34f4b51622a3541a9b4f4871fddfa5cab36e --- /dev/null +++ b/patches/opm-grid-2018.04.patch @@ -0,0 +1,12 @@ +diff --git a/opm/grid/common/p2pcommunicator.hh b/opm/grid/common/p2pcommunicator.hh +index 1b4b327..55a2b2f 100644 +--- a/opm/grid/common/p2pcommunicator.hh ++++ b/opm/grid/common/p2pcommunicator.hh +@@ -23,6 +23,7 @@ + #include <vector> + #include <set> + #include <map> ++#include <cassert> + + #include <dune/common/version.hh> + diff --git a/patches/opm-parser-2017.04.patch b/patches/opm-parser-2017.04.patch deleted file mode 100644 index a9fbd643fe875e9c131348ed8cd926203b21ebbf..0000000000000000000000000000000000000000 --- a/patches/opm-parser-2017.04.patch +++ /dev/null @@ -1,45 +0,0 @@ -diff --git a/CMakeLists.txt b/CMakeLists.txt -index a53c6634..2a6b5147 100644 ---- a/CMakeLists.txt -+++ b/CMakeLists.txt -@@ -99,24 +99,8 @@ SET( CMAKE_C_FLAGS_DEBUG "-O0 -DDEBUG -ggdb3 ${CMAKE_C_FLAGS} ${CMAKE_ - SET( CMAKE_C_FLAGS_RELEASE "-O3 -DNDEBUG -mtune=native ${CMAKE_C_FLAGS} ${CMAKE_C_FLAGS_RELEASE}") - - --if (BOOST_ROOT AND NOT DEFINED Boost_NO_SYSTEM_PATHS) -- set (Boost_NO_SYSTEM_PATHS TRUE) --endif (BOOST_ROOT AND NOT DEFINED Boost_NO_SYSTEM_PATHS) -- --# if we are building shared libraries ourselves, then don't include Boost in them --if (BUILD_SHARED_LIBS) -- set(Boost_USE_STATIC_LIBS OFF) --elseif (DEFINED BUILD_SHARED_LIBS) -- set(Boost_USE_STATIC_LIBS ON) --endif () --set(Boost_USE_MULTITHREADED ON) --set(Boost_USE_STATIC_RUNTIME OFF) -- -- --add_definitions(-DOPM_PARSER_DECK_API=1) --# Requires BOOST filesystem version 3, thus 1.44 is necessary. --add_definitions(-DBOOST_FILESYSTEM_VERSION=3) --find_package(Boost 1.44.0 COMPONENTS filesystem date_time system unit_test_framework regex REQUIRED) -+# get the prerequisite Boost libraries -+find_package(Boost 1.44.0 COMPONENTS filesystem date_time system unit_test_framework regex ${OPM_PARSER_QUIET}) - include_directories(SYSTEM ${Boost_INCLUDE_DIRS}) - - include_directories(BEFORE ${PROJECT_SOURCE_DIR}) -@@ -177,9 +161,9 @@ include(OpmProject) - include(ConfigVars) - set(CMAKE_LIBRARY_OUTPUT_DIRECTORY ${CMAKE_BINARY_DIR}/lib) - set(opm-parser_NAME opm-parser) --set(opm-parser_LIBRARIES opmparser opmjson) --if(NOT BUILD_SHARED_LIBS) -- list(APPEND opm-parser_LIBRARIES cjson) --endif() -+#set(opm-parser_LIBRARIES opmparser opmjson) -+#if(NOT BUILD_SHARED_LIBS) -+# list(APPEND opm-parser_LIBRARIES cjson) -+#endif() - set(opm-parser_TARGET opmparser) - opm_cmake_config(opm-parser) diff --git a/patches/opm-parser-2017.10.patch b/patches/opm-parser-2017.10.patch deleted file mode 100644 index 77f691d9b2ea9bb01562de931c3042e9f7f09f71..0000000000000000000000000000000000000000 --- a/patches/opm-parser-2017.10.patch +++ /dev/null @@ -1,20 +0,0 @@ -diff --git a/CMakeLists.txt b/CMakeLists.txt -index d81c5dc2..964c19de 100644 ---- a/CMakeLists.txt -+++ b/CMakeLists.txt -@@ -91,14 +91,7 @@ else () - add_subdirectory(external/cjson) - endif() - --# if building shared libraries, then don't include Boost in them --if (BUILD_SHARED_LIBS) -- set(Boost_USE_STATIC_LIBS OFF) -- set(Boost_USE_STATIC_RUNTIME OFF) --else () --# if using dynamic boost, the header file must generate a main() function -- set(Boost_USE_STATIC_LIBS ON) --endif () -+add_definitions(-DBOOST_TEST_DYN_LINK) - - find_package(Boost 1.44.0 - COMPONENTS filesystem diff --git a/test/porousmediumflow/2p/implicit/CMakeLists.txt b/test/porousmediumflow/2p/implicit/CMakeLists.txt index b18c82cdbfb9e845d52c5cee9ba8922cf4253df7..069135d12e5e6bf1150eda5ab6d8bda8df308c8d 100644 --- a/test/porousmediumflow/2p/implicit/CMakeLists.txt +++ b/test/porousmediumflow/2p/implicit/CMakeLists.txt @@ -1,5 +1,6 @@ add_subdirectory(adaptive) add_subdirectory(boxdfm) +add_subdirectory(cornerpoint) add_subdirectory(fracture) add_subdirectory(incompressible) add_subdirectory(nonisothermal) diff --git a/test/porousmediumflow/2p/implicit/cornerpoint/CMakeLists.txt b/test/porousmediumflow/2p/implicit/cornerpoint/CMakeLists.txt new file mode 100644 index 0000000000000000000000000000000000000000..6ad7ac5bbccbd2b8a63232956a3b6decd624dadf --- /dev/null +++ b/test/porousmediumflow/2p/implicit/cornerpoint/CMakeLists.txt @@ -0,0 +1,21 @@ +dune_symlink_to_source_files(FILES "test_2p_cornerpoint.input") +dune_symlink_to_source_files(FILES grids) + +dune_add_test(NAME test_2p_cornerpoint + SOURCES test_2p_cornerpoint.cc + CMAKE_GUARD HAVE_OPM_GRID + COMPILE_DEFINITIONS HAVE_ECL_INPUT=1 + COMMAND ${CMAKE_SOURCE_DIR}/bin/testing/runtest.py + CMD_ARGS --script fuzzy + --files ${CMAKE_SOURCE_DIR}/test/references/cc2pcornerpoint-reference.vtu + ${CMAKE_CURRENT_BINARY_DIR}/cc2pcornerpoint-00005.vtu + --command "${CMAKE_CURRENT_BINARY_DIR}/test_2p_cornerpoint -Problem.Name cc2pcornerpoint") + +set(CMAKE_BUILD_TYPE Release) + +#install sources +install(FILES +problem.hh +spatialparams.hh +test_2p_cornerpoint.cc +DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/dumux/test/implicit/2p/cornerpoint) diff --git a/test/porousmediumflow/2p/implicit/cornerpoint/grids/hardsmall.grdecl b/test/porousmediumflow/2p/implicit/cornerpoint/grids/hardsmall.grdecl new file mode 100644 index 0000000000000000000000000000000000000000..5ac8ca6029c0449f3b5330dea05feba185ba0f28 --- /dev/null +++ b/test/porousmediumflow/2p/implicit/cornerpoint/grids/hardsmall.grdecl @@ -0,0 +1,38 @@ +MAPUNITS +'METRES' +/ + +SPECGRID +2 1 2 1 F +/ + +COORD +0 0 0 0 0 2 +1 0 0 1 0 2 +2 0 0 2 0 2 +0 1 0 0 1 2 +1 1 0 1 1 2 +2 1 0 2 1 2 +/ + +ZCORN +0 0 0 0 0 0 0 0 +1 1 0.8 1.2 1 1 1.2 1.2 +1 1 0.8 1.2 1 1 1.2 1.2 +2 2 2.2 2 2 2 1.9 2 +/ + +ACTNUM +1 1 +1 1 +/ + +PORO +0.1 0.1 +0.3 0.3 +/ + +PERMX +100 100 +500 500 +/ diff --git a/test/porousmediumflow/2p/implicit/cornerpoint/problem.hh b/test/porousmediumflow/2p/implicit/cornerpoint/problem.hh new file mode 100644 index 0000000000000000000000000000000000000000..dd80e513066eeed9819cf73b4be454581355eab1 --- /dev/null +++ b/test/porousmediumflow/2p/implicit/cornerpoint/problem.hh @@ -0,0 +1,260 @@ +// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- +// vi: set et ts=4 sw=4 sts=4: +/***************************************************************************** + * See the file COPYING for full copying permissions. * + * * + * This program 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 2 of the License, or * + * (at your option) any later version. * + * * + * This program 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 this program. If not, see <http://www.gnu.org/licenses/>. * + *****************************************************************************/ +/*! + * \ingroup TwoPTests + * \brief The properties for the 2p cornerpoint test + */ +#ifndef DUMUX_TWOP_CORNERPOINT_TEST_PROBLEM_HH +#define DUMUX_TWOP_CORNERPOINT_TEST_PROBLEM_HH + +#include <opm/grid/CpGrid.hpp> + +#include <dumux/io/cpgridcreator.hh> +#include <dumux/discretization/cellcentered/tpfa/properties.hh> + +#include <dumux/material/components/trichloroethene.hh> +#include <dumux/material/components/simpleh2o.hh> +#include <dumux/material/fluidsystems/1pliquid.hh> +#include <dumux/material/fluidsystems/2pimmiscible.hh> + +#include <dumux/porousmediumflow/problem.hh> +#include <dumux/porousmediumflow/2p/model.hh> +#include <dumux/porousmediumflow/2p/incompressiblelocalresidual.hh> + +#include "spatialparams.hh" + +namespace Dumux { +// forward declarations +template<class TypeTag> class TwoPCornerPointTestProblem; + +namespace Properties { +NEW_TYPE_TAG(TwoPCornerPoint, INHERITS_FROM(TwoP, CCTpfaModel)); + +// Set the grid type +SET_TYPE_PROP(TwoPCornerPoint, Grid, Dune::CpGrid); + +// Set the problem type +SET_TYPE_PROP(TwoPCornerPoint, Problem, TwoPCornerPointTestProblem<TypeTag>); + +// the local residual containing the analytic derivative methods +SET_TYPE_PROP(TwoPCornerPoint, LocalResidual, TwoPIncompressibleLocalResidual<TypeTag>); + +// Set the fluid system +SET_PROP(TwoPCornerPoint, FluidSystem) +{ + using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar); + using WettingPhase = FluidSystems::OnePLiquid<Scalar, Components::SimpleH2O<Scalar> >; + using NonwettingPhase = FluidSystems::OnePLiquid<Scalar, Components::Trichloroethene<Scalar> >; + using type = FluidSystems::TwoPImmiscible<Scalar, WettingPhase, NonwettingPhase>; +}; + +// Set the spatial parameters +SET_PROP(TwoPCornerPoint, SpatialParams) +{ +private: + using FVGridGeometry = typename GET_PROP_TYPE(TypeTag, FVGridGeometry); + using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar); +public: + using type = TwoPCornerPointTestSpatialParams<FVGridGeometry, Scalar>; +}; + +// Set the grid creator +SET_TYPE_PROP(TwoPCornerPoint, GridCreator, CpGridCreator); + +// Enable caching +SET_BOOL_PROP(TwoPCornerPoint, EnableGridVolumeVariablesCache, false); +SET_BOOL_PROP(TwoPCornerPoint, EnableGridFluxVariablesCache, false); +SET_BOOL_PROP(TwoPCornerPoint, EnableFVGridGeometryCache, false); +} // end namespace Properties + +/*! + * \ingroup TwoPTests + * \brief The incompressible 2p cornerpoint test problem. + */ +template<class TypeTag> +class TwoPCornerPointTestProblem : public PorousMediumFlowProblem<TypeTag> +{ + using ParentType = PorousMediumFlowProblem<TypeTag>; + using GridView = typename GET_PROP_TYPE(TypeTag, GridView); + using Element = typename GridView::template Codim<0>::Entity; + using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar); + using FluidSystem = typename GET_PROP_TYPE(TypeTag, FluidSystem); + using PrimaryVariables = typename GET_PROP_TYPE(TypeTag, PrimaryVariables); + using ElementVolumeVariables = typename GET_PROP_TYPE(TypeTag, GridVolumeVariables)::LocalView; + using FVGridGeometry = typename GET_PROP_TYPE(TypeTag, FVGridGeometry); + using FVElementGeometry = typename FVGridGeometry::LocalView; + using SubControlVolumeFace = typename FVElementGeometry::SubControlVolumeFace; + using SubControlVolume = typename FVElementGeometry::SubControlVolume; + using BoundaryTypes = typename GET_PROP_TYPE(TypeTag, BoundaryTypes); + using GlobalPosition = typename Element::Geometry::GlobalCoordinate; + using NumEqVector = typename GET_PROP_TYPE(TypeTag, NumEqVector); + using Indices = typename GET_PROP_TYPE(TypeTag, ModelTraits)::Indices; + using GridCreator = typename GET_PROP_TYPE(TypeTag, GridCreator); + enum { dimWorld = GridView::dimensionworld }; + +public: + TwoPCornerPointTestProblem(std::shared_ptr<const FVGridGeometry> fvGridGeometry) + : ParentType(fvGridGeometry) + { + gravity_ = {0, 0, 9.81}; + injectionElement_ = getParam<int>("Problem.InjectionElement"); + injectionRate_ = getParam<Scalar>("Problem.InjectionRate"); + } + + /*! + * \brief Specifies which kind of boundary condition should be + * used for which equation on a given boundary segment. + * + * \param element The finite element + * \param scvf The sub control volume face + */ + BoundaryTypes boundaryTypes(const Element &element, + const SubControlVolumeFace &scvf) const + { + BoundaryTypes bcTypes; + + // set no-flux on top and bottom, hydrostatic on the rest + // use the scvf normal to decide + const auto& normal = scvf.unitOuterNormal(); + using std::abs; + if (abs(normal[dimWorld-1]) > 0.5) + bcTypes.setAllNeumann(); + else + bcTypes.setAllDirichlet(); + + return bcTypes; + } + + /*! + * \brief Evaluates the boundary conditions for a Dirichlet + * boundary segment + * + * \param values Stores the Dirichlet values for the conservation equations in + * \f$ [ \textnormal{unit of primary variable} ] \f$ + * \param globalPos The global position + */ + PrimaryVariables dirichletAtPos(const GlobalPosition &globalPos) const + { + return initialAtPos(globalPos); + } + + /*! + * \brief Evaluate the boundary conditions for a neumann + * boundary segment. + * + * \param values Stores the Neumann values for the conservation equations in + * \f$ [ \textnormal{unit of conserved quantity} / (m^(dim-1) \cdot s )] \f$ + * \param globalPos The position of the integration point of the boundary segment. + * + * For this method, the \a values parameter stores the mass flux + * in normal direction of each phase. Negative values mean influx. + */ + NumEqVector neumannAtPos(const GlobalPosition &globalPos) const + { + return NumEqVector(0.0); + } + + //! \copydoc Dumux::FVProblem::source() + NumEqVector source(const Element &element, + const FVElementGeometry& fvGeometry, + const ElementVolumeVariables& elemVolVars, + const SubControlVolume &scv) const + { + NumEqVector values(0.0); + + int eIdx = GridCreator::grid().leafGridView().indexSet().index(element); + if (eIdx == injectionElement_) + values[FluidSystem::phase1Idx] = injectionRate_/element.geometry().volume(); + + return values; + } + + const GlobalPosition& gravity() const + { + return gravity_; + } + + /*! + * \brief Evaluates the initial values for a control volume + * + * \param values Stores the initial values for the conservation equations in + * \f$ [ \textnormal{unit of primary variables} ] \f$ + * \param globalPos The global position + */ + PrimaryVariables initialAtPos(const GlobalPosition &globalPos) const + { + PrimaryVariables values; + + // hydrostatic pressure + Scalar densityW = 1000; + values[Indices::pressureIdx] = 1e5 + densityW*(this->gravity()*globalPos); + values[Indices::saturationIdx] = 0.0; + + return values; + } + + /*! + * \brief Returns the temperature \f$\mathrm{[K]}\f$ for an isothermal problem. + * + * This is not specific to the discretization. By default it just + * throws an exception so it must be overloaded by the problem if + * no energy equation is used. + */ + Scalar temperature() const + { + return 293.15; + } + + /*! + * \brief Append all quantities of interest which can be derived + * from the solution of the current time step to the VTK + * writer. + */ + template<class VTKWriter> + void addFieldsToWriter(VTKWriter& vtk) + { + const auto numElements = this->fvGridGeometry().gridView().size(0); + + permX_.resize(numElements); + permZ_.resize(numElements); + + vtk.addField(permX_, "PERMX [mD]"); + vtk.addField(permZ_, "PERMZ [mD]"); + + const auto& gridView = this->fvGridGeometry().gridView(); + for (const auto& element : elements(gridView)) + { + const auto eIdx = this->fvGridGeometry().elementMapper().index(element); + + // transfer output to mD = 9.86923e-16 m^2 + permX_[eIdx] = this->spatialParams().permeabilityX(eIdx)/9.86923e-16; + permZ_[eIdx] = this->spatialParams().permeabilityZ(eIdx)/9.86923e-16; + } + } + +private: + GlobalPosition gravity_; + int injectionElement_; + Scalar injectionRate_; + std::vector<Scalar> permX_, permZ_; +}; + +} // end namespace Dumux + +#endif diff --git a/test/porousmediumflow/2p/implicit/cornerpoint/spatialparams.hh b/test/porousmediumflow/2p/implicit/cornerpoint/spatialparams.hh new file mode 100644 index 0000000000000000000000000000000000000000..66daf1148bbe29aa4f6113a3699ab4696221fd6a --- /dev/null +++ b/test/porousmediumflow/2p/implicit/cornerpoint/spatialparams.hh @@ -0,0 +1,203 @@ +// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- +// vi: set et ts=4 sw=4 sts=4: +/***************************************************************************** + * See the file COPYING for full copying permissions. * + * * + * This program 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 2 of the License, or * + * (at your option) any later version. * + * * + * This program 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 this program. If not, see <http://www.gnu.org/licenses/>. * + *****************************************************************************/ +/*! + * \ingroup TwoPTests + * \brief The spatial params for the 2p cornerpoint test + */ +#ifndef DUMUX_TWOP_CORNERPOINT_TEST_SPATIAL_PARAMS_HH +#define DUMUX_TWOP_CORNERPOINT_TEST_SPATIAL_PARAMS_HH + +#include <dumux/io/cpgridcreator.hh> +#include <dumux/material/spatialparams/fv.hh> +#include <dumux/material/fluidmatrixinteractions/2p/regularizedvangenuchten.hh> +#include <dumux/material/fluidmatrixinteractions/2p/efftoabslaw.hh> + +namespace Dumux { + +/*! + * \ingroup TwoPTests + * \brief The spatial params for the incompressible 2p test + */ +template<class FVGridGeometry, class Scalar> +class TwoPCornerPointTestSpatialParams +: public FVSpatialParams<FVGridGeometry, Scalar, TwoPCornerPointTestSpatialParams<FVGridGeometry, Scalar>> +{ + using GridView = typename FVGridGeometry::GridView; + using Element = typename GridView::template Codim<0>::Entity; + using FVElementGeometry = typename FVGridGeometry::LocalView; + using SubControlVolume = typename FVElementGeometry::SubControlVolume; + using ThisType = TwoPCornerPointTestSpatialParams<FVGridGeometry, Scalar>; + using ParentType = FVSpatialParams<FVGridGeometry, Scalar, ThisType>; + using GridCreator = CpGridCreator; + + static constexpr int dimWorld = GridView::dimensionworld; + using GlobalPosition = typename Element::Geometry::GlobalCoordinate; + + using EffectiveLaw = RegularizedVanGenuchten<Scalar>; + + using DimWorldMatrix = Dune::FieldMatrix<Scalar, dimWorld, dimWorld>; + +public: + using MaterialLaw = EffToAbsLaw<EffectiveLaw>; + using MaterialLawParams = typename MaterialLaw::Params; + using PermeabilityType = DimWorldMatrix; + + TwoPCornerPointTestSpatialParams(std::shared_ptr<const FVGridGeometry> fvGridGeometry) + : ParentType(fvGridGeometry) + { + homogeneous_ = getParam<bool>("Problem.Homogeneous"); + + const std::vector<int>& globalCell = GridCreator::grid().globalCell(); + + if (GridCreator::deck().hasKeyword("PORO")) { + std::cout << "Found PORO..." << std::endl; + std::vector<double> eclVector = GridCreator::deck().getKeyword("PORO").getRawDoubleData(); + porosity_.resize(globalCell.size()); + + for (size_t i = 0; i < globalCell.size(); ++i) { + if (homogeneous_) + porosity_[i] = 0.2; + else + porosity_[i] = eclVector[globalCell[i]]; + } + } + + if (GridCreator::deck().hasKeyword("PERMX")) { + std::cout << "Found PERMX..." << std::endl; + std::vector<double> eclVector = GridCreator::deck().getKeyword("PERMX").getRawDoubleData(); + permX_.resize(globalCell.size()); + + for (size_t i = 0; i < globalCell.size(); ++i) { + // assume that values are given in mD = 9.86923e-16 m^2 + if (homogeneous_) + permX_[i] = 9.86923e-16; + else + permX_[i] = 9.86923e-16*eclVector[globalCell[i]]; + } + } + + if (GridCreator::deck().hasKeyword("PERMZ")) { + std::cout << "Found PERMZ..." << std::endl; + std::vector<double> eclVector = GridCreator::deck().getKeyword("PERMZ").getRawDoubleData(); + permZ_.resize(globalCell.size()); + + for (size_t i = 0; i < globalCell.size(); ++i) { + // assume that values are given in mD = 9.86923e-16 m^2 + if (homogeneous_) + permZ_[i] = 9.86923e-16; + else + permZ_[i] = 9.86923e-16*eclVector[globalCell[i]]; + } + } + else { + std::cout << "PERMZ not found, set equal to PERMX..." << std::endl; + permZ_ = permX_; + } + + // parameters for the Van Genuchten law + materialParams_.setVgAlpha(0.00045); + materialParams_.setVgn(7.3); + } + + /*! + * \brief Function for defining the (intrinsic) permeability \f$[m^2]\f$. + * In this test, we use element-wise distributed permeabilities. + * + * \param element The current element + * \param scv The sub-control volume inside the element. + * \param elemSol The solution at the dofs connected to the element. + * \return permeability + */ + template<class ElementSolution> + PermeabilityType permeability(const Element& element, + const SubControlVolume& scv, + const ElementSolution& elemSol) const + { + int eIdx = GridCreator::grid().leafGridView().indexSet().index(element); + + PermeabilityType K(0); + K[0][0] = K[1][1] = permX_[eIdx]; + K[2][2] = permZ_[eIdx]; + + return K; + } + + /*! + * \brief Returns the porosity \f$[-]\f$ + * + * \param element The current element + * \param scv The sub-control volume inside the element. + * \param elemSol The solution at the dofs connected to the element. + * \return porosity + */ + template<class ElementSolution> + Scalar porosity(const Element& element, + const SubControlVolume& scv, + const ElementSolution& elemSol) const + { + int eIdx = GridCreator::grid().leafGridView().indexSet().index(element); + + return porosity_[eIdx]; + } + + /*! + * \brief Returns the parameter object for the material law. + * + * \param element The current element + * \param scv The sub-control volume inside the element. + * \param elemSol The solution at the dofs connected to the element. + * \return the material parameters object + */ + template<class ElementSolution> + const MaterialLawParams& materialLawParams(const Element& element, + const SubControlVolume& scv, + const ElementSolution& elemSol) const + { + return materialParams_; + } + + /*! + * \brief Function for defining which phase is to be considered as the wetting phase. + * + * \return the wetting phase index + * \param globalPos The global position + */ + template<class FluidSystem> + int wettingPhaseAtPos(const GlobalPosition& globalPos) const + { + return FluidSystem::phase0Idx; + } + + Scalar permeabilityX(int eIdx) + { return permX_[eIdx]; } + + Scalar permeabilityZ(int eIdx) + { return permZ_[eIdx]; } + +private: + MaterialLawParams materialParams_; + std::vector<Scalar> porosity_; + std::vector<Scalar> permX_; + std::vector<Scalar> permZ_; + bool homogeneous_; +}; + +} // end namespace Dumux + +#endif diff --git a/test/porousmediumflow/2p/implicit/cornerpoint/test_2p_cornerpoint.cc b/test/porousmediumflow/2p/implicit/cornerpoint/test_2p_cornerpoint.cc new file mode 100644 index 0000000000000000000000000000000000000000..b291a8eef661ea81e9c9da28389ff61cd4556b5f --- /dev/null +++ b/test/porousmediumflow/2p/implicit/cornerpoint/test_2p_cornerpoint.cc @@ -0,0 +1,246 @@ +// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- +// vi: set et ts=4 sw=4 sts=4: +/***************************************************************************** + * See the file COPYING for full copying permissions. * + * * + * This program 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 2 of the License, or * + * (at your option) any later version. * + * * + * This program 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 this program. If not, see <http://www.gnu.org/licenses/>. * + *****************************************************************************/ +/*! + * \file + * + * \brief test for the two-phase porousmedium flow model on a cornerpoint grid + */ +#include <config.h> +// As of versions Opm 2018.04 and Dune 2.6, using opm-grid and dune-uggrid +// together results in a compiler error. The following lines enforce that +// dune-uggrid isn't considered. +#if HAVE_UG +#undef HAVE_UG +#define HAVE_UG 0 +#endif + +#include <ctime> +#include <iostream> + +#include <dune/common/parallel/mpihelper.hh> +#include <dune/common/timer.hh> +#include <dune/grid/io/file/dgfparser/dgfexception.hh> +#include <dune/grid/io/file/vtk.hh> +#include <dune/istl/io.hh> + +#include "problem.hh" + +#include <dumux/common/properties.hh> +#include <dumux/common/parameters.hh> +#include <dumux/common/valgrind.hh> +#include <dumux/common/dumuxmessage.hh> +#include <dumux/common/defaultusagemessage.hh> + +#include <dumux/linear/amgbackend.hh> +#include <dumux/nonlinear/newtonsolver.hh> + +#include <dumux/assembly/fvassembler.hh> +#include <dumux/assembly/diffmethod.hh> + +#include <dumux/discretization/methods.hh> + +#include <dumux/io/vtkoutputmodule.hh> + +/*! + * \brief Provides an interface for customizing error messages associated with + * reading in parameters. + * + * \param progName The name of the program, that was tried to be started. + * \param errorMsg The error message that was issued by the start function. + * Comprises the thing that went wrong and a general help message. + */ +void usage(const char *progName, const std::string &errorMsg) +{ + if (errorMsg.size() > 0) { + std::string errorMessageOut = "\nUsage: "; + errorMessageOut += progName; + errorMessageOut += " [options]\n"; + errorMessageOut += errorMsg; + errorMessageOut += "\n\nThe list of mandatory arguments for this program is:\n" + "\t-TimeManager.TEnd End of the simulation [s] \n" + "\t-TimeManager.DtInitial Initial timestep size [s] \n" + "\t-Grid.LowerLeft Lower left corner coordinates\n" + "\t-Grid.UpperRight Upper right corner coordinates\n" + "\t-Grid.Cells Number of cells in respective coordinate directions\n" + "\t definition in DGF format\n" + "\t-SpatialParams.LensLowerLeft coordinates of the lower left corner of the lens [m] \n" + "\t-SpatialParams.LensUpperRight coordinates of the upper right corner of the lens [m] \n" + "\t-SpatialParams.Permeability Permeability of the domain [m^2] \n" + "\t-SpatialParams.PermeabilityLens Permeability of the lens [m^2] \n"; + + std::cout << errorMessageOut + << "\n"; + } +} + +int main(int argc, char** argv) try +{ + using namespace Dumux; + + // define the type tag for this problem + using TypeTag = TTAG(TwoPCornerPoint); + + // initialize MPI, finalize is done automatically on exit + const auto& mpiHelper = Dune::MPIHelper::instance(argc, argv); + + // print dumux start message + if (mpiHelper.rank() == 0) + DumuxMessage::print(/*firstCall=*/true); + + // parse command line arguments and input file + Parameters::init(argc, argv, usage); + + // try to create a grid (from the given grid file or the input file) + using GridCreator = typename GET_PROP_TYPE(TypeTag, GridCreator); + GridCreator::makeGrid(); + GridCreator::loadBalance(); + + //////////////////////////////////////////////////////////// + // run instationary non-linear problem on this grid + //////////////////////////////////////////////////////////// + + // we compute on the leaf grid view + const auto& leafGridView = GridCreator::grid().leafGridView(); + + // create the finite volume grid geometry + using FVGridGeometry = typename GET_PROP_TYPE(TypeTag, FVGridGeometry); + auto fvGridGeometry = std::make_shared<FVGridGeometry>(leafGridView); + fvGridGeometry->update(); + + // the problem (initial and boundary conditions) + using Problem = typename GET_PROP_TYPE(TypeTag, Problem); + auto problem = std::make_shared<Problem>(fvGridGeometry); + + // the solution vector + using SolutionVector = typename GET_PROP_TYPE(TypeTag, SolutionVector); + SolutionVector x(fvGridGeometry->numDofs()); + problem->applyInitialSolution(x); + auto xOld = x; + + // the grid variables + using GridVariables = typename GET_PROP_TYPE(TypeTag, GridVariables); + auto gridVariables = std::make_shared<GridVariables>(problem, fvGridGeometry); + gridVariables->init(x, xOld); + + // get some time loop parameters + using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar); + const auto tEnd = getParam<Scalar>("TimeLoop.TEnd"); + const auto maxDt = getParam<Scalar>("TimeLoop.MaxTimeStepSize"); + auto dt = getParam<Scalar>("TimeLoop.DtInitial"); + + // check if we are about to restart a previously interrupted simulation + Scalar restartTime = 0; + if (Parameters::getTree().hasKey("Restart") || Parameters::getTree().hasKey("TimeLoop.Restart")) + restartTime = getParam<Scalar>("TimeLoop.Restart"); + + // intialize the vtk output module + using VtkOutputFields = typename GET_PROP_TYPE(TypeTag, VtkOutputFields); + + VtkOutputModule<TypeTag> vtkWriter(*problem, *fvGridGeometry, *gridVariables, + x, problem->name(), "", Dune::VTK::conforming); + + VtkOutputFields::init(vtkWriter); //!< Add model specific output fields + problem->addFieldsToWriter(vtkWriter); //!< Add some more problem dependent fields + vtkWriter.write(0.0); + + // instantiate time loop + auto timeLoop = std::make_shared<TimeLoop<Scalar>>(restartTime, dt, tEnd); + timeLoop->setMaxTimeStepSize(maxDt); + + // the assembler with time loop for instationary problem + using Assembler = FVAssembler<TypeTag, DiffMethod::numeric>; + auto assembler = std::make_shared<Assembler>(problem, fvGridGeometry, gridVariables, timeLoop); + + // the linear solver + using LinearSolver = AMGBackend<TypeTag>; + auto linearSolver = std::make_shared<LinearSolver>(leafGridView, fvGridGeometry->dofMapper()); + + // the non-linear solver + using NewtonSolver = Dumux::NewtonSolver<Assembler, LinearSolver>; + NewtonSolver nonLinearSolver(assembler, linearSolver); + + // time loop + timeLoop->start(); do + { + // set previous solution for storage evaluations + assembler->setPreviousSolution(xOld); + + // solve the non-linear system with time step control + nonLinearSolver.solve(x, *timeLoop); + + // make the new solution the old solution + xOld = x; + gridVariables->advanceTimeStep(); + + // advance to the time loop to the next step + timeLoop->advanceTimeStep(); + + // write vtk output + vtkWriter.write(timeLoop->time()); + + // report statistics of this time step + timeLoop->reportTimeStep(); + + // set new dt as suggested by the Newton solver + timeLoop->setTimeStepSize(nonLinearSolver.suggestTimeStepSize(timeLoop->timeStepSize())); + + } while (!timeLoop->finished()); + + // output some Newton statistics + nonLinearSolver.report(); + + timeLoop->finalize(leafGridView.comm()); + + //////////////////////////////////////////////////////////// + // finalize, print dumux message to say goodbye + //////////////////////////////////////////////////////////// + + // print dumux end message + if (mpiHelper.rank() == 0) + { + Parameters::print(); + DumuxMessage::print(/*firstCall=*/false); + } + + return 0; +} // end main +catch (Dumux::ParameterException &e) +{ + std::cerr << std::endl << e << " ---> Abort!" << std::endl; + return 1; +} +catch (Dune::DGFException & e) +{ + std::cerr << "DGF exception thrown (" << e << + "). Most likely, the DGF file name is wrong " + "or the DGF file is corrupted, " + "e.g. missing hash at end of file or wrong number (dimensions) of entries." + << " ---> Abort!" << std::endl; + return 2; +} +catch (Dune::Exception &e) +{ + std::cerr << "Dune reported error: " << e << " ---> Abort!" << std::endl; + return 3; +} +catch (...) +{ + std::cerr << "Unknown exception thrown! ---> Abort!" << std::endl; + return 4; +} diff --git a/test/porousmediumflow/2p/implicit/cornerpoint/test_2p_cornerpoint.input b/test/porousmediumflow/2p/implicit/cornerpoint/test_2p_cornerpoint.input new file mode 100644 index 0000000000000000000000000000000000000000..b4e0ef3e375034edd8372993ecbe524b4fc39e14 --- /dev/null +++ b/test/porousmediumflow/2p/implicit/cornerpoint/test_2p_cornerpoint.input @@ -0,0 +1,12 @@ +[TimeLoop] +DtInitial = 1e3 # [s] +TEnd = 1e4 # [s] + +[Grid] +File = ./grids/hardsmall.grdecl + +[Problem] +Name = cc2pcornerpoint # name passed to the output routines +Homogeneous = 0 # use a homogeneous permeablity/porosity distribution +InjectionElement = 0 # index of the element where injection should take place +InjectionRate = 0.1 # injection rate in kg/s diff --git a/test/references/cc2pcornerpoint-reference.vtu b/test/references/cc2pcornerpoint-reference.vtu index 4a5c8f7cd814f55abe5811656dd80e227abeb177..fc413fa086e360c9a32975159afeffa16df5f807 100644 --- a/test/references/cc2pcornerpoint-reference.vtu +++ b/test/references/cc2pcornerpoint-reference.vtu @@ -2,39 +2,36 @@ <VTKFile type="UnstructuredGrid" version="0.1" byte_order="LittleEndian"> <UnstructuredGrid> <Piece NumberOfCells="4" NumberOfPoints="22"> - <CellData Scalars="Sn"> - <DataArray type="Float32" Name="Sn" NumberOfComponents="1" format="ascii"> - 0.714401 0.349245 0.297318 0.0336473 + <CellData Scalars="porosity"> + <DataArray type="Float32" Name="porosity" NumberOfComponents="1" format="ascii"> + 0.1 0.1 0.3 0.3 </DataArray> - <DataArray type="Float32" Name="Sw" NumberOfComponents="1" format="ascii"> - 0.285599 0.650755 0.702682 0.966353 + <DataArray type="Float32" Name="pc" NumberOfComponents="1" format="ascii"> + 2614.11 2091.05 2024.13 1414.34 </DataArray> - <DataArray type="Float32" Name="pn" NumberOfComponents="1" format="ascii"> - 186999 121275 123879 117922 + <DataArray type="Float32" Name="Sw" NumberOfComponents="1" format="ascii"> + 0.285556 0.652038 0.702397 0.96918 </DataArray> <DataArray type="Float32" Name="pw" NumberOfComponents="1" format="ascii"> - 184385 119183 121856 116490 - </DataArray> - <DataArray type="Float32" Name="pc" NumberOfComponents="1" format="ascii"> - 2614.03 2092.72 2023.74 1432.07 + 184557 119320 121961 116546 </DataArray> <DataArray type="Float32" Name="rhoW" NumberOfComponents="1" format="ascii"> 1000 1000 1000 1000 </DataArray> - <DataArray type="Float32" Name="rhoN" NumberOfComponents="1" format="ascii"> - 1460 1460 1460 1460 - </DataArray> <DataArray type="Float32" Name="mobW" NumberOfComponents="1" format="ascii"> - 22.587 247.769 312.188 933.954 + 22.5775 249.234 311.805 943.527 </DataArray> - <DataArray type="Float32" Name="mobN" NumberOfComponents="1" format="ascii"> - 989.77 245.533 177.857 2.08447 + <DataArray type="Float32" Name="Sn" NumberOfComponents="1" format="ascii"> + 0.714444 0.347962 0.297603 0.0308204 </DataArray> - <DataArray type="Float32" Name="porosity" NumberOfComponents="1" format="ascii"> - 0.1 0.1 0.3 0.3 + <DataArray type="Float32" Name="pn" NumberOfComponents="1" format="ascii"> + 187171 121411 123985 117961 </DataArray> - <DataArray type="Float32" Name="temperature" NumberOfComponents="1" format="ascii"> - 293.15 293.15 293.15 293.15 + <DataArray type="Float32" Name="rhoN" NumberOfComponents="1" format="ascii"> + 1460 1460 1460 1460 + </DataArray> + <DataArray type="Float32" Name="mobN" NumberOfComponents="1" format="ascii"> + 989.879 243.734 178.199 1.74053 </DataArray> <DataArray type="Float32" Name="process rank" NumberOfComponents="1" format="ascii"> 0 0 0 0