From e01f7122945c3e1b50f8499b25af77cd58694caa Mon Sep 17 00:00:00 2001 From: Arne Morten Kvarving Date: Mon, 16 Sep 2019 10:40:08 +0200 Subject: [PATCH] changed: ewoms/common -> opm/models/utils --- examples/co2injection_flash_ecfv.cpp | 2 +- examples/co2injection_flash_ni_ecfv.cpp | 2 +- examples/co2injection_flash_ni_vcfv.cpp | 2 +- examples/co2injection_flash_vcfv.cpp | 2 +- examples/co2injection_immiscible_ecfv.cpp | 2 +- examples/co2injection_immiscible_ni_ecfv.cpp | 2 +- examples/co2injection_immiscible_ni_vcfv.cpp | 2 +- examples/co2injection_immiscible_vcfv.cpp | 2 +- examples/co2injection_ncp_ecfv.cpp | 2 +- examples/co2injection_ncp_ni_ecfv.cpp | 2 +- examples/co2injection_ncp_ni_vcfv.cpp | 2 +- examples/co2injection_ncp_vcfv.cpp | 2 +- examples/co2injection_pvs_ecfv.cpp | 2 +- examples/co2injection_pvs_ni_ecfv.cpp | 2 +- examples/co2injection_pvs_ni_vcfv.cpp | 2 +- examples/co2injection_pvs_vcfv.cpp | 2 +- examples/cuvette_pvs.cpp | 2 +- examples/diffusion_flash.cpp | 2 +- examples/diffusion_ncp.cpp | 2 +- examples/diffusion_pvs.cpp | 2 +- examples/finger_immiscible_ecfv.cpp | 2 +- examples/fracture_discretefracture.cpp | 2 +- examples/groundwater_immiscible.cpp | 2 +- examples/infiltration_pvs.cpp | 2 +- examples/lens_immiscible_ecfv_ad.cpp | 2 +- examples/lens_immiscible_ecfv_ad_23.cpp | 2 +- examples/lens_immiscible_ecfv_ad_cu1.cpp | 2 +- examples/lens_immiscible_ecfv_ad_cu2.cpp | 2 +- examples/lens_immiscible_vcfv_ad.cpp | 2 +- examples/lens_immiscible_vcfv_fd.cpp | 2 +- examples/lens_richards_ecfv.cpp | 2 +- examples/lens_richards_vcfv.cpp | 2 +- examples/obstacle_immiscible.cpp | 2 +- examples/obstacle_ncp.cpp | 2 +- examples/obstacle_pvs.cpp | 2 +- examples/outflow_pvs.cpp | 2 +- examples/powerinjection_darcy_ad.cpp | 2 +- examples/powerinjection_darcy_fd.cpp | 2 +- examples/powerinjection_forchheimer_ad.cpp | 2 +- examples/powerinjection_forchheimer_fd.cpp | 2 +- examples/reservoir_blackoil_ecfv.cpp | 2 +- examples/reservoir_blackoil_vcfv.cpp | 2 +- examples/reservoir_ncp_ecfv.cpp | 2 +- examples/reservoir_ncp_vcfv.cpp | 2 +- examples/tutorial1.cpp | 2 +- examples/waterair_pvs_ni.cpp | 2 +- opm/models/parallel/threadmanager.hh | 4 +- opm/models/utils/alignedallocator.hh | 226 +++ opm/models/utils/basicproperties.hh | 191 +++ opm/models/utils/genericguard.hh | 93 ++ opm/models/utils/parametersystem.hh | 1202 ++++++++++++++++ opm/models/utils/pffgridvector.hh | 142 ++ opm/models/utils/prefetch.hh | 54 + opm/models/utils/propertysystem.hh | 1209 +++++++++++++++++ opm/models/utils/quadraturegeometries.hh | 154 +++ opm/models/utils/signum.hh | 45 + opm/models/utils/simulator.hh | 957 +++++++++++++ opm/models/utils/start.hh | 443 ++++++ opm/models/utils/timer.hh | 217 +++ opm/models/utils/timerguard.hh | 58 + opm/simulators/linalg/bicgstabsolver.hh | 4 +- .../linalg/istlpreconditionerwrappers.hh | 4 +- opm/simulators/linalg/istlsolverwrappers.hh | 4 +- opm/simulators/linalg/linearsolverreport.hh | 4 +- opm/simulators/linalg/parallelbasebackend.hh | 6 +- opm/simulators/linalg/superlubackend.hh | 2 +- tests/models/test_propertysystem.cpp | 2 +- 67 files changed, 5052 insertions(+), 61 deletions(-) create mode 100644 opm/models/utils/alignedallocator.hh create mode 100644 opm/models/utils/basicproperties.hh create mode 100644 opm/models/utils/genericguard.hh create mode 100644 opm/models/utils/parametersystem.hh create mode 100644 opm/models/utils/pffgridvector.hh create mode 100644 opm/models/utils/prefetch.hh create mode 100644 opm/models/utils/propertysystem.hh create mode 100644 opm/models/utils/quadraturegeometries.hh create mode 100644 opm/models/utils/signum.hh create mode 100644 opm/models/utils/simulator.hh create mode 100644 opm/models/utils/start.hh create mode 100644 opm/models/utils/timer.hh create mode 100644 opm/models/utils/timerguard.hh diff --git a/examples/co2injection_flash_ecfv.cpp b/examples/co2injection_flash_ecfv.cpp index 389cbb76b..8eac6c027 100644 --- a/examples/co2injection_flash_ecfv.cpp +++ b/examples/co2injection_flash_ecfv.cpp @@ -32,7 +32,7 @@ #include #endif -#include +#include #include #include #include "problems/co2injectionflash.hh" diff --git a/examples/co2injection_flash_ni_ecfv.cpp b/examples/co2injection_flash_ni_ecfv.cpp index e8f8b3f86..97ee0ce99 100644 --- a/examples/co2injection_flash_ni_ecfv.cpp +++ b/examples/co2injection_flash_ni_ecfv.cpp @@ -29,7 +29,7 @@ #include "config.h" #include -#include +#include #include #include #include "problems/co2injectionflash.hh" diff --git a/examples/co2injection_flash_ni_vcfv.cpp b/examples/co2injection_flash_ni_vcfv.cpp index 3e628d663..530524fbd 100644 --- a/examples/co2injection_flash_ni_vcfv.cpp +++ b/examples/co2injection_flash_ni_vcfv.cpp @@ -29,7 +29,7 @@ #include "config.h" #include -#include +#include #include #include #include "problems/co2injectionflash.hh" diff --git a/examples/co2injection_flash_vcfv.cpp b/examples/co2injection_flash_vcfv.cpp index 45545d429..8feaae797 100644 --- a/examples/co2injection_flash_vcfv.cpp +++ b/examples/co2injection_flash_vcfv.cpp @@ -32,7 +32,7 @@ #include #endif -#include +#include #include #include #include "problems/co2injectionflash.hh" diff --git a/examples/co2injection_immiscible_ecfv.cpp b/examples/co2injection_immiscible_ecfv.cpp index 9a1c101e8..96306f795 100644 --- a/examples/co2injection_immiscible_ecfv.cpp +++ b/examples/co2injection_immiscible_ecfv.cpp @@ -28,7 +28,7 @@ */ #include "config.h" -#include +#include #include #include diff --git a/examples/co2injection_immiscible_ni_ecfv.cpp b/examples/co2injection_immiscible_ni_ecfv.cpp index 26a3a9183..7aeb3306a 100644 --- a/examples/co2injection_immiscible_ni_ecfv.cpp +++ b/examples/co2injection_immiscible_ni_ecfv.cpp @@ -28,7 +28,7 @@ */ #include "config.h" -#include +#include #include #include #include "problems/co2injectionproblem.hh" diff --git a/examples/co2injection_immiscible_ni_vcfv.cpp b/examples/co2injection_immiscible_ni_vcfv.cpp index 4cdfd8283..a52f5d8db 100644 --- a/examples/co2injection_immiscible_ni_vcfv.cpp +++ b/examples/co2injection_immiscible_ni_vcfv.cpp @@ -28,7 +28,7 @@ */ #include "config.h" -#include +#include #include #include #include "problems/co2injectionproblem.hh" diff --git a/examples/co2injection_immiscible_vcfv.cpp b/examples/co2injection_immiscible_vcfv.cpp index 2ad04c1ec..bfc8cc513 100644 --- a/examples/co2injection_immiscible_vcfv.cpp +++ b/examples/co2injection_immiscible_vcfv.cpp @@ -28,7 +28,7 @@ */ #include "config.h" -#include +#include #include #include diff --git a/examples/co2injection_ncp_ecfv.cpp b/examples/co2injection_ncp_ecfv.cpp index d1c573c3b..5207a29c9 100644 --- a/examples/co2injection_ncp_ecfv.cpp +++ b/examples/co2injection_ncp_ecfv.cpp @@ -28,7 +28,7 @@ */ #include "config.h" -#include +#include #include #include #include "problems/co2injectionproblem.hh" diff --git a/examples/co2injection_ncp_ni_ecfv.cpp b/examples/co2injection_ncp_ni_ecfv.cpp index a5e4b22a0..c346d05f7 100644 --- a/examples/co2injection_ncp_ni_ecfv.cpp +++ b/examples/co2injection_ncp_ni_ecfv.cpp @@ -28,7 +28,7 @@ */ #include "config.h" -#include +#include #include #include #include "problems/co2injectionproblem.hh" diff --git a/examples/co2injection_ncp_ni_vcfv.cpp b/examples/co2injection_ncp_ni_vcfv.cpp index 8d02be162..4516f8813 100644 --- a/examples/co2injection_ncp_ni_vcfv.cpp +++ b/examples/co2injection_ncp_ni_vcfv.cpp @@ -28,7 +28,7 @@ */ #include "config.h" -#include +#include #include #include #include "problems/co2injectionproblem.hh" diff --git a/examples/co2injection_ncp_vcfv.cpp b/examples/co2injection_ncp_vcfv.cpp index 9d947f59a..c0c43e0fe 100644 --- a/examples/co2injection_ncp_vcfv.cpp +++ b/examples/co2injection_ncp_vcfv.cpp @@ -28,7 +28,7 @@ */ #include "config.h" -#include +#include #include #include diff --git a/examples/co2injection_pvs_ecfv.cpp b/examples/co2injection_pvs_ecfv.cpp index 95e8d6dac..c0dd4a1d7 100644 --- a/examples/co2injection_pvs_ecfv.cpp +++ b/examples/co2injection_pvs_ecfv.cpp @@ -28,7 +28,7 @@ */ #include "config.h" -#include +#include #include #include #include "problems/co2injectionproblem.hh" diff --git a/examples/co2injection_pvs_ni_ecfv.cpp b/examples/co2injection_pvs_ni_ecfv.cpp index bb29a4f90..824b104e3 100644 --- a/examples/co2injection_pvs_ni_ecfv.cpp +++ b/examples/co2injection_pvs_ni_ecfv.cpp @@ -28,7 +28,7 @@ */ #include "config.h" -#include +#include #include #include #include "problems/co2injectionproblem.hh" diff --git a/examples/co2injection_pvs_ni_vcfv.cpp b/examples/co2injection_pvs_ni_vcfv.cpp index fdbd23fb8..404502c93 100644 --- a/examples/co2injection_pvs_ni_vcfv.cpp +++ b/examples/co2injection_pvs_ni_vcfv.cpp @@ -28,7 +28,7 @@ */ #include "config.h" -#include +#include #include #include #include "problems/co2injectionproblem.hh" diff --git a/examples/co2injection_pvs_vcfv.cpp b/examples/co2injection_pvs_vcfv.cpp index 8959d8b28..f5fe91cab 100644 --- a/examples/co2injection_pvs_vcfv.cpp +++ b/examples/co2injection_pvs_vcfv.cpp @@ -28,7 +28,7 @@ */ #include "config.h" -#include +#include #include #include diff --git a/examples/cuvette_pvs.cpp b/examples/cuvette_pvs.cpp index cdf52ae6a..fdd88d6f8 100644 --- a/examples/cuvette_pvs.cpp +++ b/examples/cuvette_pvs.cpp @@ -27,7 +27,7 @@ */ #include "config.h" -#include +#include #include #include "problems/cuvetteproblem.hh" diff --git a/examples/diffusion_flash.cpp b/examples/diffusion_flash.cpp index bdcd158d2..aaae2da68 100644 --- a/examples/diffusion_flash.cpp +++ b/examples/diffusion_flash.cpp @@ -27,7 +27,7 @@ */ #include "config.h" -#include +#include #include #include "problems/diffusionproblem.hh" diff --git a/examples/diffusion_ncp.cpp b/examples/diffusion_ncp.cpp index b72fd2657..d91aa3c3c 100644 --- a/examples/diffusion_ncp.cpp +++ b/examples/diffusion_ncp.cpp @@ -27,7 +27,7 @@ */ #include "config.h" -#include +#include #include #include "problems/diffusionproblem.hh" diff --git a/examples/diffusion_pvs.cpp b/examples/diffusion_pvs.cpp index f2f310e98..ebf4096c1 100644 --- a/examples/diffusion_pvs.cpp +++ b/examples/diffusion_pvs.cpp @@ -27,7 +27,7 @@ */ #include "config.h" -#include +#include #include #include "problems/diffusionproblem.hh" diff --git a/examples/finger_immiscible_ecfv.cpp b/examples/finger_immiscible_ecfv.cpp index d6de1ac8b..c79d04464 100644 --- a/examples/finger_immiscible_ecfv.cpp +++ b/examples/finger_immiscible_ecfv.cpp @@ -27,7 +27,7 @@ */ #include "config.h" -#include +#include #include #include #include "problems/fingerproblem.hh" diff --git a/examples/fracture_discretefracture.cpp b/examples/fracture_discretefracture.cpp index 34a81868d..b8094997e 100644 --- a/examples/fracture_discretefracture.cpp +++ b/examples/fracture_discretefracture.cpp @@ -27,7 +27,7 @@ */ #include "config.h" -#include +#include #include "problems/fractureproblem.hh" int main(int argc, char **argv) diff --git a/examples/groundwater_immiscible.cpp b/examples/groundwater_immiscible.cpp index b9f00aa7e..c9a634896 100644 --- a/examples/groundwater_immiscible.cpp +++ b/examples/groundwater_immiscible.cpp @@ -27,7 +27,7 @@ */ #include "config.h" -#include +#include #include #include "problems/groundwaterproblem.hh" diff --git a/examples/infiltration_pvs.cpp b/examples/infiltration_pvs.cpp index 7b7e89763..747bfa14b 100644 --- a/examples/infiltration_pvs.cpp +++ b/examples/infiltration_pvs.cpp @@ -27,7 +27,7 @@ */ #include "config.h" -#include +#include #include #include "problems/infiltrationproblem.hh" diff --git a/examples/lens_immiscible_ecfv_ad.cpp b/examples/lens_immiscible_ecfv_ad.cpp index 3dbc96a57..9a91e1ab7 100644 --- a/examples/lens_immiscible_ecfv_ad.cpp +++ b/examples/lens_immiscible_ecfv_ad.cpp @@ -30,7 +30,7 @@ #include "lens_immiscible_ecfv_ad.hh" -#include +#include int main(int argc, char **argv) { diff --git a/examples/lens_immiscible_ecfv_ad_23.cpp b/examples/lens_immiscible_ecfv_ad_23.cpp index d9266e463..9c79a7587 100644 --- a/examples/lens_immiscible_ecfv_ad_23.cpp +++ b/examples/lens_immiscible_ecfv_ad_23.cpp @@ -82,7 +82,7 @@ public: END_PROPERTIES -#include +#include int main(int argc, char **argv) { diff --git a/examples/lens_immiscible_ecfv_ad_cu1.cpp b/examples/lens_immiscible_ecfv_ad_cu1.cpp index bf1e99341..bcfc7142c 100644 --- a/examples/lens_immiscible_ecfv_ad_cu1.cpp +++ b/examples/lens_immiscible_ecfv_ad_cu1.cpp @@ -36,7 +36,7 @@ #include "lens_immiscible_ecfv_ad.hh" -#include +#include // fake forward declaration to prevent esoteric compiler warning int mainCU1(int argc, char **argv); diff --git a/examples/lens_immiscible_ecfv_ad_cu2.cpp b/examples/lens_immiscible_ecfv_ad_cu2.cpp index 71e9de88b..c6d1d5df8 100644 --- a/examples/lens_immiscible_ecfv_ad_cu2.cpp +++ b/examples/lens_immiscible_ecfv_ad_cu2.cpp @@ -36,7 +36,7 @@ #include "lens_immiscible_ecfv_ad.hh" -#include +#include // fake forward declaration to prevent esoteric compiler warning int mainCU2(int argc, char **argv); diff --git a/examples/lens_immiscible_vcfv_ad.cpp b/examples/lens_immiscible_vcfv_ad.cpp index d51992156..7816ad4c8 100644 --- a/examples/lens_immiscible_vcfv_ad.cpp +++ b/examples/lens_immiscible_vcfv_ad.cpp @@ -28,7 +28,7 @@ */ #include "config.h" -#include +#include #include #include "problems/lensproblem.hh" diff --git a/examples/lens_immiscible_vcfv_fd.cpp b/examples/lens_immiscible_vcfv_fd.cpp index 55f539e00..efde5b6b0 100644 --- a/examples/lens_immiscible_vcfv_fd.cpp +++ b/examples/lens_immiscible_vcfv_fd.cpp @@ -28,7 +28,7 @@ */ #include "config.h" -#include +#include #include #include "problems/lensproblem.hh" diff --git a/examples/lens_richards_ecfv.cpp b/examples/lens_richards_ecfv.cpp index 8ad79058d..979682d5b 100644 --- a/examples/lens_richards_ecfv.cpp +++ b/examples/lens_richards_ecfv.cpp @@ -27,7 +27,7 @@ */ #include "config.h" -#include +#include #include #include "problems/richardslensproblem.hh" diff --git a/examples/lens_richards_vcfv.cpp b/examples/lens_richards_vcfv.cpp index 990dff8d0..c8e7b5828 100644 --- a/examples/lens_richards_vcfv.cpp +++ b/examples/lens_richards_vcfv.cpp @@ -27,7 +27,7 @@ */ #include "config.h" -#include +#include #include #include "problems/richardslensproblem.hh" diff --git a/examples/obstacle_immiscible.cpp b/examples/obstacle_immiscible.cpp index a545a8f76..426cc6e77 100644 --- a/examples/obstacle_immiscible.cpp +++ b/examples/obstacle_immiscible.cpp @@ -28,7 +28,7 @@ */ #include "config.h" -#include +#include #include #include "problems/obstacleproblem.hh" diff --git a/examples/obstacle_ncp.cpp b/examples/obstacle_ncp.cpp index 0c5f50294..e8acc2ef9 100644 --- a/examples/obstacle_ncp.cpp +++ b/examples/obstacle_ncp.cpp @@ -27,7 +27,7 @@ */ #include "config.h" -#include +#include #include #include "problems/obstacleproblem.hh" diff --git a/examples/obstacle_pvs.cpp b/examples/obstacle_pvs.cpp index d35b5a34d..221d787ba 100644 --- a/examples/obstacle_pvs.cpp +++ b/examples/obstacle_pvs.cpp @@ -29,7 +29,7 @@ */ #include "config.h" -#include +#include #include #include "problems/obstacleproblem.hh" diff --git a/examples/outflow_pvs.cpp b/examples/outflow_pvs.cpp index 4c495b4d1..0b8168c53 100644 --- a/examples/outflow_pvs.cpp +++ b/examples/outflow_pvs.cpp @@ -27,7 +27,7 @@ */ #include "config.h" -#include +#include #include #include "problems/outflowproblem.hh" diff --git a/examples/powerinjection_darcy_ad.cpp b/examples/powerinjection_darcy_ad.cpp index 452d62d75..da1b8ecee 100644 --- a/examples/powerinjection_darcy_ad.cpp +++ b/examples/powerinjection_darcy_ad.cpp @@ -27,7 +27,7 @@ */ #include "config.h" -#include +#include #include #include "problems/powerinjectionproblem.hh" diff --git a/examples/powerinjection_darcy_fd.cpp b/examples/powerinjection_darcy_fd.cpp index 65cae1c95..5724e9c15 100644 --- a/examples/powerinjection_darcy_fd.cpp +++ b/examples/powerinjection_darcy_fd.cpp @@ -27,7 +27,7 @@ */ #include "config.h" -#include +#include #include #include "problems/powerinjectionproblem.hh" diff --git a/examples/powerinjection_forchheimer_ad.cpp b/examples/powerinjection_forchheimer_ad.cpp index 6284c548b..c514d3ece 100644 --- a/examples/powerinjection_forchheimer_ad.cpp +++ b/examples/powerinjection_forchheimer_ad.cpp @@ -27,7 +27,7 @@ */ #include "config.h" -#include +#include #include #include "problems/powerinjectionproblem.hh" diff --git a/examples/powerinjection_forchheimer_fd.cpp b/examples/powerinjection_forchheimer_fd.cpp index 56d3a5520..fcebe1d42 100644 --- a/examples/powerinjection_forchheimer_fd.cpp +++ b/examples/powerinjection_forchheimer_fd.cpp @@ -27,7 +27,7 @@ */ #include "config.h" -#include +#include #include #include "problems/powerinjectionproblem.hh" diff --git a/examples/reservoir_blackoil_ecfv.cpp b/examples/reservoir_blackoil_ecfv.cpp index 3832d4351..09fa42890 100644 --- a/examples/reservoir_blackoil_ecfv.cpp +++ b/examples/reservoir_blackoil_ecfv.cpp @@ -28,7 +28,7 @@ */ #include "config.h" -#include +#include #include #include #include "problems/reservoirproblem.hh" diff --git a/examples/reservoir_blackoil_vcfv.cpp b/examples/reservoir_blackoil_vcfv.cpp index d10d5bbe5..ff5e5b996 100644 --- a/examples/reservoir_blackoil_vcfv.cpp +++ b/examples/reservoir_blackoil_vcfv.cpp @@ -27,7 +27,7 @@ */ #include "config.h" -#include +#include #include #include #include "problems/reservoirproblem.hh" diff --git a/examples/reservoir_ncp_ecfv.cpp b/examples/reservoir_ncp_ecfv.cpp index b78af5455..9d9470687 100644 --- a/examples/reservoir_ncp_ecfv.cpp +++ b/examples/reservoir_ncp_ecfv.cpp @@ -27,7 +27,7 @@ */ #include "config.h" -#include +#include #include #include #include "problems/reservoirproblem.hh" diff --git a/examples/reservoir_ncp_vcfv.cpp b/examples/reservoir_ncp_vcfv.cpp index 73aec4e75..0570fc2e0 100644 --- a/examples/reservoir_ncp_vcfv.cpp +++ b/examples/reservoir_ncp_vcfv.cpp @@ -28,7 +28,7 @@ */ #include "config.h" -#include +#include #include #include #include "problems/reservoirproblem.hh" diff --git a/examples/tutorial1.cpp b/examples/tutorial1.cpp index e8509e4ba..a65e78f0f 100644 --- a/examples/tutorial1.cpp +++ b/examples/tutorial1.cpp @@ -27,7 +27,7 @@ * immisciblility. */ #include "config.h" /*@\label{tutorial1:include-begin}@*/ -#include /*@\label{tutorial1:include-end}@*/ +#include /*@\label{tutorial1:include-end}@*/ #include "tutorial1problem.hh" /*@\label{tutorial1:include-problem-header}@*/ int main(int argc, char **argv) diff --git a/examples/waterair_pvs_ni.cpp b/examples/waterair_pvs_ni.cpp index ce77a8e18..52ea06b05 100644 --- a/examples/waterair_pvs_ni.cpp +++ b/examples/waterair_pvs_ni.cpp @@ -27,7 +27,7 @@ */ #include "config.h" -#include +#include #include #include "problems/waterairproblem.hh" diff --git a/opm/models/parallel/threadmanager.hh b/opm/models/parallel/threadmanager.hh index ca7791612..7f89df0c2 100644 --- a/opm/models/parallel/threadmanager.hh +++ b/opm/models/parallel/threadmanager.hh @@ -31,8 +31,8 @@ #include #endif -#include -#include +#include +#include #include diff --git a/opm/models/utils/alignedallocator.hh b/opm/models/utils/alignedallocator.hh new file mode 100644 index 000000000..d5a44d292 --- /dev/null +++ b/opm/models/utils/alignedallocator.hh @@ -0,0 +1,226 @@ +// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- +// vi: set et ts=4 sw=4 sts=4: +/* + 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 2 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 . + + Consult the COPYING file in the top-level source directory of this + module for the precise wording of the license and the list of + copyright holders. +*/ +/*! + * \file + * \brief This is a stand-alone version of boost::alignment::aligned_allocator from Boost + * 1.58 + * + * The file has been modified to assume a C++-2011 compatible compiler on a POSIX + * operating system to remove the boost dependencies which the original version + * contained. The original copyright notice for this file is: + * +
+ (c) 2014 Glen Joseph Fernandes
+ glenjofe at gmail dot com
+
+ Distributed under the Boost Software
+ License, Version 1.0.
+ http://boost.org/LICENSE_1_0.txt
+
+*/ +#ifndef EWOMS_ALIGNED_ALLOCATOR_HH +#define EWOMS_ALIGNED_ALLOCATOR_HH + +#include +#include +#include +#include + +namespace Opm { + +namespace detail { +constexpr inline bool is_alignment(std::size_t value) noexcept +{ + return (value > 0) && ((value & (value - 1)) == 0); +} + +template +struct is_alignment_constant + : std::integral_constant 0) && ((N & (N - 1)) == 0)> +{}; + +template +struct min_size + : std::integral_constant +{ }; + +template +struct offset_object +{ + char offset; + T object; +}; + +template +struct alignment_of + : min_size) - sizeof(T)>::type +{}; + +template +struct max_align + : std::integral_constant B) ? A : B> +{}; + +template +struct max_count_of + : std::integral_constant(0) / sizeof(T)> +{}; + +using std::addressof; +} + +inline void* aligned_alloc(std::size_t alignment, + std::size_t size) noexcept +{ + assert(detail::is_alignment(alignment)); + if (alignment < sizeof(void*)) { + alignment = sizeof(void*); + } + void* p; + if (::posix_memalign(&p, alignment, size) != 0) { + p = 0; + } + return p; +} + +inline void aligned_free(void* ptr) + noexcept +{ + ::free(ptr); +} + + +template +class aligned_allocator { + static_assert(detail::is_alignment_constant::value, "Alignment must be powers of two!"); + +public: + typedef T value_type; + typedef T* pointer; + typedef const T* const_pointer; + typedef void* void_pointer; + typedef const void* const_void_pointer; + typedef std::size_t size_type; + typedef std::ptrdiff_t difference_type; + typedef T& reference; + typedef const T& const_reference; + +private: + typedef detail::max_align::value> MaxAlign; + +public: + template + struct rebind { + typedef aligned_allocator other; + }; + + aligned_allocator() + noexcept = default; + + template + aligned_allocator(const aligned_allocator&) noexcept { + } + + pointer address(reference value) const + noexcept { + return detail::addressof(value); + } + + const_pointer address(const_reference value) const + noexcept { + return detail::addressof(value); + } + + pointer allocate(size_type size, + const_void_pointer = 0) { + void* p = aligned_alloc(MaxAlign::value, + sizeof(T) * size); + if (!p && size > 0) { + throw std::bad_alloc(); + } + return static_cast(p); + } + + void deallocate(pointer ptr, size_type) { + aligned_free(ptr); + } + + constexpr size_type max_size() const + noexcept { + return detail::max_count_of::value; + } + + template + void construct(U* ptr, Args&&... args) { + void* p = ptr; + ::new(p) U(std::forward(args)...); + } + + template + void construct(U* ptr) { + void* p = ptr; + ::new(p) U(); + } + + template + void destroy(U* ptr) { + (void)ptr; + ptr->~U(); + } +}; + +template +class aligned_allocator { + static_assert(detail::is_alignment_constant::value, + "The specified alignment is not a power of two!"); + +public: + typedef void value_type; + typedef void* pointer; + typedef const void* const_pointer; + + template + struct rebind { + typedef aligned_allocator other; + }; +}; + +template +inline bool operator==(const aligned_allocator&, const aligned_allocator&) noexcept +{ + return true; +} + +template +inline bool operator!=(const aligned_allocator&, const aligned_allocator&) noexcept +{ + return false; +} +} + +#endif diff --git a/opm/models/utils/basicproperties.hh b/opm/models/utils/basicproperties.hh new file mode 100644 index 000000000..ab760162a --- /dev/null +++ b/opm/models/utils/basicproperties.hh @@ -0,0 +1,191 @@ +// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- +// vi: set et ts=4 sw=4 sts=4: +/* + 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 2 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 . + + Consult the COPYING file in the top-level source directory of this + module for the precise wording of the license and the list of + copyright holders. +*/ +/*! + * \file + * + * \brief Defines a type tags and some fundamental properties all models. + */ +#ifndef EWOMS_BASIC_PROPERTIES_HH +#define EWOMS_BASIC_PROPERTIES_HH + +#include + +#include +#include +#include + +#if HAVE_DUNE_FEM +#include +#endif + +#include + +BEGIN_PROPERTIES + +/////////////////////////////////// +// Type tag definitions: +// +// NumericModel +// | +// +-> ImplicitModel +/////////////////////////////////// + +//! Type tag for all models. +NEW_TYPE_TAG(NumericModel, INHERITS_FROM(ParameterSystem)); + +//! Type tag for all fully coupled models. +NEW_TYPE_TAG(ImplicitModel, INHERITS_FROM(NumericModel)); + +/////////////////////////////////// +// Property names which are always available: +// +// Scalar +/////////////////////////////////// + +//! Property to specify the type of scalar values. +NEW_PROP_TAG(Scalar); + +//! Property which provides a Dune::ParameterTree. +NEW_PROP_TAG(ParameterTree); + +//! Property which defines the group that is queried for parameters by default +NEW_PROP_TAG(ModelParameterGroup); + +//! Property which provides a Vanguard (manages grids) +NEW_PROP_TAG(Vanguard); + +NEW_PROP_TAG(GridView); + +#if HAVE_DUNE_FEM +NEW_PROP_TAG(GridPart); +#endif + +//! Property which tells the Vanguard how often the grid should be refined +//! after creation. +NEW_PROP_TAG(GridGlobalRefinements); + +//! Property provides the name of the file from which the additional runtime +//! parameters should to be loaded from +NEW_PROP_TAG(ParameterFile); + +/*! + * \brief Print all properties on startup? + * + * 0 means 'no', 1 means 'yes', 2 means 'print only to logfiles'. The + * default is 2. + */ +NEW_PROP_TAG(PrintProperties); + +/*! + * \brief Print all parameters on startup? + * + * 0 means 'no', 1 means 'yes', 2 means 'print only to logfiles'. The + * default is 2. + */ +NEW_PROP_TAG(PrintParameters); + +//! The default value for the simulation's end time +NEW_PROP_TAG(EndTime); + +//! The default value for the simulation's initial time step size +NEW_PROP_TAG(InitialTimeStepSize); + +//! The default value for the simulation's restart time +NEW_PROP_TAG(RestartTime); + +//! The name of the file with a number of forced time step lengths +NEW_PROP_TAG(PredeterminedTimeStepsFile); + +/////////////////////////////////// +// Values for the properties +/////////////////////////////////// + +//! Set the default type of scalar values to double +SET_TYPE_PROP(NumericModel, Scalar, double); + +//! Set the ParameterTree property +SET_PROP(NumericModel, ParameterTree) +{ + typedef Dune::ParameterTree type; + + static Dune::ParameterTree& tree() + { + static Dune::ParameterTree obj_; + return obj_; + } +}; + +//! use the global group as default for the model's parameter group +SET_STRING_PROP(NumericModel, ModelParameterGroup, ""); + +//! Use the DgfVanguard by default +SET_TYPE_PROP(NumericModel, Vanguard, Opm::DgfVanguard); + +//! Set a value for the GridFile property +SET_STRING_PROP(NumericModel, GridFile, ""); + +#if HAVE_DUNE_FEM +SET_PROP(NumericModel, GridPart) +{ + typedef typename GET_PROP_TYPE(TypeTag, Grid) Grid; + typedef Dune::Fem::AdaptiveLeafGridPart type; +}; + +SET_TYPE_PROP(NumericModel, GridView, typename GET_PROP_TYPE(TypeTag, GridPart)::GridViewType); +#else +//! Use the leaf grid view by default. +//! +//! Except for spatial refinement, there is rarly a reason to use +//! anything else... +SET_TYPE_PROP(NumericModel, GridView, typename GET_PROP_TYPE(TypeTag, Grid)::LeafGridView); +#endif + +//! Set a value for the ParameterFile property +SET_STRING_PROP(NumericModel, ParameterFile, ""); + +//! Set the number of refinement levels of the grid to 0. This does not belong +//! here, strictly speaking. +SET_INT_PROP(NumericModel, GridGlobalRefinements, 0); + +//! By default, print the properties on startup +SET_INT_PROP(NumericModel, PrintProperties, 2); + +//! By default, print the values of the run-time parameters on startup +SET_INT_PROP(NumericModel, PrintParameters, 2); + +//! The default value for the simulation's end time +SET_SCALAR_PROP(NumericModel, EndTime, -1e35); + +//! The default value for the simulation's initial time step size +SET_SCALAR_PROP(NumericModel, InitialTimeStepSize, -1e35); + +//! The default value for the simulation's restart time +SET_SCALAR_PROP(NumericModel, RestartTime, -1e35); + +//! By default, do not force any time steps +SET_STRING_PROP(NumericModel, PredeterminedTimeStepsFile, ""); + + +END_PROPERTIES + +#endif diff --git a/opm/models/utils/genericguard.hh b/opm/models/utils/genericguard.hh new file mode 100644 index 000000000..7e9e87983 --- /dev/null +++ b/opm/models/utils/genericguard.hh @@ -0,0 +1,93 @@ +// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- +// vi: set et ts=4 sw=4 sts=4: +/* + 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 2 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 . + + Consult the COPYING file in the top-level source directory of this + module for the precise wording of the license and the list of + copyright holders. +*/ +/*! + * \file + * + * \copydoc Opm::GenericGuard + */ +#ifndef EWOMS_GENERIC_GUARD_HH +#define EWOMS_GENERIC_GUARD_HH + +namespace Opm { +/*! + * \ingroup Common + * + * \brief A simple class which makes sure that a cleanup function is called once the + * object is destroyed. + * + * This class is particularly useful in conjunction with lambdas for code that might + * throw exceptions. + */ +template +class GenericGuard +{ +public: + GenericGuard(Callback& callback) + : callback_(callback) + , isEnabled_(true) + { } + + // allow moves + GenericGuard(GenericGuard&& other) + : callback_(other.callback_) + , isEnabled_(other.isEnabled_) + { + other.isEnabled_ = false; + } + + // disable copies + GenericGuard(const GenericGuard& other) = delete; + + ~GenericGuard() + { + if (isEnabled_) + callback_(); + } + + /*! + * \brief Specify whether the guard object is "on duty" or not. + * + * If the guard object is destroyed while it is "off-duty", the cleanup callback is + * not called. At construction, guards are on duty. + */ + void setEnabled(bool value) + { isEnabled_ = value; } + + /*! + * \brief Returns whether the guard object is "on duty" or not. + */ + bool enabled() const + { return isEnabled_; } + +private: + Callback& callback_; + bool isEnabled_; +}; + +template +GenericGuard make_guard(Callback& callback) +{ return GenericGuard(callback); } + +} // namespace Opm + +#endif diff --git a/opm/models/utils/parametersystem.hh b/opm/models/utils/parametersystem.hh new file mode 100644 index 000000000..a1f5f6af0 --- /dev/null +++ b/opm/models/utils/parametersystem.hh @@ -0,0 +1,1202 @@ +// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- +// vi: set et ts=4 sw=4 sts=4: +/* + 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 2 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 . + + Consult the COPYING file in the top-level source directory of this + module for the precise wording of the license and the list of + copyright holders. +*/ +/*! + * \file + * + * \brief This file provides the infrastructure to retrieve run-time parameters + * + * Internally, runtime parameters are implemented using + * Dune::ParameterTree with the default value taken from the property + * system. + */ +#ifndef EWOMS_PARAMETER_SYSTEM_HH +#define EWOMS_PARAMETER_SYSTEM_HH + +#include + +#include +#include + +#if HAVE_QUAD +#include +#endif // HAVE_QUAD + +#include +#include + +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include + +/*! + * \ingroup Parameter + * + * \brief Register a run-time parameter. + * + * In OPM, parameters can only be used after they have been + * registered. + * + * Example: + * + * \code + * // Registers a run-time parameter "UpwindWeight" which has type + * // "Scalar" and the description "Relative weight of the upwind node." + * EWOMS_REGISTER_PARAM(TypeTag, Scalar, UpwindWeight, + * "Relative weight of the upwind node."); + * \endcode + */ +#define EWOMS_REGISTER_PARAM(TypeTag, ParamType, ParamName, Description) \ + ::Opm::Parameters::registerParam( \ + #ParamName, #ParamName, Description) + +/*! + * \ingroup Parameter + * + * \brief Indicate that a given parameter should not be mentioned in the help message + * + * This allows to deal with unused parameters + */ +#define EWOMS_HIDE_PARAM(TypeTag, ParamName) \ + ::Opm::Parameters::hideParam(#ParamName) + +/*! + * \ingroup Parameter + * + * \brief Indicate that all parameters are registered for a given type tag. + * + * If \c EWOMS_REGISTER_PARAM is called after the invokation of + * \c END_PARAM_REGISTRATION, a std::logic_error exception + * will be thrown. + */ +#define EWOMS_END_PARAM_REGISTRATION(TypeTag) \ + ::Opm::Parameters::endParamRegistration() + +/*! + * \ingroup Parameter + * + * \brief Retrieve a runtime parameter. + * + * The default value is specified via the property system. + * + * Example: + * + * \code + * // Retrieves scalar value UpwindWeight, default + * // is taken from the property UpwindWeight + * EWOMS_GET_PARAM(TypeTag, Scalar, UpwindWeight); + * \endcode + */ +#define EWOMS_GET_PARAM(TypeTag, ParamType, ParamName) \ + (::Opm::Parameters::get(#ParamName, \ + #ParamName)) + +//!\cond SKIP_THIS +#define EWOMS_GET_PARAM_(TypeTag, ParamType, ParamName) \ + (::Opm::Parameters::get( \ + #ParamName, #ParamName, \ + /*errorIfNotRegistered=*/false)) + +/*! + * \ingroup Parameter + * + * \brief Retrieves the lists of parameters specified at runtime and their values. + * + * The two arguments besides the TypeTag are assumed to be STL containers which store + * std::pair. + */ +#define EWOMS_GET_PARAM_LISTS(TypeTag, UsedParamList, UnusedParamList) \ + (::Opm::Parameters::getLists(UsedParamList, UnusedParamList)) + +//!\cond SKIP_THIS +#define EWOMS_RESET_PARAMS_(TypeTag) \ + (::Opm::Parameters::reset()) + +/*! + * \ingroup Parameter + * + * \brief Returns true if a parameter has been specified at runtime, false + * otherwise. + * + * If the parameter in question has not been registered, this throws an exception. + */ +#define EWOMS_PARAM_IS_SET(TypeTag, ParamType, ParamName) \ + (::Opm::Parameters::isSet(#ParamName, \ + #ParamName)) + +namespace Opm { +namespace Parameters { + +struct ParamInfo +{ + std::string paramName; + std::string paramTypeName; + std::string typeTagName; + std::string propertyName; + std::string usageString; + std::string compileTimeValue; + bool isHidden; + + bool operator==(const ParamInfo& other) const + { + return other.paramName == paramName + && other.paramTypeName == paramTypeName + && other.typeTagName == typeTagName + && other.propertyName == propertyName + && other.usageString == usageString + && other.compileTimeValue == compileTimeValue; + } +}; + +// forward declaration +template +const ParamType get(const char *propTagName, + const char *paramName, + bool errorIfNotRegistered = true); + +class ParamRegFinalizerBase_ +{ +public: + virtual ~ParamRegFinalizerBase_() + {} + virtual void retrieve() = 0; +}; + +template +class ParamRegFinalizer_ : public ParamRegFinalizerBase_ +{ +public: + ParamRegFinalizer_(const std::string& paramName) + : paramName_(paramName) + {} + + virtual void retrieve() override + { + // retrieve the parameter once to make sure that its value does + // not contain a syntax error. + ParamType __attribute__((unused)) dummy = + get(/*propTagName=*/paramName_.data(), + paramName_.data(), + /*errorIfNotRegistered=*/true); + } + +private: + std::string paramName_; +}; +} // namespace Parameters + +} // namespace Opm + +BEGIN_PROPERTIES + +// type tag which is supposed to spliced in or inherited from if the +// parameter system is to be used +NEW_TYPE_TAG(ParameterSystem); + +NEW_PROP_TAG(ParameterMetaData); + + +//! Set the ParameterMetaData property +SET_PROP(ParameterSystem, ParameterMetaData) +{ + typedef Dune::ParameterTree type; + + static Dune::ParameterTree& tree() + { return *storage_().tree; } + + static std::map& mutableRegistry() + { return storage_().registry; } + + static const std::map& registry() + { return storage_().registry; } + + static std::list > ®istrationFinalizers() + { return storage_().finalizers; } + + static bool& registrationOpen() + { return storage_().registrationOpen; } + + static void clear() + { + storage_().tree.reset(new Dune::ParameterTree()); + storage_().finalizers.clear(); + storage_().registrationOpen = true; + storage_().registry.clear(); + } + +private: + // this is not pretty, but handling these attributes as static variables inside + // member functions of the ParameterMetaData property class triggers a bug in clang + // 3.5's address sanitizer which causes these variables to be initialized multiple + // times... + struct Storage_ { + Storage_() + { + tree.reset(new Dune::ParameterTree()); + registrationOpen = true; + } + + std::unique_ptr tree; + std::map registry; + std::list > finalizers; + bool registrationOpen; + }; + static Storage_& storage_() { + static Storage_ obj; + return obj; + } +}; + + +END_PROPERTIES + +namespace Opm { + +namespace Parameters { +// function prototype declarations +void printParamUsage_(std::ostream& os, const ParamInfo& paramInfo); +void getFlattenedKeyList_(std::list& dest, + const Dune::ParameterTree& tree, + const std::string& prefix = ""); + +inline std::string breakLines_(const std::string& msg, + int indentWidth, + int maxWidth) +{ + std::string result; + int startInPos = 0; + int inPos = 0; + int lastBreakPos = 0; + int ttyPos = 0; + for (; inPos < int(msg.size()); ++ inPos, ++ ttyPos) { + if (msg[inPos] == '\n') { + result += msg.substr(startInPos, inPos - startInPos + 1); + startInPos = inPos + 1; + lastBreakPos = startInPos + 1; + + // we need to use -1 here because ttyPos is incremented after the loop body + ttyPos = -1; + continue; + } + + if (std::isspace(msg[inPos])) + lastBreakPos = inPos; + + if (ttyPos >= maxWidth) { + if (lastBreakPos > startInPos) { + result += msg.substr(startInPos, lastBreakPos - startInPos); + startInPos = lastBreakPos + 1; + lastBreakPos = startInPos; + inPos = startInPos; + } + else { + result += msg.substr(startInPos, inPos - startInPos); + startInPos = inPos; + lastBreakPos = startInPos; + inPos = startInPos; + } + + result += "\n"; + for (int i = 0; i < indentWidth; ++i) + result += " "; + ttyPos = indentWidth; + } + } + + result += msg.substr(startInPos); + + return result; +} + +inline int getTtyWidth_() +{ + int ttyWidth = 10*1000; // effectively do not break lines at all. + if (isatty(STDOUT_FILENO)) { +#if defined TIOCGWINSZ + // This is a bit too linux specific, IMO. let's do it anyway + struct winsize ttySize; + ioctl(STDOUT_FILENO, TIOCGWINSZ, &ttySize); + ttyWidth = std::max(80, ttySize.ws_col); +#else + // default for systems that do not implement the TIOCGWINSZ ioctl + ttyWidth = 100; +#endif + } + + return ttyWidth; +} + +inline void printParamUsage_(std::ostream& os, const ParamInfo& paramInfo) +{ + std::string paramMessage, paramType, paramDescription; + + int ttyWidth = getTtyWidth_(); + + // convert the CamelCase name to a command line --parameter-name. + std::string cmdLineName = "-"; + const std::string camelCaseName = paramInfo.paramName; + for (unsigned i = 0; i < camelCaseName.size(); ++i) { + if (isupper(camelCaseName[i])) + cmdLineName += "-"; + cmdLineName += static_cast(std::tolower(camelCaseName[i])); + } + + // assemble the printed output + paramMessage = " "; + paramMessage += cmdLineName; + + // add the =VALUE_TYPE part + bool isString = false; + if (paramInfo.paramTypeName == Dune::className() + || paramInfo.paramTypeName == "const char *") + { + paramMessage += "=STRING"; + isString = true; + } + else if (paramInfo.paramTypeName == Dune::className() + || paramInfo.paramTypeName == Dune::className() + || paramInfo.paramTypeName == Dune::className() +#if HAVE_QUAD + || paramInfo.paramTypeName == Dune::className() +#endif // HAVE_QUAD + ) + paramMessage += "=SCALAR"; + else if (paramInfo.paramTypeName == Dune::className() + || paramInfo.paramTypeName == Dune::className() + || paramInfo.paramTypeName == Dune::className() + || paramInfo.paramTypeName == Dune::className()) + paramMessage += "=INTEGER"; + else if (paramInfo.paramTypeName == Dune::className()) + paramMessage += "=BOOLEAN"; + else if (paramInfo.paramTypeName.empty()) { + // the parameter is a flag. Do nothing! + } + else { + // unknown type + paramMessage += "=VALUE"; + } + + // fill up the up help string to the 50th character + paramMessage += " "; + while (paramMessage.size() < 50) + paramMessage += " "; + + + // append the parameter usage string. + paramMessage += paramInfo.usageString; + + // add the default value + if (!paramInfo.paramTypeName.empty()) { + if (paramMessage.back() != '.') + paramMessage += '.'; + paramMessage += " Default: "; + if (paramInfo.paramTypeName == "bool") { + if (paramInfo.compileTimeValue == "0") + paramMessage += "false"; + else + paramMessage += "true"; + } + else if (isString) { + paramMessage += "\""; + paramMessage += paramInfo.compileTimeValue; + paramMessage += "\""; + } + else + paramMessage += paramInfo.compileTimeValue; + } + + paramMessage = breakLines_(paramMessage, /*indent=*/52, ttyWidth); + paramMessage += "\n"; + + // print everything + os << paramMessage; +} + +inline void getFlattenedKeyList_(std::list& dest, + const Dune::ParameterTree& tree, + const std::string& prefix) +{ + // add the keys of the current sub-structure + auto keyIt = tree.getValueKeys().begin(); + const auto& keyEndIt = tree.getValueKeys().end(); + for (; keyIt != keyEndIt; ++keyIt) { + std::string newKey(prefix); + newKey += *keyIt; + dest.push_back(newKey); + } + + // recursively add all substructure keys + auto subStructIt = tree.getSubKeys().begin(); + const auto& subStructEndIt = tree.getSubKeys().end(); + for (; subStructIt != subStructEndIt; ++subStructIt) { + std::string newPrefix(prefix); + newPrefix += *subStructIt; + newPrefix += "."; + + getFlattenedKeyList_(dest, tree.sub(*subStructIt), newPrefix); + } +} + +// print the values of a list of parameters +template +void printParamList_(std::ostream& os, const std::list& keyList, bool printDefaults = false) +{ + typedef typename GET_PROP(TypeTag, ParameterMetaData) ParamsMeta; + + const Dune::ParameterTree& tree = ParamsMeta::tree(); + + auto keyIt = keyList.begin(); + const auto& keyEndIt = keyList.end(); + for (; keyIt != keyEndIt; ++keyIt) { + const auto& paramInfo = ParamsMeta::registry().at(*keyIt); + const std::string& defaultValue = paramInfo.compileTimeValue; + std::string value = defaultValue; + if (tree.hasKey(*keyIt)) + value = tree.get(*keyIt, ""); + os << *keyIt << "=\"" << value << "\""; + if (printDefaults) + os << " # default: \"" << defaultValue << "\""; + os << "\n"; + } +} + +//! \endcond + +/*! + * \ingroup Parameter + * \brief Print a usage message for all run-time parameters. + * + * \param helpPreamble The string that is printed after the error message and before the + * list of parameters. + * \param errorMsg The error message to be printed, if any + * \param os The \c std::ostream which should be used. + */ +template +void printUsage(const std::string& helpPreamble, + const std::string& errorMsg = "", + std::ostream& os = std::cerr) +{ + typedef typename GET_PROP(TypeTag, ParameterMetaData) ParamsMeta; + + if (errorMsg != "") { + os << errorMsg << "\n" + << "\n"; + } + + os << breakLines_(helpPreamble, /*indent=*/2, /*maxWidth=*/getTtyWidth_()); + os << "\n"; + + os << "Recognized options:\n"; + + if (!helpPreamble.empty()) { + ParamInfo pInfo; + pInfo.paramName = "h,--help"; + pInfo.usageString = "Print this help message and exit"; + printParamUsage_(os, pInfo); + } + + auto paramIt = ParamsMeta::registry().begin(); + const auto& paramEndIt = ParamsMeta::registry().end(); + for (; paramIt != paramEndIt; ++paramIt) { + if (!paramIt->second.isHidden) + printParamUsage_(os, paramIt->second); + } +} + +/// \cond 0 +inline int noPositionalParameters_(std::set& seenParams OPM_UNUSED, + std::string& errorMsg, + int argc OPM_UNUSED, + const char** argv, + int paramIdx, + int posParamIdx OPM_UNUSED) +{ + errorMsg = std::string("Illegal parameter \"")+argv[paramIdx]+"\"."; + return 0; +} + +/// \endcond + + +inline void removeLeadingSpace_(std::string& s) +{ + unsigned i; + for (i = 0; i < s.size(); ++ i) + if (!std::isspace(s[i])) + break; + s = s.substr(i); +} + +inline std::string transformKey_(const std::string& s, + bool capitalizeFirstLetter = true, + const std::string& errorPrefix = "") +{ + std::string result; + + if (s.empty()) + throw std::runtime_error(errorPrefix+"Empty parameter names are invalid"); + + if (!std::isalpha(s[0])) + throw std::runtime_error(errorPrefix+"Parameter name '" + s + "' is invalid: First character must be a letter"); + + if (capitalizeFirstLetter) + result += static_cast(std::toupper(s[0])); + else + result += s[0]; + + for (unsigned i = 1; i < s.size(); ++i) { + if (s[i] == '-') { + ++ i; + if (s.size() <= i || !std::isalpha(s[i])) + throw std::runtime_error(errorPrefix+"Invalid parameter name '" + s + "'"); + result += static_cast(std::toupper(s[i])); + } + else if (!std::isalnum(s[i])) + throw std::runtime_error(errorPrefix+"Invalid parameter name '" + s + "'"); + else + result += s[i]; + } + + return result; +} + +inline std::string parseKey_(std::string& s) +{ + unsigned i; + for (i = 0; i < s.size(); ++ i) + if (std::isspace(s[i]) || s[i] == '=') + break; + + std::string ret = s.substr(0, i); + s = s.substr(i); + return ret; +} + +// parse a quoted string +inline std::string parseQuotedValue_(std::string& s, const std::string& errorPrefix) +{ + if (s.empty() || s[0] != '"') + throw std::runtime_error(errorPrefix+"Expected quoted string"); + + std::string result; + unsigned i = 1; + for (; i < s.size(); ++i) { + // handle escape characters + if (s[i] == '\\') { + ++ i; + if (s.size() <= i) + throw std::runtime_error(errorPrefix+"Unexpected end of quoted string"); + + if (s[i] == 'n') + result += '\n'; + else if (s[i] == 'r') + result += '\r'; + else if (s[i] == 't') + result += '\t'; + else if (s[i] == '"') + result += '"'; + else if (s[i] == '\\') + result += '\\'; + else + throw std::runtime_error(errorPrefix+"Unknown escape character '\\" + s[i] + "'"); + } + else if (s[i] == '"') + break; + else + result += s[i]; + } + + s = s.substr(i+1); + return result; +} + +inline std::string parseUnquotedValue_(std::string& s, const std::string& errorPrefix OPM_UNUSED) +{ + unsigned i; + for (i = 0; i < s.size(); ++ i) + if (std::isspace(s[i])) + break; + + std::string ret = s.substr(0, i); + s = s.substr(i); + return ret; +} + +/*! + * \ingroup Parameter + * \brief Parse the parameters provided on the command line. + * + * This function does some basic syntax checks. + * + * \param argc The number of parameters passed by the operating system to the + * main() function + * \param argv The array of strings passed by the operating system to the main() + * function + * \param helpPreamble If non-empty, the --help and -h parameters will be recognized and + * the content of the string will be printed before the list of + * command line parameters + * \return Empty string if everything worked out. Otherwise the thing that could + * not be read. + */ +template +std::string parseCommandLineOptions(int argc, + const char **argv, + const std::string& helpPreamble = "", + const PositionalArgumentCallback& posArgCallback = noPositionalParameters_) +{ + Dune::ParameterTree& paramTree = GET_PROP(TypeTag, ParameterMetaData)::tree(); + + // handle the "--help" parameter + if (!helpPreamble.empty()) { + for (int i = 1; i < argc; ++i) { + if (std::string("-h") == argv[i] + || std::string("--help") == argv[i]) { + printUsage(helpPreamble, /*errorMsg=*/"", std::cout); + return "Help called"; + } + } + } + + std::set seenKeys; + int numPositionalParams = 0; + for (int i = 1; i < argc; ++i) { + // All non-positional command line options need to start with '-' + if (strlen(argv[i]) < 4 + || argv[i][0] != '-' + || argv[i][1] != '-') + { + std::string errorMsg; + int numHandled = posArgCallback(seenKeys, errorMsg, argc, argv, i, numPositionalParams); + + if (numHandled < 1) { + std::ostringstream oss; + + if (!helpPreamble.empty()) + printUsage(helpPreamble, errorMsg, std::cerr); + + return errorMsg; + } + else { + ++ numPositionalParams; + i += numHandled - 1; + continue; + } + } + + std::string paramName, paramValue; + + // read a --my-opt=abc option. This gets transformed + // into the parameter "MyOpt" with the value being + // "abc" + + // There is nothing after the '-' + if (argv[i][2] == 0 || !std::isalpha(argv[i][2])) { + std::ostringstream oss; + oss << "Parameter name of argument " << i + << " ('" << argv[i] << "') " + << "is invalid because it does not start with a letter."; + + if (!helpPreamble.empty()) + printUsage(helpPreamble, oss.str(), std::cerr); + + return oss.str(); + } + + // copy everything after the "--" into a separate string + std::string s(argv[i] + 2); + + // parse argument + paramName = transformKey_(parseKey_(s), /*capitalizeFirst=*/true); + if (seenKeys.count(paramName) > 0) { + std::string msg = + std::string("Parameter '")+paramName+"' specified multiple times as a " + "command line parameter"; + + if (!helpPreamble.empty()) + printUsage(helpPreamble, msg, std::cerr); + return msg; + } + seenKeys.insert(paramName); + + if (s.empty() || s[0] != '=') { + std::string msg = + std::string("Parameter '")+paramName+"' is missing a value. " + +" Please use "+argv[i]+"=value."; + + if (!helpPreamble.empty()) + printUsage(helpPreamble, msg, std::cerr); + return msg; + } + + paramValue = s.substr(1); + + // Put the key=value pair into the parameter tree + paramTree[paramName] = paramValue; + } + return ""; +} + +/*! + * \ingroup Parameter + * \brief Read the parameters from an INI-style file. + * + * This function does some basic syntax checks. + */ +template +void parseParameterFile(const std::string& fileName, bool overwrite = true) +{ + Dune::ParameterTree& paramTree = GET_PROP(TypeTag, ParameterMetaData)::tree(); + + std::set seenKeys; + std::ifstream ifs(fileName); + unsigned curLineNum = 0; + while (ifs) { + // string and file processing in c++ is quite blunt! + std::string curLine; + std::getline(ifs, curLine); + curLineNum += 1; + std::string errorPrefix = fileName+":"+std::to_string(curLineNum)+": "; + + // strip leading white space + removeLeadingSpace_(curLine); + + // ignore empty and comment lines + if (curLine.empty() || curLine[0] == '#' || curLine[0] == ';') + continue; + + // TODO (?): support for parameter groups. + + // find the "key" of the key=value pair + std::string key = parseKey_(curLine); + std::string canonicalKey = transformKey_(key, /*capitalizeFirst=*/true, errorPrefix); + + if (seenKeys.count(canonicalKey) > 0) + throw std::runtime_error(errorPrefix+"Parameter '"+canonicalKey+"' seen multiple times in the same file"); + seenKeys.insert(canonicalKey); + + // deal with the equals sign + removeLeadingSpace_(curLine); + if (curLine.empty() || curLine[0] != '=') + std::runtime_error(errorPrefix+"Syntax error, expecting 'key=value'"); + + curLine = curLine.substr(1); + removeLeadingSpace_(curLine); + + if (curLine.empty() || curLine[0] == '#' || curLine[0] == ';') + std::runtime_error(errorPrefix+"Syntax error, expecting 'key=value'"); + + // get the value + std::string value; + if (curLine[0] == '"') + value = parseQuotedValue_(curLine, errorPrefix); + else + value = parseUnquotedValue_(curLine, errorPrefix); + + // ignore trailing comments + removeLeadingSpace_(curLine); + if (!curLine.empty() && curLine[0] != '#' && curLine[0] != ';') + std::runtime_error(errorPrefix+"Syntax error, expecting 'key=value'"); + + // all went well, add the parameter to the database object + if (overwrite || !paramTree.hasKey(canonicalKey)) + paramTree[canonicalKey] = value; + } +} + +/*! + * \ingroup Parameter + * \brief Print values of the run-time parameters. + * + * \param os The \c std::ostream on which the message should be printed + */ +template +void printValues(std::ostream& os = std::cout) +{ + typedef typename GET_PROP(TypeTag, ParameterMetaData) ParamsMeta; + + const Dune::ParameterTree& tree = ParamsMeta::tree(); + + std::list runTimeAllKeyList; + std::list runTimeKeyList; + std::list unknownKeyList; + + getFlattenedKeyList_(runTimeAllKeyList, tree); + auto keyIt = runTimeAllKeyList.begin(); + const auto& keyEndIt = runTimeAllKeyList.end(); + for (; keyIt != keyEndIt; ++keyIt) { + if (ParamsMeta::registry().find(*keyIt) == ParamsMeta::registry().end()) { + // key was not registered by the program! + unknownKeyList.push_back(*keyIt); + } + else { + // the key was specified at run-time + runTimeKeyList.push_back(*keyIt); + } + } + + // loop over all registered parameters + std::list compileTimeKeyList; + auto paramInfoIt = ParamsMeta::registry().begin(); + const auto& paramInfoEndIt = ParamsMeta::registry().end(); + for (; paramInfoIt != paramInfoEndIt; ++paramInfoIt) { + // check whether the key was specified at run-time + const auto& keyName = paramInfoIt->first; + if (tree.hasKey(keyName)) + continue; + else + compileTimeKeyList.push_back(keyName); + } + + // report the values of all registered (and unregistered) + // parameters + if (runTimeKeyList.size() > 0) { + os << "# [known parameters which were specified at run-time]\n"; + printParamList_(os, runTimeKeyList, /*printDefaults=*/true); + } + + if (compileTimeKeyList.size() > 0) { + os << "# [parameters which were specified at compile-time]\n"; + printParamList_(os, compileTimeKeyList, /*printDefaults=*/false); + } + + if (unknownKeyList.size() > 0) { + os << "# [unused run-time specified parameters]\n"; + auto unusedKeyIt = unknownKeyList.begin(); + const auto& unusedKeyEndIt = unknownKeyList.end(); + for (; unusedKeyIt != unusedKeyEndIt; ++unusedKeyIt) { + os << *unusedKeyIt << "=\"" << tree.get(*unusedKeyIt, "") << "\"\n" << std::flush; + } + } +} + +/*! + * \ingroup Parameter + * \brief Print the list of unused run-time parameters. + * + * \param os The \c std::ostream on which the message should be printed + * + * \return true if something was printed + */ +template +bool printUnused(std::ostream& os = std::cout) +{ + typedef typename GET_PROP(TypeTag, ParameterMetaData) ParamsMeta; + + const Dune::ParameterTree& tree = ParamsMeta::tree(); + std::list runTimeAllKeyList; + std::list unknownKeyList; + + getFlattenedKeyList_(runTimeAllKeyList, tree); + auto keyIt = runTimeAllKeyList.begin(); + const auto& keyEndIt = runTimeAllKeyList.end(); + for (; keyIt != keyEndIt; ++keyIt) { + if (ParamsMeta::registry().find(*keyIt) == ParamsMeta::registry().end()) { + // key was not registered by the program! + unknownKeyList.push_back(*keyIt); + } + } + + if (unknownKeyList.size() > 0) { + os << "# [unused run-time specified parameters]\n"; + auto unusedKeyIt = unknownKeyList.begin(); + const auto& unusedKeyEndIt = unknownKeyList.end(); + for (; unusedKeyIt != unusedKeyEndIt; ++unusedKeyIt) { + os << *unusedKeyIt << "=\"" << tree.get(*unusedKeyIt, "") << "\"\n" << std::flush; + } + return true; + } + return false; +} + +//! \cond SKIP_THIS +template +class Param +{ + typedef typename GET_PROP(TypeTag, ParameterMetaData) ParamsMeta; + +public: + template + static const ParamType get(const char *propTagName, + const char *paramName, + bool errorIfNotRegistered = true) + { + return retrieve_(propTagName, + paramName, + errorIfNotRegistered); + } + + static void clear() + { + ParamsMeta::clear(); + } + + template + static bool isSet(const char *propTagName OPM_OPTIM_UNUSED, + const char *paramName OPM_OPTIM_UNUSED, + bool errorIfNotRegistered = true) + { + +#ifndef NDEBUG + // make sure that the parameter is used consistently. since + // this is potentially quite expensive, it is only done if + // debugging code is not explicitly turned off. + check_(Dune::className(), propTagName, paramName); +#endif + + if (errorIfNotRegistered) { + if (ParamsMeta::registrationOpen()) + throw std::runtime_error("Parameters can only checked after _all_ of them have " + "been registered."); + + if (ParamsMeta::registry().find(paramName) == ParamsMeta::registry().end()) + throw std::runtime_error("Accessing parameter "+std::string(paramName) + +" without prior registration is not allowed."); + } + + std::string canonicalName(paramName); + + // check whether the parameter is in the parameter tree + return ParamsMeta::tree().hasKey(canonicalName); + } + + +private: + struct Blubb + { + std::string propertyName; + std::string paramTypeName; + std::string groupName; + + Blubb& operator=(const Blubb& b) + { + propertyName = b.propertyName; + paramTypeName = b.paramTypeName; + groupName = b.groupName; + return *this; + } + }; + + static void check_(const std::string& paramTypeName, + const std::string& propertyName, + const char *paramName) + { + typedef std::unordered_map StaticData; + static StaticData staticData; + + typename StaticData::iterator it = staticData.find(paramName); + Blubb *b; + if (it == staticData.end()) { + Blubb a; + a.propertyName = propertyName; + a.paramTypeName = paramTypeName; + staticData[paramName] = a; + b = &staticData[paramName]; + } + else + b = &(it->second); + + if (b->propertyName != propertyName) { + throw std::logic_error("GET_*_PARAM for parameter '"+std::string(paramName) + +"' called for at least two different properties ('" + +b->propertyName+"' and '"+propertyName+"')"); + } + + if (b->paramTypeName != paramTypeName) { + throw std::logic_error("GET_*_PARAM for parameter '"+std::string(paramName) + +"' called with at least two different types (" + +b->paramTypeName+" and "+paramTypeName+")"); + } + } + + template + static const ParamType retrieve_(const char OPM_OPTIM_UNUSED *propTagName, + const char *paramName, + bool errorIfNotRegistered = true) + { +#ifndef NDEBUG + // make sure that the parameter is used consistently. since + // this is potentially quite expensive, it is only done if + // debugging code is not explicitly turned off. + check_(Dune::className(), propTagName, paramName); +#endif + + if (errorIfNotRegistered) { + if (ParamsMeta::registrationOpen()) + throw std::runtime_error("Parameters can only retieved after _all_ of them have " + "been registered."); + + if (ParamsMeta::registry().find(paramName) == ParamsMeta::registry().end()) + throw std::runtime_error("Accessing parameter "+std::string(paramName) + +" without prior registration is not allowed."); + } + + // prefix the parameter name by the model's GroupName. E.g. If + // the model specifies its group name to be 'Stokes', in an + // INI file this would result in something like: + // + // [Stokes] + // NewtonWriteConvergence = true + std::string canonicalName(paramName); + + // retrieve actual parameter from the parameter tree + const ParamType defaultValue = GET_PROP_VALUE_(TypeTag, PropTag); + return ParamsMeta::tree().template get(canonicalName, defaultValue ); + } +}; + +template +const ParamType get(const char *propTagName, const char *paramName, bool errorIfNotRegistered) +{ + return Param::template get(propTagName, + paramName, + errorIfNotRegistered); +} + +template +void getLists(Container& usedParams, Container& unusedParams) +{ + usedParams.clear(); + unusedParams.clear(); + + typedef typename GET_PROP(TypeTag, ParameterMetaData) ParamsMeta; + if (ParamsMeta::registrationOpen()) + throw std::runtime_error("Parameter lists can only retieved after _all_ of them have " + "been registered."); + + // get all parameter keys + std::list allKeysList; + const auto& paramTree = ParamsMeta::tree(); + getFlattenedKeyList_(allKeysList, paramTree); + + for (const auto& key : allKeysList) { + if (ParamsMeta::registry().find(key) == ParamsMeta::registry().end()) { + // key was not registered + unusedParams.emplace_back(key, paramTree[key]); + } + else { + // key was registered + usedParams.emplace_back(key, paramTree[key]); + } + } +} + +template +void reset() +{ + return Param::clear(); +} + +template +bool isSet(const char *propTagName, const char *paramName, bool errorIfNotRegistered = true) +{ + return Param::template isSet(propTagName, + paramName, + errorIfNotRegistered); +} + +template +void registerParam(const char *paramName, const char *propertyName, const char *usageString) +{ + typedef typename GET_PROP(TypeTag, ParameterMetaData) ParamsMeta; + if (!ParamsMeta::registrationOpen()) + throw std::logic_error("Parameter registration was already closed before " + "the parameter '"+std::string(paramName)+"' was registered."); + + ParamsMeta::registrationFinalizers().emplace_back( + new ParamRegFinalizer_(paramName)); + + ParamInfo paramInfo; + paramInfo.paramName = paramName; + paramInfo.paramTypeName = Dune::className(); + std::string tmp = Dune::className(); + tmp.replace(0, strlen("Opm::Properties::TTag::"), ""); + paramInfo.propertyName = propertyName; + paramInfo.usageString = usageString; + std::ostringstream oss; + oss << GET_PROP_VALUE_(TypeTag, PropTag); + paramInfo.compileTimeValue = oss.str(); + paramInfo.isHidden = false; + if (ParamsMeta::registry().find(paramName) != ParamsMeta::registry().end()) { + // allow to register a parameter twice, but only if the + // parameter name, type and usage string are exactly the same. + if (ParamsMeta::registry().at(paramName) == paramInfo) + return; + throw std::logic_error("Parameter "+std::string(paramName) + +" registered twice with non-matching characteristics."); + } + + ParamsMeta::mutableRegistry()[paramName] = paramInfo; +} + +template +void hideParam(const char *paramName) +{ + // make sure that a property with the parameter name exists. we cannot check if a + // parameter exists at compile time, so this will only be caught at runtime + static const auto defaultValue OPM_UNUSED = GET_PROP_VALUE_(TypeTag, PropTag); + + typedef typename GET_PROP(TypeTag, ParameterMetaData) ParamsMeta; + if (!ParamsMeta::registrationOpen()) + throw std::logic_error("Parameter '"+std::string(paramName)+"' declared as hidden" + " when parameter registration was already closed."); + + auto paramInfoIt = ParamsMeta::mutableRegistry().find(paramName); + if (paramInfoIt == ParamsMeta::mutableRegistry().end()) + throw std::logic_error("Tried to declare unknown parameter '" + +std::string(paramName)+"' hidden."); + + auto& paramInfo = paramInfoIt->second; + paramInfo.isHidden = true; +} + +template +void endParamRegistration() +{ + typedef typename GET_PROP(TypeTag, ParameterMetaData) ParamsMeta; + if (!ParamsMeta::registrationOpen()) + throw std::logic_error("Parameter registration was already closed. It is only possible " + "to close it once."); + + ParamsMeta::registrationOpen() = false; + + // loop over all parameters and retrieve their values to make sure + // that there is no syntax error + auto pIt = ParamsMeta::registrationFinalizers().begin(); + const auto& pEndIt = ParamsMeta::registrationFinalizers().end(); + for (; pIt != pEndIt; ++pIt) + (*pIt)->retrieve(); + ParamsMeta::registrationFinalizers().clear(); +} +//! \endcond + +} // namespace Parameters +} // namespace Opm + +#endif diff --git a/opm/models/utils/pffgridvector.hh b/opm/models/utils/pffgridvector.hh new file mode 100644 index 000000000..46e65757d --- /dev/null +++ b/opm/models/utils/pffgridvector.hh @@ -0,0 +1,142 @@ +// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- +// vi: set et ts=4 sw=4 sts=4: +/* + 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 2 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 . + + Consult the COPYING file in the top-level source directory of this + module for the precise wording of the license and the list of + copyright holders. +*/ +/*! + * \file + * \copydoc PffGridVector + */ +#ifndef EWOMS_PFF_GRID_VECTOR_HH +#define EWOMS_PFF_GRID_VECTOR_HH + +#include + +#include +#include + +#include + +namespace Opm { +/*! + * \brief A random-access container which stores data attached to a grid's degrees of + * freedom in a prefetch friendly manner. + * + * This container often reduces the number of cache faults considerably, thus improving + * performance. On the flipside data cannot be written to on an individual basis and it + * requires significantly more memory than a plain array. PffVector stands for "PreFetch + * Friendly Grid Vector". + */ +template +class PffGridVector +{ + typedef typename GridView::template Codim<0>::Entity Element; + +#if DUNE_VERSION_NEWER(DUNE_GRID, 2,6) + typedef Dune::MultipleCodimMultipleGeomTypeMapper ElementMapper; +#else + typedef Dune::MultipleCodimMultipleGeomTypeMapper ElementMapper; +#endif + +public: + PffGridVector(const GridView& gridView, const DofMapper& dofMapper) + : gridView_(gridView) +#if DUNE_VERSION_NEWER(DUNE_GRID, 2,6) + , elementMapper_(gridView_, Dune::mcmgElementLayout()) +#else + , elementMapper_(gridView_) +#endif + , dofMapper_(dofMapper) + { } + + template + void update(const DistFn& distFn) + { + unsigned numElements = gridView_.size(/*codim=*/0); + unsigned numLocalDofs = computeNumLocalDofs_(); + + elemData_.resize(numElements); + data_.resize(numLocalDofs); + + // update the pointers for the element data: for this, we need to loop over the + // whole grid and update a stencil for each element + Data *curElemDataPtr = &data_[0]; + Stencil stencil(gridView_, dofMapper_); + auto elemIt = gridView_.template begin(); + const auto& elemEndIt = gridView_.template end(); + for (; elemIt != elemEndIt; ++elemIt) { + // set the DOF data pointer for the current element + const auto& elem = *elemIt; + unsigned elemIdx = elementMapper_.index(elem); + elemData_[elemIdx] = curElemDataPtr; + + stencil.update(elem); + unsigned numDof = stencil.numDof(); + for (unsigned localDofIdx = 0; localDofIdx < numDof; ++ localDofIdx) + distFn(curElemDataPtr[localDofIdx], stencil, localDofIdx); + + // update the element data pointer to make it point to the beginning of the + // data for DOFs of the next element + curElemDataPtr += numDof; + } + } + + void prefetch(const Element& elem) const + { + unsigned elemIdx = elementMapper_.index(elem); + + // we use 0 as the temporal locality, because it is reasonable to assume that an + // entry will only be accessed once. + Opm::prefetch(elemData_[elemIdx]); + } + + const Data& get(const Element& elem, unsigned localDofIdx) const + { + unsigned elemIdx = elementMapper_.index(elem); + return elemData_[elemIdx][localDofIdx]; + } + +private: + unsigned computeNumLocalDofs_() const + { + unsigned result = 0; + + // loop over the whole grid and sum up the number of local DOFs of all Stencils + Stencil stencil(gridView_, dofMapper_); + auto elemIt = gridView_.template begin(); + const auto& elemEndIt = gridView_.template end(); + for (; elemIt != elemEndIt; ++elemIt) { + stencil.update(*elemIt); + result += stencil.numDof(); + } + + return result; + } + + GridView gridView_; + ElementMapper elementMapper_; + const DofMapper& dofMapper_; + std::vector data_; + std::vector elemData_; +}; + +} // namespace Opm + +#endif diff --git a/opm/models/utils/prefetch.hh b/opm/models/utils/prefetch.hh new file mode 100644 index 000000000..5f2c77552 --- /dev/null +++ b/opm/models/utils/prefetch.hh @@ -0,0 +1,54 @@ +// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- +// vi: set et ts=4 sw=4 sts=4: +/* + 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 2 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 . + + Consult the COPYING file in the top-level source directory of this + module for the precise wording of the license and the list of + copyright holders. +*/ +/*! + * \file + * \copydoc prefetch + */ +#ifndef EWOMS_PREFETCH_HH +#define EWOMS_PREFETCH_HH + +namespace Opm { +/*! + * \brief Template function which emits prefetch instructions for a range of memory + * + * This function does not change the semantics of the code, but used correctly it will + * improve performace because the number of cache misses will be reduced. + */ +template +void prefetch(const T& val, unsigned n = 1) +{ +#if __clang__ || __GNUC__ + // this value is architecture specific, but a cache line size of 64 bytes seems to be + // used by all contemporary architectures. + static const int cacheLineSize = 64; + + const char *beginPtr = reinterpret_cast(&val); + const char *endPtr = reinterpret_cast(&val + n); + for (; beginPtr < endPtr; beginPtr += cacheLineSize) + __builtin_prefetch(beginPtr, writeOnly, temporalLocality); +#endif +} + +} // namespace Opm + +#endif diff --git a/opm/models/utils/propertysystem.hh b/opm/models/utils/propertysystem.hh new file mode 100644 index 000000000..a6a3368ed --- /dev/null +++ b/opm/models/utils/propertysystem.hh @@ -0,0 +1,1209 @@ +// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- +// vi: set et ts=4 sw=4 sts=4: +/* + 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 2 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 . + + Consult the COPYING file in the top-level source directory of this + module for the precise wording of the license and the list of + copyright holders. +*/ +/*! + * \file + * \brief Provides the magic behind the eWoms property system. + * + * Properties allow to associate arbitrary data types to + * identifiers. A property is always defined on a pair (\c TypeTag, + * \c PropertyTag) where \c TypeTag is the identifier for the object the + * property is defined for and \c PropertyTag is an unique identifier of + * the property. + * + * Type tags are hierarchic and inherit properties defined on their + * ancesters. At each level, properties defined on lower levels can be + * overwritten or even made undefined. + * + * Properties may use other properties for the respective type tag and + * these properties can also be defined on an arbitrary level of the + * hierarchy. The only restriction on this is that cycles are not + * allowed when defining properties. + */ +#ifndef EWOMS_PROPERTIES_HH +#define EWOMS_PROPERTIES_HH + +#include +#include + +#include + +#include // required for 'is_base_of' +#include +#include +#include +#include +#include +#include +#include +#include +#include + +//! \cond SKIP_THIS + +namespace Opm { +namespace Properties { + +#define EWOMS_GET_HEAD_(Arg1, ...) Arg1 + +#if !defined NO_PROPERTY_INTROSPECTION + +//! Internal macro which is only required if the property introspection is enabled +#define PROP_INFO_(EffTypeTagName, PropKind, PropTagName, ...) \ + template <> \ + struct PropertyInfo \ + { \ + static int init() { \ + PropertyRegistryKey key( \ + /*effTypeTagName=*/ Dune::className(), \ + /*kind=*/PropKind, \ + /*name=*/propertyName(), \ + /*value=*/#__VA_ARGS__, \ + /*file=*/__FILE__, \ + /*line=*/__LINE__); \ + PropertyRegistry::addKey(key); \ + return 0; \ + } \ + static std::string propertyName() { return #PropTagName; } \ + }; \ + namespace fooPropInfo_ ## EffTypeTagName { \ + static const int foo_ ## PropTagName OPM_UNUSED = \ + PropertyInfo::init(); \ + } + + +//! Internal macro which is only required if the property introspection is enabled +#define TTAG_INFO_(TagName, ...) \ + template <> \ + struct TypeTagInfo \ + { \ + static int init() { \ + TypeTagRegistry::addAllChildren(); \ + return 0; \ + } \ + }; \ + static const int fooTypeTagInfo_ ## TagName OPM_UNUSED = \ + TypeTagInfo::init(); + +//! Internal macro which is only required if the property introspection is enabled +#define SPLICE_INFO_(SpliceName, ...) \ + template <> \ + struct SpliceInfo \ + { \ + static int init() { \ + TypeTagRegistry::addAllSplices(); \ + return 0; \ + } \ + }; \ + static const int fooSpliceInfo_ ## SpliceName OPM_UNUSED = \ + SpliceInfo::init(); + +#else +//! Don't do anything if introspection is disabled +#define PROP_INFO_(EffTypeTagName, PropKind, PropTagName, ...) +#define TTAG_INFO_(EffTypeTagName, ...) +#define SPLICE_INFO_(EffTypeTagName, ...) +#endif + +// some macros for simplification + +//! \endcond + +/*! + * \ingroup Properties + * \brief Indicates that property definitions follow + */ +#define BEGIN_PROPERTIES namespace Opm { namespace Properties { + +/*! + * \ingroup Properties + * \brief Indicates that all properties have been specified (for now) + */ +#define END_PROPERTIES }} + +/*! + * \ingroup Properties + * \brief Convert a type tag name to a type + * + * The main advantage of the type of a \c TypeTag is that it can be + * passed as a template argument. + */ +#define TTAG(TypeTagName) Opm::Properties::TTag::TypeTagName + +/*! + * \ingroup Properties + * \brief Makes a type out of a property tag name + * + * Again property type names can be passed as template argument. This + * is rarely needed, though. + */ +#define PTAG(PropTagName) Opm::Properties::PTag::PropTagName + +/*! + * \ingroup Properties + * \brief Define a new type tag. + * + * A type tag can inherit the properties defined on up to five parent + * type tags. Examples: + * + * \code + * // The type tag doesn't inherit any properties from other type tags + * NEW_TYPE_TAG(FooTypeTag); + * + * // BarTypeTag inherits all properties from FooTypeTag + * NEW_TYPE_TAG(BarTypeTag, INHERITS_FROM(FooTypeTag)); + * + * // FooBarTypeTag inherits the properties of FooTypeTag as well as + * // those of BarTypeTag. Properties defined on BarTypeTag have + * // preceedence over those defined for FooTypeTag: + * NEW_TYPE_TAG(FooBarTypeTag, INHERITS_FROM(FooTypeTag, BarTypeTag)); + * \endcode + */ +#define NEW_TYPE_TAG(...) \ + namespace TTag { \ + struct EWOMS_GET_HEAD_(__VA_ARGS__, dontcare) \ + : public TypeTag<__VA_ARGS__> \ + { }; \ + TTAG_INFO_(__VA_ARGS__, void) \ + } \ + extern int semicolonHack_ + +/*! + * \ingroup Properties + * \brief Define splices for a given type tag. + * + * Splices can be seen as children which can be overridden lower in + * the hierarchy. It can thus be seen as a "deferred inheritance" + * mechanism. Example: + * + * \code + * // First, define type tags for two different linear solvers: + * // BiCGStab and SuperLU. The first needs the "MaxIterations" + * // property, the second defines the "UsePivoting" property. + * NEW_TYPE_TAG(BiCGStabSolver); + * NEW_PROP_TAG(MaxIterations); + * SET_INT_PROP(BiCGStabSolver, MaxIterations, 100); + * + * NEW_TYPE_TAG(SuperLUSolver); + * NEW_PROP_TAG(UsePivoting); + * SET_BOOL_PROP(SuperLUSolver, UsePivoting, true); + * + * // The model type tag defines the splice 'LinearSolver' and sets it + * // to the 'BiCGStabSolver' type tag. + * NEW_TYPE_TAG(ModelTypeTag); + * NEW_PROP_TAG(LinearSolver); + * SET_SPLICES(ModelTypeTag, LinearSolver); + * SET_TAG_PROP(ModelTypeTag, LinearSolver, BiCGStabSolver); + * + * // The problem type tag is derived from the model type tag, but uses + * // the SuperLU solver. Since this is done using a splice, all properties + * // defined for the "SuperLUSolver" are inherited and the ones for the + * // BiCGStabSolver type tag are undefined + * NEW_TYPE_TAG(ProblemTypeTag, INHERITS_FROM(ModelTypeTag)); + * SET_TAG_PROP(ProblemTypeTag, LinearSolver, SuperLUSolver); + * \endcode + */ +#define SET_SPLICES(TypeTagName, ...) \ + namespace PTag { \ + template<> \ + struct Splices \ + { \ + typedef RevertedTuple<__VA_ARGS__>::type tuple; \ + }; \ + SPLICE_INFO_(TypeTagName, __VA_ARGS__) \ + } \ + extern int semicolonHack_ + +/*! + * \ingroup Properties + * \brief Syntactic sugar for NEW_TYPE_TAG. + * + * See the documentation for NEW_TYPE_TAG. + */ +#define INHERITS_FROM(...) __VA_ARGS__ + +/*! +// * \ingroup Properties + * \brief Define a property tag. + * + * A property tag is the unique identifier for a property. It may only + * be declared once in your program. There is also no hierarchy of + * property tags as for type tags. + * + * Examples: + * + * \code + * NEW_PROP_TAG(blubbPropTag); + * NEW_PROP_TAG(blabbPropTag); + * \endcode + */ +#define NEW_PROP_TAG(PTagName) \ + namespace PTag { \ + struct PTagName; } extern int semicolonHack_ + +//! \cond SKIP_THIS +#define SET_PROP_(EffTypeTagName, PropKind, PropTagName, ...) \ + template \ + struct Property; \ + PROP_INFO_(EffTypeTagName, \ + /*kind=*/PropKind, \ + PropTagName, \ + /*value=*/__VA_ARGS__) \ + template \ + struct Property +//! \endcond + +/*! + * \ingroup Properties + * \brief Set a property for a specific type tag. + * + * After this macro, you must to specify a complete body of a class + * template, including the trailing semicolon. If you need to retrieve + * another property within the class body, you can use \c TypeTag as the + * argument for the type tag for the \c GET_PROP macro. + * + * Example: + * + * \code + * SET_PROP(FooTypeTag, blubbPropTag) + * { + * static int value = 10; + * static int calculate(int arg) + * { calculateInternal_(arg); } + * + * private: + * // retrieve the blabbProp property for the TypeTag the + * // property is defined on. Note that blabbProb does not need to + * // be defined on FooTypeTag, but can also be defined for some + * // derived type tag. + * typedef typename GET_PROP(TypeTag, blabbProp) blabb; + * + * static int calculateInternal_(int arg) + * { return arg * blabb::value; }; + * \endcode + * }; + */ +#define SET_PROP(EffTypeTagName, PropTagName) \ + template \ + struct Property; \ + PROP_INFO_(EffTypeTagName, \ + /*kind=*/"opaque", \ + PropTagName, \ + /*value=*/"") \ + template \ + struct Property + +/*! + * \ingroup Properties + * \brief Explicitly unset a property for a type tag. + * + * This means that the property will not be inherited from the type + * tag's parents. + * + * Example: + * + * \code + * // make the blabbPropTag property undefined for the BarTypeTag. + * UNSET_PROP(BarTypeTag, blabbPropTag); + * \endcode + */ +#define UNSET_PROP(EffTypeTagName, PropTagName) \ + template <> \ + struct PropertyUnset; \ + PROP_INFO_(EffTypeTagName, \ + /*kind=*/"withdraw", \ + PropTagName, \ + /*value=*/) \ + template <> \ + struct PropertyUnset \ + : public PropertyExplicitlyUnset \ + {} + +/*! + * \ingroup Properties + * \brief Set a property to a simple constant integer value. + * + * The constant can be accessed by the \c value attribute. + */ +#define SET_INT_PROP(EffTypeTagName, PropTagName, /*Value*/...) \ + SET_PROP_(EffTypeTagName, \ + /*kind=*/"int ", \ + PropTagName, \ + /*value=*/__VA_ARGS__) \ + { \ + typedef int type; \ + static const int value = __VA_ARGS__; \ + } + +/*! + * \ingroup Properties + * \brief Set a property to a simple constant boolean value. + * + * The constant can be accessed by the \c value attribute. + */ +#define SET_BOOL_PROP(EffTypeTagName, PropTagName, /*Value*/...) \ + SET_PROP_(EffTypeTagName, \ + /*kind=*/"bool ", \ + PropTagName, \ + /*value=*/__VA_ARGS__) \ + { \ + typedef bool type; \ + static const bool value = __VA_ARGS__; \ + } + +/*! + * \ingroup Properties + * \brief Set a property which defines a type. + * + * The type can be accessed by the \c type attribute. + */ +#define SET_TYPE_PROP(EffTypeTagName, PropTagName, /*Value*/...) \ + SET_PROP_(EffTypeTagName, \ + /*kind=*/"type ", \ + PropTagName, \ + /*value=*/__VA_ARGS__) \ + { \ + typedef __VA_ARGS__ type; \ + } + +/*! + * \ingroup Properties + * \brief This macro provides the boiler plate code to declare a static constant member + * outside of a property definition. + * + * This macro is fairly low-level and there should rarely be a need to use it. + */ +#define PROP_STATIC_CONST_MEMBER_DEFINITION_PREFIX_(EffTypeTagName, PropTagName) \ + template \ + const typename Property::type \ + Property + +/*! + * \ingroup Properties + * \brief This macro provides the boiler plate code to declare a static member outside of + * a property definition. + * + * This macro is fairly low-level and there should rarely be a need to use it. + */ +#define PROP_STATIC_MEMBER_DEFINITION_PREFIX_(EffTypeTagName, PropTagName) \ + template \ + typename Property::type \ + Property + +/*! + * \ingroup Properties + * \brief Set a property to a simple constant scalar value. + * + * The constant can be accessed by the \c value attribute. In order to + * use this macro, the property tag \c Scalar needs to be defined for + * the type tag. + */ +#define SET_SCALAR_PROP(EffTypeTagName, PropTagName, ...) \ + SET_PROP_(EffTypeTagName, \ + /*kind=*/"scalar", \ + PropTagName, \ + /*value=*/__VA_ARGS__) \ + { \ + typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar; \ + public: \ + typedef Scalar type; \ + static const Scalar value; \ + }; \ + PROP_STATIC_CONST_MEMBER_DEFINITION_PREFIX_(EffTypeTagName, PropTagName)::value(__VA_ARGS__) + +/*! + * \ingroup Properties + * \brief Set a property to a simple constant string value. + * + * The constant can be accessed by the \c value attribute and is of + * type std::string. + */ +#define SET_STRING_PROP(EffTypeTagName, PropTagName, ...) \ + SET_PROP_(EffTypeTagName, \ + /*kind=*/"string", \ + PropTagName, \ + /*value=*/__VA_ARGS__) \ + { \ + public: \ + typedef std::string type; \ + static const std::string value; \ + }; \ + PROP_STATIC_CONST_MEMBER_DEFINITION_PREFIX_(EffTypeTagName, PropTagName)::value(__VA_ARGS__) + +/*! + * \ingroup Properties + * \brief Define a property containing a type tag. + * + * This is convenient for splices. + */ +#define SET_TAG_PROP(EffTypeTagName, PropTagName, ValueTypeTagName) \ + SET_PROP_(EffTypeTagName, \ + /*kind=*/"tag ", \ + PropTagName, \ + /*value=*/TTAG(ValueTypeTagName)) \ + { \ + typedef TTAG(ValueTypeTagName) type; \ + } + + +/*! + * \ingroup Properties + * \brief Retrieve a property for a type tag. + * + * If you use \c GET_PROP within a template and want to refer to some + * type (including the property itself), \c GET_PROP must be preceeded by + * the '\c typename' keyword. + */ +#define GET_PROP(TypeTag, PropTagName) \ + ::Opm::Properties::GetProperty::p +//!\cond SKIP_THIS +#define GET_PROP_(TypeTag, PropTag) \ + ::Opm::Properties::GetProperty::p +//!\endcond + +/*! + * \ingroup Properties + * \brief Access the \c value attribute of a property for a type tag. + * + * This is just for convenience and equivalent to + * GET_PROP(TypeTag, PropTag)::value . If the property doesn't + * have an attribute named \c value, this yields a compiler error. + */ +#define GET_PROP_VALUE(TypeTag, PropTagName) \ + ::Opm::Properties::GetProperty::p::value +//!\cond SKIP_THIS +#define GET_PROP_VALUE_(TypeTag, PropTag) \ + ::Opm::Properties::GetProperty::p::value +//!\endcond + +/*! + * \ingroup Properties + * \brief Access the \c type attribute of a property for a type tag. + * + * This is just for convenience and equivalent to + * GET_PROP(TypeTag, PropTag)::type. If the property doesn't + * have an attribute named \c type, this yields a compiler error. Also, + * if you use this macro within a template, it must be preceeded by + * the \c typename keyword. + */ +#define GET_PROP_TYPE(TypeTag, PropTagName) \ + ::Opm::Properties::GetProperty::p::type +//!\cond SKIP_THIS +#define GET_PROP_TYPE_(TypeTag, PropTag) \ + ::Opm::Properties::GetProperty::p::type +//!\endcond + +#if !defined NO_PROPERTY_INTROSPECTION +/*! + * \ingroup Properties + * \brief Return a human readable diagnostic message how exactly a + * property was defined. + * + * This is only enabled if the \c NO_PROPERTY_INTROSPECTION macro is not + * defined. + * + * Example: + * + * \code + * int main() + * { + * std::cout << PROP_DIAGNOSTIC(FooBarTypeTag, blabbPropTag) << "\n"; + * }; + * \endcode + */ +#define PROP_DIAGNOSTIC(TypeTag, PropTagName) \ + ::Opm::Properties::getDiagnostic(#PropTagName) + +#else +#define PROP_DIAGNOSTIC(TypeTag, PropTagName) "Property introspection disabled by macro NO_PROPERTY_INTROSPECTION." +#endif + +////////////////////////////////////////////// +// some serious template kung fu. Don't look at it too closely, it +// might damage your brain! +////////////////////////////////////////////// + +//! \cond SKIP_THIS + +namespace PTag {} +namespace TTag {} + +#if !defined NO_PROPERTY_INTROSPECTION + +namespace TTag +{ +template +struct TypeTagInfo +{}; +} + +namespace PTag +{ +template +struct SpliceInfo +{}; +} + +template +struct PropertyInfo +{}; + +class PropertyRegistryKey +{ +public: + PropertyRegistryKey() + {} + + PropertyRegistryKey(const std::string& effTypeTagName, + const std::string& propertyKind, + const std::string& propertyName, + const std::string& propertyValue, + const std::string& fileDefined, + int lineDefined) + : effTypeTagName_(effTypeTagName) + , propertyKind_(propertyKind) + , propertyName_(propertyName) + , propertyValue_(propertyValue) + , fileDefined_(fileDefined) + , lineDefined_(lineDefined) + { } + + // copy constructor + PropertyRegistryKey(const PropertyRegistryKey&) = default; + PropertyRegistryKey& operator=(const PropertyRegistryKey&) = default; + + const std::string& effTypeTagName() const + { return effTypeTagName_; } + const std::string& propertyKind() const + { return propertyKind_; } + const std::string& propertyName() const + { return propertyName_; } + const std::string& propertyValue() const + { return propertyValue_; } + const std::string& fileDefined() const + { return fileDefined_; } + int lineDefined() const + { return lineDefined_; } + +private: + std::string effTypeTagName_; + std::string propertyKind_; + std::string propertyName_; + std::string propertyValue_; + std::string fileDefined_; + int lineDefined_; +}; + + +template +struct GetProperty; + +class TypeTagRegistry +{ +public: + struct SpliceRegistryEntry { + SpliceRegistryEntry(const std::string& name) + { name_ = name; } + + std::string propertyName() const + { return name_; } + + private: + std::string name_; + }; + + typedef std::list > SpliceList; + typedef std::map SpliceListMap; + + typedef std::list ChildrenList; + typedef std::map ChildrenListMap; + + // this method adds the children of a type tag to the registry if the registry does + // not yet contain an entry for that type tag + template + static void addAllChildren() + { + std::string typeTagName = Dune::className(); + + if (children_().find(typeTagName) == children_().end()) + addChildren(); + } + + + // end of recursion. the last argument is not a child, but 'void' + // which is required for the macro magic... + template + static void addChildren() + {} + + // the last argument is not a child, but 'void' which is required + // for the macro magic... + template + static void addChildren() + { + std::string typeTagName = Dune::className(); + + children_()[typeTagName].emplace_front(Dune::className()); + + addChildren(); + } + + // this method adds the splices of a type tag to the registry if the registry does + // not yet contain an entry for that type tag + template + static void addAllSplices() + { + std::string typeTagName = Dune::className(); + + if (splices_().find(typeTagName) == splices_().end()) + addSplices(); + } + + + // end of recursion. the last argument is not a child, but 'void' + // which is required for the macro magic... + template + static void addSplices() + { } + + // the last argument is not a child, but 'void' which is required + // for the macro magic... + template + static void addSplices() + { + const std::string& typeTagName = Dune::className(); + const std::string& propName = + PropertyInfo::template GetEffectiveTypeTag_::type, Splice1>::propertyName(); + + splices_()[typeTagName].emplace_front(new SpliceRegistryEntry(propName)); + + addSplices(); + } + + static const SpliceList& splices(const std::string& typeTagName) + { return splices_()[typeTagName]; } + + static const ChildrenList& children(const std::string& typeTagName) + { return children_()[typeTagName]; } + +private: + static SpliceListMap& splices_() + { + static SpliceListMap data; + return data; + } + static ChildrenListMap& children_() + { + static ChildrenListMap data; + return data; + } +}; + +class PropertyRegistry +{ + typedef Opm::Properties::TypeTagRegistry TypeTagRegistry; + +public: + typedef std::map KeyList; + typedef std::map KeyListMap; + + static void addKey(const PropertyRegistryKey& key) + { + keys_()[key.effTypeTagName()][key.propertyName()] = key; + } + + static const std::string& getSpliceTypeTagName(const std::string& typeTagName, + const std::string& propertyName) + { + const auto& keyIt = keys_().find(typeTagName); + const auto& keyEndIt = keys_().end(); + if (keyIt == keyEndIt) + throw std::runtime_error("Unknown type tag key '"+typeTagName+"'"); + + // check whether the propery is directly defined for the type tag currently + // checked, ... + const auto& propIt = keyIt->second.find(propertyName); + const auto& propEndIt = keyIt->second.end(); + if (propIt != propEndIt) + return propIt->second.propertyValue(); + + // ..., if not, check all splices, ... + typedef typename TypeTagRegistry::SpliceList SpliceList; + const SpliceList& splices = TypeTagRegistry::splices(typeTagName); + SpliceList::const_iterator spliceIt = splices.begin(); + for (; spliceIt != splices.end(); ++spliceIt) { + const auto& spliceTypeTagName = + PropertyRegistry::getSpliceTypeTagName(typeTagName, + (*spliceIt)->propertyName()); + + if (spliceTypeTagName == "") + continue; + + const auto& tmp = getSpliceTypeTagName(spliceTypeTagName, propertyName); + if (tmp != "") + return tmp; + } + + // .. if still not, check all normal children. + typedef TypeTagRegistry::ChildrenList ChildrenList; + const ChildrenList& children = TypeTagRegistry::children(typeTagName); + ChildrenList::const_iterator ttagIt = children.begin(); + for (; ttagIt != children.end(); ++ttagIt) { + const auto& tmp = getSpliceTypeTagName(*ttagIt, propertyName); + if (tmp != "") + return tmp; + } + + // if the property was not defined on a given type tag, return + // the empty string. + static std::string tmp(""); + return tmp; + } + + static const PropertyRegistryKey& getKey(const std::string& effTypeTagName, + const std::string& propertyName) + { return keys_()[effTypeTagName][propertyName]; } + + static const KeyList& getKeys(const std::string& effTypeTagName) + { return keys_()[effTypeTagName]; } + +private: + static KeyListMap& keys_() + { + static KeyListMap data; + return data; + } +}; + +#endif // !defined NO_PROPERTY_INTROSPECTION + +struct PropertyUndefined { }; +class PropertyExplicitlyUnset {}; + +template +struct Property : public PropertyUndefined +{}; + +template +struct PropertyUnset : public PropertyUndefined +{}; + +template +struct propertyExplicitlyUnset +{ + const static bool value = + std::is_base_of + >::value; +}; + +template +class propertyExplicitlyUnsetOnTree +{ + static const bool explicitlyUnset = propertyExplicitlyUnset::value; + + template + struct unsetOnAllChildren + { static const bool value = true; }; + + template + struct unsetOnAllChildren > + { static const bool value = + propertyExplicitlyUnsetOnTree::value + && unsetOnAllChildren >::value; }; + +public: + static const bool value = + (explicitlyUnset || (!Tree::isLeaf && unsetOnAllChildren::value)); +}; + +template +struct propertyExplicitlyUnsetOnTree +{ + const static bool value = std::true_type::value; +}; + +template +struct propertyDefinedOnSelf +{ + const static bool value = + ! std::is_base_of >::value; +}; + +// template class to revert the order or a std::tuple's arguments. This is required to +// make the properties of children defined on the right overwrite the properties of the +// children on the left. See https://sydius.me/2011/07/reverse-tuple-in-c/ +template +class RevertedTuple +{ +private: + template + struct RevertedTupleOuter + { + template + struct RevertedTupleInner: RevertedTupleOuter::template RevertedTupleInner { }; + }; + + template + struct RevertedTupleOuter<0, All...> + { + template + struct RevertedTupleInner { + typedef std::tuple type; + }; + }; + +public: + typedef typename RevertedTupleOuter::template RevertedTupleInner::type type; +}; + +template +class TypeTag +{ +public: + typedef SelfT SelfType; + + typedef typename RevertedTuple::type ChildrenTuple; + static const bool isLeaf = std::is_same >::value; +}; + +namespace PTag { +// this class needs to be located in the PTag namespace for reasons +// you don't really want to know... +template +struct Splices +{ + typedef typename std::tuple<> tuple; +}; +} // namespace PTag + +template +struct GetProperty +{ + // find the type tag for which the property is defined + template ::value> + struct GetEffectiveTypeTag_ + { typedef typename CurTree::SelfType type; }; + + template + struct SearchTypeTagList_; + + template + struct SearchTypeTagList_FirstThenRemaining_; + + template + struct SearchSpliceList_; + + template + struct SearchSpliceList_FirstThenRemaining_; + + // find the first type tag in a tuple for which the property is + // defined + template + struct SearchTypeTagTuple_ + { typedef void type; }; + + template + struct SearchTypeTagTuple_ > + { typedef typename SearchTypeTagList_::type type; }; + + template + struct SearchTypeTagList_ + { typedef void type; }; + + template + struct SearchTypeTagList_ + { + typedef typename SearchTypeTagList_FirstThenRemaining_< + typename GetEffectiveTypeTag_::type, + RemainingElements...>::type type; + }; + + template + struct SearchTypeTagList_FirstThenRemaining_ + { typedef EffectiveTypeTag type; }; + + template + struct SearchTypeTagList_FirstThenRemaining_ + { typedef typename SearchTypeTagList_::type type; }; + + // find the first type tag in a tuple of splices for which the + // property is defined + template + struct SearchSpliceTuple_ + { typedef void type; }; + + template + struct SearchSpliceTuple_ > + { typedef typename SearchSpliceList_::type type; }; + + template + struct SearchSpliceList_ + { typedef void type; }; + + template + struct SearchSpliceList_ + { + typedef typename SearchSpliceList_FirstThenRemaining_< + typename GetEffectiveTypeTag_::p::type>::type, + RemainingSplices...>::type type; + }; + + template + struct SearchSpliceList_FirstThenRemaining_ + { typedef EffectiveTypeTag type; }; + + template + struct SearchSpliceList_FirstThenRemaining_ + { typedef typename SearchSpliceList_::type type; }; + + // find the splice or the child type tag for which the property is defined + template ::tuple >::type > + struct SearchSplicesThenChildren_ + { typedef SpliceTypeTag type; }; + + template + struct SearchSplicesThenChildren_ + { typedef typename SearchTypeTagTuple_::type type; }; + + template + struct GetEffectiveTypeTag_ + { typedef typename SearchSplicesThenChildren_::type type; }; + +public: + typedef Property::type, PropertyTag> p; +}; + +#if !defined NO_PROPERTY_INTROSPECTION +inline int myReplaceAll_(std::string& s, + const std::string& pattern, + const std::string& replacement); +inline int myReplaceAll_(std::string& s, + const std::string& pattern, + const std::string& replacement) +{ + size_t pos; + int i = 0; + while ( (pos = s.find(pattern)) != s.npos) { + s.replace(pos, pattern.size(), replacement); + ++i; + }; + return i; +} + +inline std::string canonicalTypeTagNameToName_(const std::string& canonicalName); +inline std::string canonicalTypeTagNameToName_(const std::string& canonicalName) +{ + std::string result(canonicalName); + myReplaceAll_(result, "Opm::Properties::TTag::", "TTAG("); + myReplaceAll_(result, "::", ""); + result += ")"; + return result; +} + +inline bool getDiagnostic_(const std::string& typeTagName, + const std::string& propTagName, + std::string& result, + const std::string indent) +{ + const PropertyRegistryKey *key = 0; + + const typename PropertyRegistry::KeyList& keys = + PropertyRegistry::getKeys(typeTagName); + typename PropertyRegistry::KeyList::const_iterator it = keys.begin(); + for (; it != keys.end(); ++it) { + if (it->second.propertyName() == propTagName) { + key = &it->second; + break; + }; + } + + if (key) { + std::ostringstream oss; + oss << indent + << key->propertyKind() << " " + << key->propertyName() << " defined on '" + << canonicalTypeTagNameToName_(key->effTypeTagName()) << "' at " + << key->fileDefined() << ":" << key->lineDefined() << "\n"; + result = oss.str(); + return true; + } + + // print properties defined on children + typedef typename TypeTagRegistry::ChildrenList ChildrenList; + const ChildrenList& children = TypeTagRegistry::children(typeTagName); + ChildrenList::const_iterator ttagIt = children.begin(); + std::string newIndent = indent + " "; + for (; ttagIt != children.end(); ++ttagIt) { + if (getDiagnostic_(*ttagIt, propTagName, result, newIndent)) { + result.insert(0, indent + "Inherited from " + canonicalTypeTagNameToName_(typeTagName) + "\n"); + return true; + } + } + + return false; +} + +template +const std::string getDiagnostic(std::string propTagName) +{ + std::string result; + + std::string TypeTagName(Dune::className()); + + propTagName.replace(0, strlen("PTag("), ""); + auto n = propTagName.length(); + propTagName.replace(n - 1, 1, ""); + //TypeTagName.replace(0, strlen("Opm::Properties::TTag::"), ""); + + return result; +} + +inline void print_(const std::string& rootTypeTagName, + const std::string& curTypeTagName, + const std::string& splicePropName, + std::ostream& os, + const std::string indent, + std::set& printedProperties) +{ + if (indent == "") { + os << indent << "###########\n"; + os << indent << "# Properties\n"; + os << indent << "###########\n"; + os << indent << "Properties for " << canonicalTypeTagNameToName_(curTypeTagName) << ":"; + } + else if (splicePropName != "") + os << indent << "Inherited from splice " << splicePropName << " (set to " << canonicalTypeTagNameToName_(curTypeTagName) << "):"; + else + os << indent << "Inherited from " << canonicalTypeTagNameToName_(curTypeTagName) << ":"; + const typename PropertyRegistry::KeyList& keys = + PropertyRegistry::getKeys(curTypeTagName); + typename PropertyRegistry::KeyList::const_iterator it = keys.begin(); + bool somethingPrinted = false; + for (; it != keys.end(); ++it) { + const PropertyRegistryKey& key = it->second; + if (printedProperties.count(key.propertyName()) > 0) + continue; // property already printed + if (!somethingPrinted) { + os << "\n"; + somethingPrinted = true; + } + os << indent << " " + << key.propertyKind() << " " << key.propertyName(); + if (key.propertyKind() != "opaque") { + std::string s(key.propertyValue()); + myReplaceAll_(s, "typename ", ""); + if (myReplaceAll_(s, "::Opm::Properties::TTag::", "TTAG(")) + s += ')'; + myReplaceAll_(s, "::Opm::Properties::PTag::", ""); + myReplaceAll_(s, "::Opm::Properties::GetProperty<", "GET_PROP("); + myReplaceAll_(s, ">::p::", ")::"); + myReplaceAll_(s, "GET_PROP(TypeTag, Scalar)::type", "Scalar"); + + os << " = '" << s << "'"; + } + os << " defined at " << key.fileDefined() + << ":" << key.lineDefined() + << "\n"; + printedProperties.insert(key.propertyName()); + }; + if (!somethingPrinted) + os << " (none)\n"; + // print properties defined on splices or children + std::string newIndent = indent + " "; + + // first, iterate over the splices, ... + typedef TypeTagRegistry::SpliceList SpliceList; + const SpliceList& splices = TypeTagRegistry::splices(curTypeTagName); + SpliceList::const_iterator spliceIt = splices.begin(); + for (; spliceIt != splices.end(); ++ spliceIt) { + const auto& spliceTypeTagName = PropertyRegistry::getSpliceTypeTagName(rootTypeTagName, + (*spliceIt)->propertyName()); + print_(rootTypeTagName, spliceTypeTagName, (*spliceIt)->propertyName(), os, newIndent, printedProperties); + } + + // ... then, over the children + typedef typename TypeTagRegistry::ChildrenList ChildrenList; + const ChildrenList& children = TypeTagRegistry::children(curTypeTagName); + ChildrenList::const_iterator ttagIt = children.begin(); + for (; ttagIt != children.end(); ++ttagIt) { + print_(rootTypeTagName, *ttagIt, /*splicePropName=*/"", os, newIndent, printedProperties); + } +} + +template +void printValues(std::ostream& os = std::cout) +{ + std::set printedProps; + print_(Dune::className(), Dune::className(), /*splicePropertyName=*/"", os, /*indent=*/"", printedProps); +} +#else // !defined NO_PROPERTY_INTROSPECTION +template +void printValues(std::ostream& os = std::cout) +{ + os << + "The eWoms property system was compiled with the macro\n" + "NO_PROPERTY_INTROSPECTION defined.\n" + "No diagnostic messages this time, sorry.\n"; +} + +template +const std::string getDiagnostic(std::string propTagName) +{ + std::string result; + result = + "The eWoms property system was compiled with the macro\n" + "NO_PROPERTY_INTROSPECTION defined.\n" + "No diagnostic messages this time, sorry.\n"; + return result; +} + +#endif // !defined NO_PROPERTY_INTROSPECTION + +//! \endcond + +} // namespace Properties +} // namespace Opm + +#endif diff --git a/opm/models/utils/quadraturegeometries.hh b/opm/models/utils/quadraturegeometries.hh new file mode 100644 index 000000000..3024ebf8e --- /dev/null +++ b/opm/models/utils/quadraturegeometries.hh @@ -0,0 +1,154 @@ +// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- +// vi: set et ts=4 sw=4 sts=4: +/* + 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 2 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 . + + Consult the COPYING file in the top-level source directory of this + module for the precise wording of the license and the list of + copyright holders. +*/ +/*! + * \file + * \copydoc Opm::QuadrialteralQuadratureGeometry + */ +#ifndef EWOMS_QUADRATURE_GEOMETRIES_HH +#define EWOMS_QUADRATURE_GEOMETRIES_HH + +#include +#include +#include + +namespace Opm { +/*! + * \brief Quadrature geometry for quadrilaterals. + */ +template +class QuadrialteralQuadratureGeometry +{ +public: + enum { numCorners = (1 << dim) }; + + typedef Dune::FieldVector LocalPosition; + typedef Dune::FieldVector GlobalPosition; + + Dune::GeometryType type() const + { return Dune::GeometryType(/*topologyId=*/(1 << dim) - 1, dim); } + + template + void setCorners(const CornerContainer& corners, unsigned nCorners) + { + unsigned cornerIdx; + for (cornerIdx = 0; cornerIdx < nCorners; ++cornerIdx) { + for (unsigned j = 0; j < dim; ++j) + corners_[cornerIdx][j] = corners[cornerIdx][j]; + } + assert(cornerIdx == nCorners); + + center_ = 0; + for (cornerIdx = 0; cornerIdx < nCorners; ++cornerIdx) + center_ += corners_[cornerIdx]; + center_ /= nCorners; + } + + /*! + * \brief Returns the center of weight of the polyhedron. + */ + const GlobalPosition& center() const + { return center_; } + + /*! + * \brief Convert a local coordinate into a global one. + */ + GlobalPosition global(const LocalPosition& localPos) const + { + GlobalPosition globalPos(0.0); + + for (unsigned cornerIdx = 0; cornerIdx < numCorners; ++cornerIdx) + globalPos.axpy(cornerWeight(localPos, cornerIdx), + corners_[cornerIdx]); + + return globalPos; + } + + /*! + * \brief Returns the Jacobian matrix of the local to global + * mapping at a given local position. + */ + void jacobian(Dune::FieldMatrix& jac, + const LocalPosition& localPos) const + { + jac = 0.0; + for (unsigned cornerIdx = 0; cornerIdx < numCorners; ++cornerIdx) { + for (unsigned k = 0; k < dim; ++k) { + Scalar dWeight_dk = (cornerIdx& (1 << k)) ? 1 : -1; + for (unsigned j = 0; j < dim; ++j) { + if (k != j) { + if (cornerIdx& (1 << j)) + dWeight_dk *= localPos[j]; + else + dWeight_dk *= 1 - localPos[j]; + ; + } + } + + jac[k].axpy(dWeight_dk, corners_[cornerIdx]); + } + } + } + + /*! + * \brief Return the determinant of the Jacobian of the mapping + * from local to global coordinates at a given local + * position. + */ + Scalar integrationElement(const LocalPosition& localPos) const + { + Dune::FieldMatrix jac; + jacobian(jac, localPos); + return jac.determinant(); + } + + /*! + * \brief Return the position of the corner with a given index + */ + const GlobalPosition& corner(unsigned cornerIdx) const + { return corners_[cornerIdx]; } + + /*! + * \brief Return the weight of an individual corner for the local + * to global mapping. + */ + Scalar cornerWeight(const LocalPosition& localPos, unsigned cornerIdx) const + { + GlobalPosition globalPos(0.0); + + // this code is based on the Q1 finite element code from + // dune-localfunctions + Scalar weight = 1.0; + for (unsigned j = 0; j < dim; ++j) + weight *= (cornerIdx& (1 << j)) ? localPos[j] : (1 - localPos[j]); + + return weight; + } + +private: + GlobalPosition corners_[numCorners]; + GlobalPosition center_; +}; + +} // namespace Opm + +#endif // EWOMS_QUADRATURE_GEOMETRY_HH diff --git a/opm/models/utils/signum.hh b/opm/models/utils/signum.hh new file mode 100644 index 000000000..6d51cfe5d --- /dev/null +++ b/opm/models/utils/signum.hh @@ -0,0 +1,45 @@ +// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- +// vi: set et ts=4 sw=4 sts=4: +/* + 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 2 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 . + + Consult the COPYING file in the top-level source directory of this + module for the precise wording of the license and the list of + copyright holders. +*/ +/*! + * \file + * \copydoc sgn() + */ +#ifndef EWOMS_SIGNUM_HH +#define EWOMS_SIGNUM_HH + +namespace Opm { +/*! + * \brief Template function which returns the sign of a floating point value + * + * This is a type safe and fast implementation of a sign() function for arbitrary + * floating point values. It a slightly modified variant of + * + * https://stackoverflow.com/questions/1903954/is-there-a-standard-sign-function-signum-sgn-in-c-c + */ +template +int signum(Scalar val) +{ return (0 < val) - (val < 0); } + +} // namespace Opm + +#endif diff --git a/opm/models/utils/simulator.hh b/opm/models/utils/simulator.hh new file mode 100644 index 000000000..a2ffe073c --- /dev/null +++ b/opm/models/utils/simulator.hh @@ -0,0 +1,957 @@ +// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- +// vi: set et ts=4 sw=4 sts=4: +/* + 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 2 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 . + + Consult the COPYING file in the top-level source directory of this + module for the precise wording of the license and the list of + copyright holders. +*/ +/*! + * \file + * + * \copydoc Opm::Simulator + */ +#ifndef EWOMS_SIMULATOR_HH +#define EWOMS_SIMULATOR_HH + +#include +#include + +#include +#include +#include + +#include +#include + +#include +#include +#include +#include +#include +#include + +BEGIN_PROPERTIES + +NEW_PROP_TAG(Scalar); +NEW_PROP_TAG(Vanguard); +NEW_PROP_TAG(GridView); +NEW_PROP_TAG(Model); +NEW_PROP_TAG(Problem); +NEW_PROP_TAG(EndTime); +NEW_PROP_TAG(RestartTime); +NEW_PROP_TAG(InitialTimeStepSize); +NEW_PROP_TAG(PredeterminedTimeStepsFile); + +END_PROPERTIES + +#define EWOMS_CATCH_PARALLEL_EXCEPTIONS_FATAL(code) \ + { \ + const auto& comm = Dune::MPIHelper::getCollectiveCommunication(); \ + bool exceptionThrown = false; \ + try { code; } \ + catch (const Dune::Exception& e) { \ + exceptionThrown = true; \ + std::cerr << "Process " << comm.rank() << " threw a fatal exception: " \ + << e.what() << ". Abort!" << std::endl; \ + } \ + catch (const std::exception& e) { \ + exceptionThrown = true; \ + std::cerr << "Process " << comm.rank() << " threw a fatal exception: " \ + << e.what() << ". Abort!" << std::endl; \ + } \ + catch (...) { \ + exceptionThrown = true; \ + std::cerr << "Process " << comm.rank() << " threw a fatal exception. " \ + <<" Abort!" << std::endl; \ + } \ + \ + if (comm.max(exceptionThrown)) \ + std::abort(); \ + } + +namespace Opm { + +/*! + * \ingroup Common + * + * \brief Manages the initializing and running of time dependent + * problems. + * + * This class instantiates the grid, the model and the problem to be + * simlated and runs the simulation loop. The time axis is treated as + * a sequence of "episodes" which are defined as time intervals for + * which the problem exhibits boundary conditions and source terms + * that do not depend on time. + */ +template +class Simulator +{ + typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar; + typedef typename GET_PROP_TYPE(TypeTag, Vanguard) Vanguard; + typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView; + typedef typename GET_PROP_TYPE(TypeTag, Model) Model; + typedef typename GET_PROP_TYPE(TypeTag, Problem) Problem; + +public: + // do not allow to copy simulators around + Simulator(const Simulator& ) = delete; + + Simulator(bool verbose = true) + { + Opm::TimerGuard setupTimerGuard(setupTimer_); + + setupTimer_.start(); + + const auto& comm = Dune::MPIHelper::getCollectiveCommunication(); + verbose_ = verbose && comm.rank() == 0; + + timeStepIdx_ = 0; + startTime_ = 0.0; + time_ = 0.0; + endTime_ = EWOMS_GET_PARAM(TypeTag, Scalar, EndTime); + timeStepSize_ = EWOMS_GET_PARAM(TypeTag, Scalar, InitialTimeStepSize); + + const std::string& predetTimeStepFile = + EWOMS_GET_PARAM(TypeTag, std::string, PredeterminedTimeStepsFile); + if (!predetTimeStepFile.empty()) { + std::ifstream is(predetTimeStepFile); + while (!is.eof()) { + Scalar dt; + is >> dt; + forcedTimeSteps_.push_back(dt); + } + } + + episodeIdx_ = 0; + episodeStartTime_ = 0; + episodeLength_ = std::numeric_limits::max(); + + finished_ = false; + + if (verbose_) + std::cout << "Allocating the simulation vanguard\n" << std::flush; + + int exceptionThrown = 0; + try + { vanguard_.reset(new Vanguard(*this)); } + catch (const std::exception& e) { + exceptionThrown = 1; + if (verbose_) + std::cerr << "Rank " << comm.rank() << " threw an exception: " << e.what() << std::endl; + } + + if (comm.max(exceptionThrown)) + throw std::runtime_error("Allocating the simulation vanguard failed."); + + if (verbose_) + std::cout << "Distributing the vanguard's data\n" << std::flush; + + try + { vanguard_->loadBalance(); } + catch (const std::exception& e) { + exceptionThrown = 1; + if (verbose_) + std::cerr << "Rank " << comm.rank() << " threw an exception: " << e.what() << std::endl; + } + + if (comm.max(exceptionThrown)) + throw std::runtime_error("Could not distribute the vanguard data."); + + if (verbose_) + std::cout << "Allocating the model\n" << std::flush; + model_.reset(new Model(*this)); + + if (verbose_) + std::cout << "Allocating the problem\n" << std::flush; + problem_.reset(new Problem(*this)); + + if (verbose_) + std::cout << "Initializing the model\n" << std::flush; + + try + { model_->finishInit(); } + catch (const std::exception& e) { + exceptionThrown = 1; + if (verbose_) + std::cerr << "Rank " << comm.rank() << " threw an exception: " << e.what() << std::endl; + } + + if (comm.max(exceptionThrown)) + throw std::runtime_error("Could not initialize the model."); + + if (verbose_) + std::cout << "Initializing the problem\n" << std::flush; + + try + { problem_->finishInit(); } + catch (const std::exception& e) { + exceptionThrown = 1; + if (verbose_) + std::cerr << "Rank " << comm.rank() << " threw an exception: " << e.what() << std::endl; + } + + if (comm.max(exceptionThrown)) + throw std::runtime_error("Could not initialize the problem."); + + setupTimer_.stop(); + + if (verbose_) + std::cout << "Simulator successfully set up\n" << std::flush; + } + + /*! + * \brief Registers all runtime parameters used by the simulation. + */ + static void registerParameters() + { + EWOMS_REGISTER_PARAM(TypeTag, Scalar, EndTime, + "The simulation time at which the simulation is finished [s]"); + EWOMS_REGISTER_PARAM(TypeTag, Scalar, InitialTimeStepSize, + "The size of the initial time step [s]"); + EWOMS_REGISTER_PARAM(TypeTag, Scalar, RestartTime, + "The simulation time at which a restart should be attempted [s]"); + EWOMS_REGISTER_PARAM(TypeTag, std::string, PredeterminedTimeStepsFile, + "A file with a list of predetermined time step sizes (one " + "time step per line)"); + + Vanguard::registerParameters(); + Model::registerParameters(); + Problem::registerParameters(); + } + + /*! + * \brief Return a reference to the grid manager of simulation + */ + Vanguard& vanguard() + { return *vanguard_; } + + /*! + * \brief Return a reference to the grid manager of simulation + */ + const Vanguard& vanguard() const + { return *vanguard_; } + + /*! + * \brief Return the grid view for which the simulation is done + */ + const GridView& gridView() const + { return vanguard_->gridView(); } + + /*! + * \brief Return the physical model used in the simulation + */ + Model& model() + { return *model_; } + + /*! + * \brief Return the physical model used in the simulation + */ + const Model& model() const + { return *model_; } + + /*! + * \brief Return the object which specifies the pysical setup of + * the simulation + */ + Problem& problem() + { return *problem_; } + + /*! + * \brief Return the object which specifies the pysical setup of + * the simulation + */ + const Problem& problem() const + { return *problem_; } + + /*! + * \brief Set the time of the start of the simulation. + * + * \param t The time \f$\mathrm{[s]}\f$ which should be jumped to + */ + void setStartTime(Scalar t) + { startTime_ = t; } + + /*! + * \brief Return the time of the start of the simulation. + */ + Scalar startTime() const + { return startTime_; } + + /*! + * \brief Set the current simulated time, don't change the current + * time step index. + * + * \param t The time \f$\mathrm{[s]}\f$ which should be jumped to + */ + void setTime(Scalar t) + { time_ = t; } + + /*! + * \brief Set the current simulated time and the time step index. + * + * \param t The time \f$\mathrm{[s]}\f$ which should be jumped to + * \param stepIdx The new time step index + */ + void setTime(Scalar t, unsigned stepIdx) + { + time_ = t; + timeStepIdx_ = stepIdx; + } + + /*! + * \brief Return the number of seconds of simulated time which have elapsed since the + * start time. + * + * To get the time after the time integration, you have to add + * timeStepSize() to time(). + */ + Scalar time() const + { return time_; } + + /*! + * \brief Set the time of simulated seconds at which the simulation runs. + * + * \param t The time \f$\mathrm{[s]}\f$ at which the simulation is finished + */ + void setEndTime(Scalar t) + { endTime_ = t; } + + /*! + * \brief Returns the number of (simulated) seconds which the simulation + * runs. + */ + Scalar endTime() const + { return endTime_; } + + /*! + * \brief Returns a reference to the timer object which measures the time needed to + * set up and initialize the simulation + */ + const Opm::Timer& setupTimer() const + { return setupTimer_; } + + /*! + * \brief Returns a reference to the timer object which measures the time needed to + * run the simulation + */ + const Opm::Timer& executionTimer() const + { return executionTimer_; } + Opm::Timer& executionTimer() + { return executionTimer_; } + + /*! + * \brief Returns a reference to the timer object which measures the time needed for + * pre- and postprocessing of the solutions. + */ + const Opm::Timer& prePostProcessTimer() const + { return prePostProcessTimer_; } + + /*! + * \brief Returns a reference to the timer object which measures the time needed for + * linarizing the solutions. + */ + const Opm::Timer& linearizeTimer() const + { return linearizeTimer_; } + + /*! + * \brief Returns a reference to the timer object which measures the time needed by + * the solver. + */ + const Opm::Timer& solveTimer() const + { return solveTimer_; } + + /*! + * \brief Returns a reference to the timer object which measures the time needed to + * the solutions of the non-linear system of equations. + */ + const Opm::Timer& updateTimer() const + { return updateTimer_; } + + /*! + * \brief Returns a reference to the timer object which measures the time needed to + * write the visualization output + */ + const Opm::Timer& writeTimer() const + { return writeTimer_; } + + /*! + * \brief Set the current time step size to a given value. + * + * If the step size would exceed the length of the current + * episode, the timeStep() method will take care that the step + * size won't exceed the episode or the end of the simulation, + * though. + * + * \param timeStepSize The new value for the time step size \f$\mathrm{[s]}\f$ + */ + void setTimeStepSize(Scalar value) + { timeStepSize_ = value; } + + /*! + * \brief Set the current time step index to a given value. + * + * \param timeStepIndex The new value for the time step index + */ + void setTimeStepIndex(unsigned value) + { timeStepIdx_ = value; } + + /*! + * \brief Returns the time step length \f$\mathrm{[s]}\f$ so that we + * don't miss the beginning of the next episode or cross + * the end of the simlation. + */ + Scalar timeStepSize() const + { return timeStepSize_; } + + /*! + * \brief Returns number of time steps which have been + * executed since the beginning of the simulation. + */ + int timeStepIndex() const + { return timeStepIdx_; } + + /*! + * \brief Specify whether the simulation is finished + * + * \param yesno If true the simulation is considered finished + * before the end time is reached, else it is only + * considered finished if the end time is reached. + */ + void setFinished(bool yesno = true) + { finished_ = yesno; } + + /*! + * \brief Returns true if the simulation is finished. + * + * This is the case if either setFinished(true) has been called or + * if the end time is reached. + */ + bool finished() const + { + assert(timeStepSize_ >= 0.0); + Scalar eps = + std::max(Scalar(std::abs(this->time())), timeStepSize()) + *std::numeric_limits::epsilon()*1e3; + return finished_ || (this->time()*(1.0 + eps) >= endTime()); + } + + /*! + * \brief Returns true if the simulation is finished after the + * time level is incremented by the current time step size. + */ + bool willBeFinished() const + { + static const Scalar eps = std::numeric_limits::epsilon()*1e3; + + return finished_ || (this->time() + timeStepSize_)*(1.0 + eps) >= endTime(); + } + + /*! + * \brief Aligns the time step size to the episode boundary and to + * the end time of the simulation. + */ + Scalar maxTimeStepSize() const + { + if (finished()) + return 0.0; + + return std::min(episodeMaxTimeStepSize(), + std::max(0.0, endTime() - this->time())); + } + + /*! + * \brief Change the current episode of the simulation. + * + * \param episodeStartTime Time when the episode began \f$\mathrm{[s]}\f$ + * \param episodeLength Length of the episode \f$\mathrm{[s]}\f$ + */ + void startNextEpisode(Scalar episodeStartTime, Scalar episodeLength) + { + ++episodeIdx_; + episodeStartTime_ = episodeStartTime; + episodeLength_ = episodeLength; + } + + /*! + * \brief Start the next episode, but don't change the episode + * identifier. + * + * \param len Length of the episode \f$\mathrm{[s]}\f$, infinite if not + * specified. + */ + void startNextEpisode(Scalar len = std::numeric_limits::max()) + { + ++episodeIdx_; + episodeStartTime_ = startTime_ + time_; + episodeLength_ = len; + } + + /*! + * \brief Sets the index of the current episode. + * + * Use this method with care! + */ + void setEpisodeIndex(int episodeIdx) + { episodeIdx_ = episodeIdx; } + + /*! + * \brief Returns the index of the current episode. + * + * The first episode has the index 0. + */ + int episodeIndex() const + { return episodeIdx_; } + + /*! + * \brief Returns the absolute time when the current episode + * started \f$\mathrm{[s]}\f$. + */ + Scalar episodeStartTime() const + { return episodeStartTime_; } + + /*! + * \brief Sets the length in seconds of the current episode. + * + * Use this method with care! + */ + void setEpisodeLength(Scalar dt) + { episodeLength_ = dt; } + + /*! + * \brief Returns the length of the current episode in + * simulated time \f$\mathrm{[s]}\f$. + */ + Scalar episodeLength() const + { return episodeLength_; } + + /*! + * \brief Returns true if the current episode has just been started at the + * current time. + */ + bool episodeStarts() const + { + static const Scalar eps = std::numeric_limits::epsilon()*1e3; + + return this->time() <= (episodeStartTime_ - startTime())*(1 + eps); + } + + /*! + * \brief Returns true if the current episode is finished at the + * current time. + */ + bool episodeIsOver() const + { + static const Scalar eps = std::numeric_limits::epsilon()*1e3; + + return this->time() >= (episodeStartTime_ - startTime() + episodeLength())*(1 - eps); + } + + /*! + * \brief Returns true if the current episode will be finished + * after the current time step. + */ + bool episodeWillBeOver() const + { + static const Scalar eps = std::numeric_limits::epsilon()*1e3; + + return this->time() + timeStepSize() + >= (episodeStartTime_ - startTime() + episodeLength())*(1 - eps); + } + + /*! + * \brief Aligns the time step size to the episode boundary if the + * current time step crosses the boundary of the current episode. + */ + Scalar episodeMaxTimeStepSize() const + { + // if the current episode is over and the simulation + // wants to give it some extra time, we will return + // the time step size it suggested instead of trying + // to align it to the end of the episode. + if (episodeIsOver()) + return 0.0; + + // make sure that we don't exceed the end of the + // current episode. + return std::max(0.0, + (episodeStartTime() + episodeLength()) + - (this->time() + this->startTime())); + } + + /* + * \} + */ + + /*! + * \brief Runs the simulation using a given problem class. + * + * This method makes sure that time steps sizes are aligned to + * episode boundaries, amongst other stuff. + */ + void run() + { + // create TimerGuard objects to hedge for exceptions + TimerGuard setupTimerGuard(setupTimer_); + TimerGuard executionTimerGuard(executionTimer_); + TimerGuard prePostProcessTimerGuard(prePostProcessTimer_); + TimerGuard writeTimerGuard(writeTimer_); + + setupTimer_.start(); + Scalar restartTime = EWOMS_GET_PARAM(TypeTag, Scalar, RestartTime); + if (restartTime > -1e30) { + // try to restart a previous simulation + time_ = restartTime; + + Opm::Restart res; + EWOMS_CATCH_PARALLEL_EXCEPTIONS_FATAL(res.deserializeBegin(*this, time_)); + if (verbose_) + std::cout << "Deserialize from file '" << res.fileName() << "'\n" << std::flush; + EWOMS_CATCH_PARALLEL_EXCEPTIONS_FATAL(this->deserialize(res)); + EWOMS_CATCH_PARALLEL_EXCEPTIONS_FATAL(problem_->deserialize(res)); + EWOMS_CATCH_PARALLEL_EXCEPTIONS_FATAL(model_->deserialize(res)); + EWOMS_CATCH_PARALLEL_EXCEPTIONS_FATAL(res.deserializeEnd()); + if (verbose_) + std::cout << "Deserialization done." + << " Simulator time: " << time() << humanReadableTime(time()) + << " Time step index: " << timeStepIndex() + << " Episode index: " << episodeIndex() + << "\n" << std::flush; + } + else { + // if no restart is done, apply the initial solution + if (verbose_) + std::cout << "Applying the initial solution of the \"" << problem_->name() + << "\" problem\n" << std::flush; + + Scalar oldTimeStepSize = timeStepSize_; + int oldTimeStepIdx = timeStepIdx_; + timeStepSize_ = 0.0; + timeStepIdx_ = -1; + + EWOMS_CATCH_PARALLEL_EXCEPTIONS_FATAL(model_->applyInitialSolution()); + + // write initial condition + if (problem_->shouldWriteOutput()) + EWOMS_CATCH_PARALLEL_EXCEPTIONS_FATAL(problem_->writeOutput()); + + timeStepSize_ = oldTimeStepSize; + timeStepIdx_ = oldTimeStepIdx; + } + setupTimer_.stop(); + + executionTimer_.start(); + bool episodeBegins = episodeIsOver() || (timeStepIdx_ == 0); + // do the time steps + while (!finished()) { + prePostProcessTimer_.start(); + if (episodeBegins) { + // notify the problem that a new episode has just been + // started. + EWOMS_CATCH_PARALLEL_EXCEPTIONS_FATAL(problem_->beginEpisode()); + + if (finished()) { + // the problem can chose to terminate the simulation in + // beginEpisode(), so we have handle this case. + EWOMS_CATCH_PARALLEL_EXCEPTIONS_FATAL(problem_->endEpisode()); + prePostProcessTimer_.stop(); + + break; + } + } + episodeBegins = false; + + if (verbose_) { + std::cout << "Begin time step " << timeStepIndex() + 1 << ". " + << "Start time: " << this->time() << " seconds" << humanReadableTime(this->time()) + << ", step size: " << timeStepSize() << " seconds" << humanReadableTime(timeStepSize()) + << "\n"; + } + + // pre-process the current solution + EWOMS_CATCH_PARALLEL_EXCEPTIONS_FATAL(problem_->beginTimeStep()); + + if (finished()) { + // the problem can chose to terminate the simulation in + // beginTimeStep(), so we have handle this case. + EWOMS_CATCH_PARALLEL_EXCEPTIONS_FATAL(problem_->endTimeStep()); + EWOMS_CATCH_PARALLEL_EXCEPTIONS_FATAL(problem_->endEpisode()); + prePostProcessTimer_.stop(); + + break; + } + prePostProcessTimer_.stop(); + + try { + // execute the time integration scheme + problem_->timeIntegration(); + } + catch (...) { + // exceptions in the time integration might be recoverable. clean up in + // case they are + const auto& model = problem_->model(); + prePostProcessTimer_ += model.prePostProcessTimer(); + linearizeTimer_ += model.linearizeTimer(); + solveTimer_ += model.solveTimer(); + updateTimer_ += model.updateTimer(); + + throw; + } + + const auto& model = problem_->model(); + prePostProcessTimer_ += model.prePostProcessTimer(); + linearizeTimer_ += model.linearizeTimer(); + solveTimer_ += model.solveTimer(); + updateTimer_ += model.updateTimer(); + + // post-process the current solution + prePostProcessTimer_.start(); + EWOMS_CATCH_PARALLEL_EXCEPTIONS_FATAL(problem_->endTimeStep()); + prePostProcessTimer_.stop(); + + // write the result to disk + writeTimer_.start(); + if (problem_->shouldWriteOutput()) + EWOMS_CATCH_PARALLEL_EXCEPTIONS_FATAL(problem_->writeOutput()); + writeTimer_.stop(); + + // do the next time integration + Scalar oldDt = timeStepSize(); + EWOMS_CATCH_PARALLEL_EXCEPTIONS_FATAL(problem_->advanceTimeLevel()); + + if (verbose_) { + std::cout << "Time step " << timeStepIndex() + 1 << " done. " + << "CPU time: " << executionTimer_.realTimeElapsed() << " seconds" << humanReadableTime(executionTimer_.realTimeElapsed()) + << ", end time: " << this->time() + oldDt << " seconds" << humanReadableTime(this->time() + oldDt) + << ", step size: " << oldDt << " seconds" << humanReadableTime(oldDt) + << "\n" << std::flush; + } + + // advance the simulated time by the current time step size + time_ += oldDt; + ++timeStepIdx_; + + prePostProcessTimer_.start(); + // notify the problem if an episode is finished + if (episodeIsOver()) { + // Notify the problem about the end of the current episode... + EWOMS_CATCH_PARALLEL_EXCEPTIONS_FATAL(problem_->endEpisode()); + episodeBegins = true; + } + else { + Scalar dt; + if (timeStepIdx_ < static_cast(forcedTimeSteps_.size())) + // use the next time step size from the input file + dt = forcedTimeSteps_[timeStepIdx_]; + else + // ask the problem to provide the next time step size + dt = std::min(maxTimeStepSize(), problem_->nextTimeStepSize()); + + setTimeStepSize(dt); + } + prePostProcessTimer_.stop(); + + // write restart file if mandated by the problem + writeTimer_.start(); + if (problem_->shouldWriteRestartFile()) + EWOMS_CATCH_PARALLEL_EXCEPTIONS_FATAL(serialize()); + writeTimer_.stop(); + } + executionTimer_.stop(); + + EWOMS_CATCH_PARALLEL_EXCEPTIONS_FATAL(problem_->finalize()); + } + + /*! + * \brief Given a time step size in seconds, return it in a format which is more + * easily parsable by humans. + * + * e.g. 874000.0 will become "10.12 days" + */ + static std::string humanReadableTime(Scalar timeInSeconds, bool isAmendment=true) + { + std::ostringstream oss; + oss << std::setprecision(4); + if (isAmendment) + oss << " ("; + if (timeInSeconds >= 365.25*24*60*60) { + int years = static_cast(timeInSeconds/(365.25*24*60*60)); + int days = static_cast((timeInSeconds - years*(365.25*24*60*60))/(24*60*60)); + + double accuracy = 1e-2; + double hours = + std::round(1.0/accuracy* + (timeInSeconds + - years*(365.25*24*60*60) + - days*(24*60*60))/(60*60)) + *accuracy; + + oss << years << " years, " << days << " days, " << hours << " hours"; + } + else if (timeInSeconds >= 24.0*60*60) { + int days = static_cast(timeInSeconds/(24*60*60)); + int hours = static_cast((timeInSeconds - days*(24*60*60))/(60*60)); + + double accuracy = 1e-2; + double minutes = + std::round(1.0/accuracy* + (timeInSeconds + - days*(24*60*60) + - hours*(60*60))/60) + *accuracy; + + oss << days << " days, " << hours << " hours, " << minutes << " minutes"; + } + else if (timeInSeconds >= 60.0*60) { + int hours = static_cast(timeInSeconds/(60*60)); + int minutes = static_cast((timeInSeconds - hours*(60*60))/60); + + double accuracy = 1e-2; + double seconds = + std::round(1.0/accuracy* + (timeInSeconds + - hours*(60*60) + - minutes*60)) + *accuracy; + + oss << hours << " hours, " << minutes << " minutes, " << seconds << " seconds"; + } + else if (timeInSeconds >= 60.0) { + int minutes = static_cast(timeInSeconds/60); + + double accuracy = 1e-3; + double seconds = + std::round(1.0/accuracy* + (timeInSeconds + - minutes*60)) + *accuracy; + + oss << minutes << " minutes, " << seconds << " seconds"; + } + else if (!isAmendment) + oss << timeInSeconds << " seconds"; + else + return ""; + if (isAmendment) + oss << ")"; + + return oss.str(); + } + + /*! + * \name Saving/restoring the simulation state + * \{ + */ + + /*! + * \brief This method writes the complete state of the simulation + * to the harddisk. + * + * The file will start with the prefix returned by the name() + * method, has the current time of the simulation clock in it's + * name and uses the extension .ers. (Ewoms ReStart + * file.) See Opm::Restart for details. + */ + void serialize() + { + typedef Opm::Restart Restarter; + Restarter res; + res.serializeBegin(*this); + if (gridView().comm().rank() == 0) + std::cout << "Serialize to file '" << res.fileName() << "'" + << ", next time step size: " << timeStepSize() + << "\n" << std::flush; + + this->serialize(res); + problem_->serialize(res); + model_->serialize(res); + res.serializeEnd(); + } + + /*! + * \brief Write the time manager's state to a restart file. + * + * \tparam Restarter The type of the object which takes care to serialize + * data + * \param restarter The serializer object + */ + template + void serialize(Restarter& restarter) + { + restarter.serializeSectionBegin("Simulator"); + restarter.serializeStream() + << episodeIdx_ << " " + << episodeStartTime_ << " " + << episodeLength_ << " " + << startTime_ << " " + << time_ << " " + << timeStepIdx_ << " "; + restarter.serializeSectionEnd(); + } + + /*! + * \brief Read the time manager's state from a restart file. + * + * \tparam Restarter The type of the object which takes care to deserialize + * data + * \param restarter The deserializer object + */ + template + void deserialize(Restarter& restarter) + { + restarter.deserializeSectionBegin("Simulator"); + restarter.deserializeStream() + >> episodeIdx_ + >> episodeStartTime_ + >> episodeLength_ + >> startTime_ + >> time_ + >> timeStepIdx_; + restarter.deserializeSectionEnd(); + } + +private: + std::unique_ptr vanguard_; + std::unique_ptr model_; + std::unique_ptr problem_; + + int episodeIdx_; + Scalar episodeStartTime_; + Scalar episodeLength_; + + Opm::Timer setupTimer_; + Opm::Timer executionTimer_; + Opm::Timer prePostProcessTimer_; + Opm::Timer linearizeTimer_; + Opm::Timer solveTimer_; + Opm::Timer updateTimer_; + Opm::Timer writeTimer_; + + std::vector forcedTimeSteps_; + Scalar startTime_; + Scalar time_; + Scalar endTime_; + + Scalar timeStepSize_; + int timeStepIdx_; + + bool finished_; + bool verbose_; +}; +} // namespace Opm + +#endif diff --git a/opm/models/utils/start.hh b/opm/models/utils/start.hh new file mode 100644 index 000000000..11fab3f71 --- /dev/null +++ b/opm/models/utils/start.hh @@ -0,0 +1,443 @@ +// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- +// vi: set et ts=4 sw=4 sts=4: +/* + 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 2 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 . + + Consult the COPYING file in the top-level source directory of this + module for the precise wording of the license and the list of + copyright holders. +*/ +/*! + * \file + * \brief Provides convenience routines to bring up the simulation at runtime. + */ +#ifndef EWOMS_START_HH +#define EWOMS_START_HH + +#include +// the following header is not required here, but it must be included before +// dune/common/densematrix.hh because of some c++ ideosyncrasies +#include + +#include "parametersystem.hh" + +#include +#include +#include +#include + +#include + +#include + +#include +#include +#include +#include + +#if HAVE_DUNE_FEM +#include +#endif + +#include +#include +#include +#include +#include + +#include +#include +#include +#include +#include + +#if HAVE_MPI +#include +#endif + +BEGIN_PROPERTIES + +// forward declaration of property tags +NEW_PROP_TAG(Scalar); +NEW_PROP_TAG(Simulator); +NEW_PROP_TAG(ThreadManager); +NEW_PROP_TAG(PrintProperties); +NEW_PROP_TAG(PrintParameters); +NEW_PROP_TAG(ParameterFile); +NEW_PROP_TAG(Problem); + +END_PROPERTIES + +//! \cond SKIP_THIS + +namespace Opm { +/*! + * \brief Announce all runtime parameters to the registry but do not specify them yet. + */ +template +static inline void registerAllParameters_(bool finalizeRegistration = true) +{ + typedef typename GET_PROP_TYPE(TypeTag, Simulator) Simulator; + typedef typename GET_PROP_TYPE(TypeTag, ThreadManager) ThreadManager; + + EWOMS_REGISTER_PARAM(TypeTag, std::string, ParameterFile, + "An .ini file which contains a set of run-time " + "parameters"); + EWOMS_REGISTER_PARAM(TypeTag, int, PrintProperties, + "Print the values of the compile time properties at " + "the start of the simulation"); + EWOMS_REGISTER_PARAM(TypeTag, int, PrintParameters, + "Print the values of the run-time parameters at the " + "start of the simulation"); + + Simulator::registerParameters(); + ThreadManager::registerParameters(); + + if (finalizeRegistration) + EWOMS_END_PARAM_REGISTRATION(TypeTag); +} + +/*! + * \brief Register all runtime parameters, parse the command line + * arguments and the parameter file. + * + * \param argc The number of command line arguments + * \param argv Array with the command line argument strings + */ +template +static inline int setupParameters_(int argc, + const char **argv, + bool registerParams=true, + bool allowUnused=false, + bool handleHelp = true) +{ + typedef typename GET_PROP_TYPE(TypeTag, Problem) Problem; + + // first, get the MPI rank of the current process + int myRank = 0; + + //////////////////////////////////////////////////////////// + // Register all parameters + //////////////////////////////////////////////////////////// + if (registerParams) + registerAllParameters_(); + + //////////////////////////////////////////////////////////// + // set the parameter values + //////////////////////////////////////////////////////////// + + // fill the parameter tree with the options from the command line + const auto& positionalParamCallback = Problem::handlePositionalParameter; + std::string helpPreamble = ""; + if (myRank == 0 && handleHelp) + helpPreamble = Problem::helpPreamble(argc, argv); + std::string s = + Parameters::parseCommandLineOptions(argc, + argv, + helpPreamble, + positionalParamCallback); + if (!s.empty()) + return /*status=*/1; + + std::string paramFileName = EWOMS_GET_PARAM_(TypeTag, std::string, ParameterFile); + if (paramFileName != "") { + //////////////////////////////////////////////////////////// + // add the parameters specified using an .ini file + //////////////////////////////////////////////////////////// + + // check whether the parameter file is readable. + std::ifstream tmp; + tmp.open(paramFileName.c_str()); + if (!tmp.is_open()) { + std::ostringstream oss; + if (myRank == 0) { + oss << "Parameter file \"" << paramFileName + << "\" does not exist or is not readable."; + Parameters::printUsage(argv[0], oss.str()); + } + return /*status=*/1; + } + + // read the parameter file. + Parameters::parseParameterFile(paramFileName, /*overwrite=*/false); + } + + // make sure that no unknown parameters are encountered + typedef std::pair KeyValuePair; + typedef std::list ParamList; + + ParamList usedParams; + ParamList unusedParams; + + EWOMS_GET_PARAM_LISTS(TypeTag, usedParams, unusedParams); + if (!allowUnused && !unusedParams.empty()) { + if (myRank == 0) { + if (unusedParams.size() == 1) + std::cerr << "The following explicitly specified parameter is unknown:\n"; + else + std::cerr << "The following " << unusedParams.size() + << " explicitly specified parameters are unknown:\n"; + + std::cerr << "\n"; + for (const auto& keyValue : unusedParams) + std::cerr << " " << keyValue.first << "=\"" << keyValue.second << "\"\n"; + std::cerr << "\n"; + + std::cerr << "Use\n" + << "\n" + << " " << argv[0] << " --help\n" + << "\n" + <<"to obtain the list of recognized command line parameters.\n\n"; + } + return /*status=*/1; + } + + return /*status=*/0; +} + +/*! + * \brief Resets the current TTY to a usable state if the program was aborted. + * + * This is intended to be called as part of a generic exception handler + */ +static inline void resetTerminal_() +{ + // make sure stderr and stderr do not contain any unwritten data and make sure that + // the TTY does not see any unfinished ANSI escape sequence. + std::cerr << " \r\n"; + std::cerr.flush(); + std::cout << " \r\n"; + std::cout.flush(); + + // it seems like some terminals sometimes takes their time to react, so let's + // accommodate them. + usleep(/*usec=*/500*1000); + + // this requires the 'stty' command to be available in the command search path. on + // most linux systems, is the case. (but even if the system() function fails, the + // worst thing which can happen is that the TTY stays potentially choked up...) + if (system("stty sane") != 0) + std::cout << "Executing the 'stty' command failed." + << " Terminal might be left in an undefined state!\n"; +} + +/*! + * \brief Resets the current TTY to a usable state if the program was interrupted by + * SIGABRT or SIGINT. + */ +static inline void resetTerminal_(int signum) +{ + // first thing to do when a nuke hits: restore the default signal handler + signal(signum, SIG_DFL); + +#if HAVE_MPI + int rank = 0; + MPI_Comm_rank(MPI_COMM_WORLD, &rank); + if (rank != 0) { + // re-raise the signal + raise(signum); + + return; + } +#endif + + if (isatty(fileno(stdout)) && isatty(fileno(stdin))) { + std::cout << "\n\nReceived signal " << signum + << " (\"" << strsignal(signum) << "\")." + << " Trying to reset the terminal.\n"; + + resetTerminal_(); + } + + // after we did our best to clean the pedestrian way, re-raise the signal + raise(signum); +} +//! \endcond + +/*! + * \ingroup Common + * + * \brief Provides a main function which reads in parameters from the + * command line and a parameter file and runs the simulation + * + * \tparam TypeTag The type tag of the problem which needs to be solved + * + * \param argc The number of command line arguments + * \param argv The array of the command line arguments + */ +template +static inline int start(int argc, char **argv) +{ + typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar; + typedef typename GET_PROP_TYPE(TypeTag, Simulator) Simulator; + typedef typename GET_PROP_TYPE(TypeTag, Problem) Problem; + typedef typename GET_PROP_TYPE(TypeTag, ThreadManager) ThreadManager; + + // set the signal handlers to reset the TTY to a well defined state on unexpected + // program aborts + if (isatty(STDIN_FILENO)) { + signal(SIGINT, resetTerminal_); + signal(SIGHUP, resetTerminal_); + signal(SIGABRT, resetTerminal_); + signal(SIGFPE, resetTerminal_); + signal(SIGSEGV, resetTerminal_); + signal(SIGPIPE, resetTerminal_); + signal(SIGTERM, resetTerminal_); + } + + Opm::resetLocale(); + + int myRank = 0; + try + { + int paramStatus = setupParameters_(argc, const_cast(argv)); + if (paramStatus == 1) + return 1; + if (paramStatus == 2) + return 0; + + ThreadManager::init(); + + // initialize MPI, finalize is done automatically on exit +#if HAVE_DUNE_FEM + Dune::Fem::MPIManager::initialize(argc, argv); + myRank = Dune::Fem::MPIManager::rank(); +#else + myRank = Dune::MPIHelper::instance(argc, argv).rank(); +#endif + + // read the initial time step and the end time + Scalar endTime = EWOMS_GET_PARAM(TypeTag, Scalar, EndTime); + if (endTime < -1e50) { + if (myRank == 0) + Parameters::printUsage(argv[0], + "Mandatory parameter '--end-time' not specified!"); + return 1; + } + + Scalar initialTimeStepSize = EWOMS_GET_PARAM(TypeTag, Scalar, InitialTimeStepSize); + if (initialTimeStepSize < -1e50) { + if (myRank == 0) + Parameters::printUsage(argv[0], + "Mandatory parameter '--initial-time-step-size' " + "not specified!"); + return 1; + } + + if (myRank == 0) { +#ifdef EWOMS_VERSION + std::string versionString = EWOMS_VERSION; +#else + std::string versionString = ""; +#endif + const std::string briefDescription = Problem::briefDescription(); + if (!briefDescription.empty()) { + std::string tmp = Parameters::breakLines_(briefDescription, + /*indentWidth=*/0, + Parameters::getTtyWidth_()); + std::cout << tmp << std::endl << std::endl; + } + else + std::cout << "eWoms " << versionString + << " will now start the trip. " + << "Please sit back, relax and enjoy the ride.\n" + << std::flush; + } + + // print the parameters if requested + int printParams = EWOMS_GET_PARAM(TypeTag, int, PrintParameters); + if (myRank == 0) { + std::string endParametersSeparator("# [end of parameters]\n"); + if (printParams) { + bool printSeparator = false; + if (printParams == 1 || !isatty(fileno(stdout))) { + Opm::Parameters::printValues(); + printSeparator = true; + } + else + // always print the list of specified but unused parameters + printSeparator = + printSeparator || + Opm::Parameters::printUnused(); + if (printSeparator) + std::cout << endParametersSeparator; + } + else + // always print the list of specified but unused parameters + if (Opm::Parameters::printUnused()) + std::cout << endParametersSeparator; + } + + // print the properties if requested + int printProps = EWOMS_GET_PARAM(TypeTag, int, PrintProperties); + if (printProps && myRank == 0) { + if (printProps == 1 || !isatty(fileno(stdout))) + Opm::Properties::printValues(); + } + + // instantiate and run the concrete problem. make sure to + // deallocate the problem and before the time manager and the + // grid + Simulator simulator; + simulator.run(); + + if (myRank == 0) { + std::cout << "eWoms reached the destination. If it is not the one that was intended, " + << "change the booking and try again.\n" + << std::flush; + } + return 0; + } +#if ! DUNE_VERSION_NEWER(DUNE_COMMON, 2,5) + catch (Dune::Exception& e) + { + if (myRank == 0) { + std::cout << "Dune reported an error: " << e.what() << std::endl << std::flush; + + std::cout << "Trying to reset TTY.\n"; + resetTerminal_(); + } + + return 2; + } +#endif + catch (std::exception& e) + { + if (myRank == 0) { + std::cout << e.what() << ". Abort!\n" << std::flush; + + std::cout << "Trying to reset TTY.\n"; + resetTerminal_(); + } + + return 1; + } + catch (...) + { + if (myRank == 0) { + std::cout << "Unknown exception thrown!\n" << std::flush; + + std::cout << "Trying to reset TTY.\n"; + resetTerminal_(); + } + + return 3; + } +} + +} // namespace Opm + +#endif diff --git a/opm/models/utils/timer.hh b/opm/models/utils/timer.hh new file mode 100644 index 000000000..855c88dd8 --- /dev/null +++ b/opm/models/utils/timer.hh @@ -0,0 +1,217 @@ +// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- +// vi: set et ts=4 sw=4 sts=4: +/* + 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 2 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 . + + Consult the COPYING file in the top-level source directory of this + module for the precise wording of the license and the list of + copyright holders. +*/ +/*! + * \file + * + * \copydoc Opm::Timer + */ +#ifndef EWOMS_TIMER_HH +#define EWOMS_TIMER_HH + +#include + +#if HAVE_MPI +#include +#endif + +namespace Opm { +/*! + * \ingroup Common + * + * \brief Provides an encapsulation to measure the system time + * + * This means the wall clock time used by the simulation, the CPU time + * used by all threads of a single process and the CPU time used by + * the overall simulation. (i.e., the time used by all threads of all + * involved processes.) + */ +class Timer +{ + struct TimeData + { + std::chrono::high_resolution_clock::time_point realtimeData; + std::clock_t cputimeData; + }; +public: + Timer() + { halt(); } + + /*! + * \brief Start counting the time resources used by the simulation. + */ + void start() + { + isStopped_ = false; + measure_(startTime_); + } + + /*! + * \brief Stop counting the time resources. + * + * Returns the wall clock time the timer was active. + */ + double stop() + { + if (!isStopped_) { + TimeData stopTime; + + measure_(stopTime); + + const auto& t1 = startTime_.realtimeData; + const auto& t2 = stopTime.realtimeData; + std::chrono::duration dt = + std::chrono::duration_cast >(t2 - t1); + + realTimeElapsed_ += dt.count(); + cpuTimeElapsed_ += + static_cast(stopTime.cputimeData + - startTime_.cputimeData)/CLOCKS_PER_SEC; + } + + isStopped_ = true; + + return realTimeElapsed_; + } + + /*! + * \brief Stop the measurement reset all timing values + */ + void halt() + { + isStopped_ = true; + cpuTimeElapsed_ = 0.0; + realTimeElapsed_ = 0.0; + } + + /*! + * \brief Make the current point in time t=0 but do not change the status of the timer. + */ + void reset() + { + cpuTimeElapsed_ = 0.0; + realTimeElapsed_ = 0.0; + + measure_(startTime_); + } + + /*! + * \brief Return the real time [s] elapsed during the periods the timer was active + * since the last reset. + */ + double realTimeElapsed() const + { + if (isStopped_) + return realTimeElapsed_; + + TimeData stopTime; + + measure_(stopTime); + + const auto& t1 = startTime_.realtimeData; + const auto& t2 = stopTime.realtimeData; + std::chrono::duration dt = + std::chrono::duration_cast >(t2 - t1); + + return realTimeElapsed_ + dt.count(); + } + + /*! + * \brief This is an alias for realTimeElapsed() + * + * Its main purpose is to make the API of the class a superset of Dune::Timer + */ + double elapsed() const + { return realTimeElapsed(); } + + /*! + * \brief Return the CPU time [s] used by all threads of the local process for the + * periods the timer was active + */ + double cpuTimeElapsed() const + { + if (isStopped_) + return cpuTimeElapsed_; + + TimeData stopTime; + + measure_(stopTime); + + const auto& t1 = startTime_.cputimeData; + const auto& t2 = stopTime.cputimeData; + + return cpuTimeElapsed_ + static_cast(t2 - t1)/CLOCKS_PER_SEC; + } + + /*! + * \brief Return the CPU time [s] used by all threads of the all processes of program + * + * The value returned only differs from cpuTimeElapsed() if MPI is used. + */ + double globalCpuTimeElapsed() const + { + double val = cpuTimeElapsed(); + double globalVal = val; + +#if HAVE_MPI + MPI_Reduce(&val, + &globalVal, + /*count=*/1, + MPI_DOUBLE, + MPI_SUM, + /*rootRank=*/0, + MPI_COMM_WORLD); +#endif + + return globalVal; + } + + /*! + * \brief Adds the time of another timer to the current one + */ + Timer& operator+=(const Timer& other) + { + realTimeElapsed_ += other.realTimeElapsed(); + cpuTimeElapsed_ += other.cpuTimeElapsed(); + + return *this; + } + +private: + // measure the current time and put it into the object passed via + // the argument. + static void measure_(TimeData& timeData) + { + // Note: On Linux -- or rather fully POSIX compliant systems -- using + // clock_gettime() would be more accurate for the CPU time. + timeData.realtimeData = std::chrono::high_resolution_clock::now(); + timeData.cputimeData = std::clock(); + } + + bool isStopped_; + double cpuTimeElapsed_; + double realTimeElapsed_; + TimeData startTime_; +}; +} // namespace Opm + +#endif diff --git a/opm/models/utils/timerguard.hh b/opm/models/utils/timerguard.hh new file mode 100644 index 000000000..9b08bbc1e --- /dev/null +++ b/opm/models/utils/timerguard.hh @@ -0,0 +1,58 @@ +// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- +// vi: set et ts=4 sw=4 sts=4: +/* + 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 2 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 . + + Consult the COPYING file in the top-level source directory of this + module for the precise wording of the license and the list of + copyright holders. +*/ +/*! + * \file + * + * \copydoc Opm::TimerGuard + */ +#ifndef EWOMS_TIMER_GUARD_HH +#define EWOMS_TIMER_GUARD_HH + +#include "timer.hh" + +namespace Opm { +/*! + * \ingroup Common + * + * \brief A simple class which makes sure that a timer gets stopped if an exception is + * thrown. + */ +class TimerGuard +{ +public: + TimerGuard(Timer& timer) + : timer_(timer) + { } + + ~TimerGuard() + { + timer_.stop(); + } + +private: + Timer& timer_; +}; + +} // namespace Opm + +#endif diff --git a/opm/simulators/linalg/bicgstabsolver.hh b/opm/simulators/linalg/bicgstabsolver.hh index c0afc139f..edfb217b7 100644 --- a/opm/simulators/linalg/bicgstabsolver.hh +++ b/opm/simulators/linalg/bicgstabsolver.hh @@ -31,8 +31,8 @@ #include "residreductioncriterion.hh" #include "linearsolverreport.hh" -#include -#include +#include +#include #include diff --git a/opm/simulators/linalg/istlpreconditionerwrappers.hh b/opm/simulators/linalg/istlpreconditionerwrappers.hh index dcbbe623e..87e3cd359 100644 --- a/opm/simulators/linalg/istlpreconditionerwrappers.hh +++ b/opm/simulators/linalg/istlpreconditionerwrappers.hh @@ -43,8 +43,8 @@ #ifndef EWOMS_ISTL_PRECONDITIONER_WRAPPERS_HH #define EWOMS_ISTL_PRECONDITIONER_WRAPPERS_HH -#include -#include +#include +#include #include diff --git a/opm/simulators/linalg/istlsolverwrappers.hh b/opm/simulators/linalg/istlsolverwrappers.hh index 701a191e8..7e78f1768 100644 --- a/opm/simulators/linalg/istlsolverwrappers.hh +++ b/opm/simulators/linalg/istlsolverwrappers.hh @@ -43,8 +43,8 @@ #ifndef EWOMS_ISTL_SOLVER_WRAPPERS_HH #define EWOMS_ISTL_SOLVER_WRAPPERS_HH -#include -#include +#include +#include #include diff --git a/opm/simulators/linalg/linearsolverreport.hh b/opm/simulators/linalg/linearsolverreport.hh index ac8c34521..a8961e33b 100644 --- a/opm/simulators/linalg/linearsolverreport.hh +++ b/opm/simulators/linalg/linearsolverreport.hh @@ -29,8 +29,8 @@ #include "convergencecriterion.hh" -#include -#include +#include +#include namespace Opm { namespace Linear { diff --git a/opm/simulators/linalg/parallelbasebackend.hh b/opm/simulators/linalg/parallelbasebackend.hh index a74c3a3ea..300812b0c 100644 --- a/opm/simulators/linalg/parallelbasebackend.hh +++ b/opm/simulators/linalg/parallelbasebackend.hh @@ -36,9 +36,9 @@ #include #include -#include -#include -#include +#include +#include +#include #include #include diff --git a/opm/simulators/linalg/superlubackend.hh b/opm/simulators/linalg/superlubackend.hh index e3b5bac70..008547d3c 100644 --- a/opm/simulators/linalg/superlubackend.hh +++ b/opm/simulators/linalg/superlubackend.hh @@ -30,7 +30,7 @@ #if HAVE_SUPERLU #include -#include +#include #include diff --git a/tests/models/test_propertysystem.cpp b/tests/models/test_propertysystem.cpp index d6933b407..3d7cc76ac 100644 --- a/tests/models/test_propertysystem.cpp +++ b/tests/models/test_propertysystem.cpp @@ -31,7 +31,7 @@ */ #include "config.h" -#include +#include #include