diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml index e3268b6c8e1ae750f3e8cb3e821e84067ce3a83e..aaa284b1a04209e1f27303a1b43670e6c8ba5dbe 100644 --- a/.gitlab-ci.yml +++ b/.gitlab-ci.yml @@ -10,7 +10,7 @@ variables: BASE_IMAGE: dorie/dune-env # Use semantic versioning (not the version of DUNE) and bump according to # to whether changes are backwards-compatible or not. - IMAGE_VERSION: "1.1" + IMAGE_VERSION: "1.2" DUNE_ENV_IMAGE: ${BASE_IMAGE}:img-v${IMAGE_VERSION} CMAKE_FLAGS: @@ -69,6 +69,7 @@ setup:dune-env-clang: prep:update-dune: &update <<: *setup stage: prep + allow_failure: true only: - master - tags @@ -97,6 +98,7 @@ build:system-tests: &build-tests script: - CMAKE_FLAGS="$CMAKE_FLAGS" $DUNECONTROL --only=dorie configure + - $DUNECONTROL --only=dorie make $MAKE_FLAGS dorie-rfg - $DUNECONTROL --only=dorie make $MAKE_FLAGS build_system_tests - $DUNECONTROL --only=dorie make doc artifacts: @@ -122,6 +124,7 @@ build:debug: &debug -DCMAKE_BUILD_TYPE=Debug -DCMAKE_CXX_FLAGS_DEBUG='-Werror'" $DUNECONTROL --only=dorie configure + - $DUNECONTROL --only=dorie make $MAKE_FLAGS dorie-rfg - $DUNECONTROL --only=dorie make $MAKE_FLAGS build_unit_tests build:debug-clang: @@ -135,6 +138,7 @@ build:debug-clang: -DCMAKE_CXX_COMPILER=clang++ -DCMAKE_CXX_FLAGS_DEBUG='-Werror'" $DUNECONTROL --only=dorie configure + - $DUNECONTROL --only=dorie make $MAKE_FLAGS dorie-rfg - $DUNECONTROL --only=dorie make $MAKE_FLAGS build_unit_tests diff --git a/CHANGELOG.md b/CHANGELOG.md index d39a614ef160a8991b37f44511e0b618335b3357..4dfbf78e3ff345779c5e3c1f3526272c2c41d208 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -40,6 +40,8 @@ * Coupling between transient models for water flow and solute transport !96 * Initial conditions generated from H5 input data !130 * Generic Python VTK file reader !143 +* Parameterizations for hydrodynamic dispersion in solute transport !141 +* Generic Python VTK file reader !143, !150 ### Changed * `Simulation` is renamed `RichardsSimulation` and moved to @@ -143,6 +145,8 @@ * Specifying scaling field `extensions` and `offset` is now optional !133 * Generalized initial condition specification in config file !129 * Structure and setup of Sphinx user docs !126 +* Switch to stable `dune-randomfield` release branch !151, !153 +* System tests for executing `dorie pfg` module !153 ### Fixed * Solver in `RichardsSimulation` was using the wrong time variable. diff --git a/CMakeLists.txt b/CMakeLists.txt index 9ccafe33c5f6a561d15a4dfc81822829937e49d3..c080a3315461e30d4883e100c3bb063e5be67aa6 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -1,5 +1,5 @@ -cmake_minimum_required(VERSION 2.8.12) -project(dorie CXX) +cmake_minimum_required(VERSION 3.10) +project(dorie C CXX) # set build type if(NOT CMAKE_BUILD_TYPE) diff --git a/README.md b/README.md index e6ebb46e0b3dd163cd23ee52702a633cb57d238b..b57487a8f41821db9754a563d701f2ec5a30950a 100644 --- a/README.md +++ b/README.md @@ -88,7 +88,7 @@ by CI tests. | [dune-functions](https://gitlab.dune-project.org/staging/dune-functions) | releases/2.6 | [dune-typetree](https://gitlab.dune-project.org/staging/dune-typetree) | releases/2.6 | [dune-pdelab](https://gitlab.dune-project.org/pdelab/dune-pdelab) | releases/2.6 -| [dune-randomfield](https://gitlab.dune-project.org/oklein/dune-randomfield) | master +| [dune-randomfield](https://gitlab.dune-project.org/oklein/dune-randomfield) | releases/2.6 #### Optional Packages @@ -103,6 +103,9 @@ by CI tests. #### Step-by-step Instructions These instructions are suitable for a clean **Ubuntu** or **macOS** setup. The main difference between the two systems is the package manager. Debian-based systems have the APT manager already built in. On Mac, we recommend installing [Homebrew](http://brew.sh/). If you prefer to use [MacPorts](https://www.macports.org/), you need to check if some of the packages require different installation options than displayed here. +Manual installations on macOS require installing HDF5 from source. This can +be tricky, but the following instructions should work on a clean system. + If you installed [Anaconda](https://conda.io/docs/user-guide/install/download.html) on your machine, you don't need to install Python or Pip. Simply skip these packages when using the package managers for installing the software. However, notice the warnings when compiling DORiE below! 1. **macOS** users need to start by installing the Apple Command Line Tools by executing @@ -125,10 +128,33 @@ If you installed [Anaconda](https://conda.io/docs/user-guide/install/download.ht **macOS:** brew update - brew install cmake doxygen gcc libpng open-mpi muparser \ + brew install cmake doxygen fftw gcc libpng open-mpi muparser \ pkg-config python3 superlu yaml-cpp - brew install hdf5 --with-mpi - brew install fftw --with-mpi + +2. **macOS only:** Install HDF5 with MPI support from source. + + 1. Download an archive of the + [HDF5 source code](https://www.hdfgroup.org/downloads/hdf5/source-code/), + and extract it. + + 2. Enter the extracted folder. In there, create a `build` directory, and + enter it: + + mkdir build && cd build + + 3. Configure your build. If you followed the instructions above, the + OpenMPI C compiler is reachable via the command `mpicc`. If not, you have + to specify a full path to it. Use the option `prefix` to specify where + you want the package to be installed. This should *not* be a + system-reserved path like `/usr/local`, and *not* be located in a + sub-directory of the source code. Execute the configuration script: + + ./../configure CC=mpicc --prefix= --enable-parallel + + 4. Build and install the library: + + make && make install + 3. The parallel linear solver of DORiE can make use of the ParMETIS package. If you want to run DORiE in parallel on multiple processes, additionally install METIS and ParMETIS: @@ -158,13 +184,21 @@ If you installed [Anaconda](https://conda.io/docs/user-guide/install/download.ht If you installed software into paths not appended to your `PATH` variable, you will have to add `CMAKE_FLAGS` to the call to make sure that CMake finds all packages. Alternatively, you can add a custom options file. See the [DUNE Installation Instructions](https://dune-project.org/doc/installation/) for details. CMake will throw an error if required packages are not found. - **Warning:** Anacoda supplies its own version of HDF5 which is typically found first by CMake. If you have Anaconda installed on your machine, add + If you installed HDF5 from source (all macOS users) or use Anaconda, + specify the path to your HDF5 installation by using the `HDF5_ROOT` + variable. On Ubuntu, add the path to the APT package, + + -DHDF5_ROOT=/usr/ + + and on macOS, add + + -DHDF5_ROOT= - HDF5_ROOT=/usr/local + to the `CMAKE_FLAGS` in the above command, replacing `` with the + path chosen as installation prefix when configuring HDF5. - to the `CMAKE_FLAGS` in the call to `dunecontrol` above. -#### Experimental Features +### Experimental Features The local operator implementing Richards equation's discretization supports multiple scheme settings. Setting these via the config file is disabled by default. You can enable this feature by reconfiguring DORiE with the CMake flag diff --git a/cmake/modules/DorieMacros.cmake b/cmake/modules/DorieMacros.cmake index b311f50a87d214742869b843b6303eab75571c2e..922b957091f206c229081f2292d78915f64de965 100644 --- a/cmake/modules/DorieMacros.cmake +++ b/cmake/modules/DorieMacros.cmake @@ -1,29 +1,51 @@ -# find all required packages -FIND_PACKAGE (HDF5 REQUIRED) +# These macros check for the following packages, yielding the respective +# targets +# +# HDF5 ... hdf5 +# yaml-cpp ... yaml-cpp +# muparser ... muparser::muparser +# +# All of these are required an will be autmatically linked to executables when +# using the DorieTesting CMake macros. + +# Add CMake policy stack +cmake_policy(PUSH) + +# Allow using `_ROOT` on CMake >3.12 +if(POLICY CMP0074) + cmake_policy(SET CMP0074 NEW) +endif() + +# HDF5 +set(HDF5_PREFER_PARALLEL TRUE) +find_package (HDF5 REQUIRED COMPONENTS C) if(NOT HDF5_IS_PARALLEL) message(SEND_ERROR "Parallel HDF5 must be installed!") endif() -add_definitions(-DHDF5_PARALLEL) -FIND_PACKAGE (FFTW REQUIRED) -FIND_PACKAGE (SuperLU REQUIRED) -FIND_PACKAGE (MPI REQUIRED) +# Define a target because the CMake module does not +add_library(hdf5 INTERFACE IMPORTED GLOBAL) + +# TODO: Use `target_XXX` functions after raising CMake requirements +# Sanitize the definitions because we cannot use the proper CMake function... +string(REPLACE -D "" HDF5_DEFINITIONS "${HDF5_DEFINITIONS}") +# Set properties directly because of a bug in CMake 3.10 +set_target_properties(hdf5 PROPERTIES + INTERFACE_LINK_LIBRARIES "${HDF5_LIBRARIES}" + INTERFACE_INCLUDE_DIRECTORIES "${HDF5_INCLUDE_DIRS}" + INTERFACE_COMPILE_DEFINITIONS "${HDF5_DEFINITIONS};H5_HAVE_PARALLEL") + +# yaml-cpp find_package (yaml-cpp 0.5.2 REQUIRED) -FIND_PACKAGE (muparser REQUIRED) -FIND_PACKAGE (METIS) -FIND_PACKAGE (ParMETIS) +# muparser +find_package (muparser REQUIRED) -include_directories(${FFTW_INCLUDES} - ${HDF5_INCLUDE_DIRS} - ${YAML_CPP_INCLUDE_DIR}) +# Report the DUNE libs +message (STATUS "DUNE Libraries: ${DUNE_LIBS}") -list (APPEND DUNE_LIBS - ${FFTW_LIBRARIES} - ${HDF5_LIBRARIES} - ${YAML_CPP_LIBRARIES} - muparser::muparser) -MESSAGE(STATUS "DUNE Libraries: ${DUNE_LIBS}") +# Remove CMake policy stack +cmake_policy(POP) # Add DORiE testing functions include(DorieTesting) diff --git a/cmake/modules/DorieTesting.cmake b/cmake/modules/DorieTesting.cmake index 92faac4ffe6f0cac0e2843b71405683576ce1c01..2be1acc42f70725c8c0dd49006a88d2966fedd9b 100644 --- a/cmake/modules/DorieTesting.cmake +++ b/cmake/modules/DorieTesting.cmake @@ -7,7 +7,8 @@ add_custom_target(unit_tests # Add the system test build target and test command add_custom_target(build_system_tests) add_custom_target(system_tests - COMMAND ctest --output-on-failure --exclude-regex ^ut.+$ + COMMAND ctest + --output-on-failure --exclude-regex "^\\(ut.+\\|python_tests\\)$" ) # create the mapping datafile as preparation for all tests @@ -53,6 +54,9 @@ endfunction() # tests within DORiE and adds flags for coverage reports. Notice that # `dune_add_test` requires more parameters than this function alone. # +# The registered executable is automatically linked to the libraries DORiE +# DORiE depends on. +# # Use this function exactly like `dune_add_test`. # function(dorie_add_unit_test) @@ -75,6 +79,8 @@ function(dorie_add_unit_test) endif() # add to build target and employ compile options + target_link_libraries(${UNIT_TEST_TARGET} + muparser::muparser hdf5 yaml-cpp spdlog) add_coverage_links(${UNIT_TEST_TARGET}) add_dependencies(build_unit_tests ${UNIT_TEST_TARGET}) endfunction() @@ -129,6 +135,9 @@ endfunction() # Use this function like `dune_add_system_test`, considering the options # given above. # +# When registering a new executable (`TARGET` is not specified), this +# executable will be linked to the libraries DORiE depends on. +# function(dorie_add_metaini_test) set(OPTIONS UNIT_TEST) set(SINGLE TARGET METAINI SCRIPT BASENAME CREATED_TARGETS) @@ -177,9 +186,16 @@ function(dorie_add_metaini_test) SCRIPT ${SYSTEM_TEST_SCRIPT} CREATED_TARGETS created_targets ) + # report created targets to parent scope set(${SYSTEM_TEST_CREATED_TARGETS} ${created_targets} PARENT_SCOPE) + # Link to dependencies + if(NOT SYSTEM_TEST_TARGET) + target_link_libraries(${created_targets} + muparser::muparser hdf5 yaml-cpp spdlog) + endif() + # add dependencies and flags if(SYSTEM_TEST_UNIT_TEST) add_coverage_links(${created_targets}) diff --git a/cmake/modules/FindFFTW.cmake b/cmake/modules/FindFFTW.cmake deleted file mode 100644 index f0415f8a79e3683af032ceb6f05491fd1830ed8f..0000000000000000000000000000000000000000 --- a/cmake/modules/FindFFTW.cmake +++ /dev/null @@ -1,114 +0,0 @@ -# - Find the FFTW library -# -# Usage: -# find_package(FFTW [REQUIRED] [QUIET] ) -# -# It sets the following variables: -# FFTW_FOUND ... true if fftw is found on the system -# FFTW_LIBRARIES ... full path to fftw library -# FFTW_INCLUDES ... fftw include directory -# -# The following variables will be checked by the function -# FFTW_USE_STATIC_LIBS ... if true, only static libraries are found -# FFTW_ROOT ... if set, the libraries are exclusively searched -# under this path -# FFTW_LIBRARY ... fftw library to use -# FFTW_INCLUDE_DIR ... fftw include directory -# - -#If environment variable FFTWDIR is specified, it has same effect as FFTW_ROOT -if( NOT FFTW_ROOT AND ENV{FFTWDIR} ) - set( FFTW_ROOT $ENV{FFTWDIR} ) -endif() - -# Check if we can use PkgConfig -find_package(PkgConfig) - -#Determine from PKG -if( PKG_CONFIG_FOUND AND NOT FFTW_ROOT ) - pkg_check_modules( PKG_FFTW QUIET "fftw3" ) -endif() - -#Check whether to search static or dynamic libs -set( CMAKE_FIND_LIBRARY_SUFFIXES_SAV ${CMAKE_FIND_LIBRARY_SUFFIXES} ) - -if( FFTW_ROOT ) - - #find libs - find_library( - FFTW_LIB - NAMES "fftw3" - PATHS ${FFTW_ROOT} - PATH_SUFFIXES "lib" "lib64" - NO_DEFAULT_PATH - ) - - find_library( - FFTWF_LIB - NAMES "fftw3_mpi" - PATHS ${FFTW_ROOT} - PATH_SUFFIXES "lib" "lib64" - NO_DEFAULT_PATH - ) - - find_library( - FFTWL_LIB - NAMES "fftw3l" - PATHS ${FFTW_ROOT} - PATH_SUFFIXES "lib" "lib64" - NO_DEFAULT_PATH - ) - - #find includes - find_path( - FFTW_INCLUDES - NAMES "fftw3.h" - PATHS ${FFTW_ROOT} - PATH_SUFFIXES "include" - NO_DEFAULT_PATH - ) - -else() - - find_library( - FFTW_LIB - NAMES "fftw3" - PATHS ${PKG_FFTW_LIBRARY_DIRS} ${LIB_INSTALL_DIR} - ) - - find_library( - FFTWF_LIB - NAMES "fftw3_mpi" - PATHS ${PKG_FFTW_LIBRARY_DIRS} ${LIB_INSTALL_DIR} - ) - - - find_library( - FFTWL_LIB - NAMES "fftw3l" - PATHS ${PKG_FFTW_LIBRARY_DIRS} ${LIB_INSTALL_DIR} - ) - - find_path( - FFTW_INCLUDES - NAMES "fftw3.h" - PATHS ${PKG_FFTW_INCLUDE_DIRS} ${INCLUDE_INSTALL_DIR} - ) - -endif( FFTW_ROOT ) - -set(FFTW_LIBRARIES ${FFTWF_LIB} ${FFTW_LIB}) - -if(FFTWL_LIB) - set(FFTW_LIBRARIES ${FFTWL_LIB} ${FFTW_LIBRARIES}) -endif() - -set( CMAKE_FIND_LIBRARY_SUFFIXES ${CMAKE_FIND_LIBRARY_SUFFIXES_SAV} ) - -MESSAGE(STATUS "Found FFTW library: ${FFTW_LIBRARIES}") - -include(FindPackageHandleStandardArgs) -find_package_handle_standard_args(FFTW DEFAULT_MSG - FFTW_INCLUDES FFTW_LIBRARIES) - -mark_as_advanced(FFTW_INCLUDES FFTW_LIBRARIES) diff --git a/cmake/modules/Findmuparser.cmake b/cmake/modules/Findmuparser.cmake index 32007f33f9162706ca63e0120570c2620020c4bc..d89d859fd81fa877671e0b7c1ddf471ca1f86dfa 100644 --- a/cmake/modules/Findmuparser.cmake +++ b/cmake/modules/Findmuparser.cmake @@ -1,44 +1,45 @@ -# - Find the muParser library +# Find the muparser library # # Usage: -# find_package(muParser [REQUIRED] [QUIET] ) +# find_package(muparser [REQUIRED] [QUIET] ) # # Hints may be given by the (environment) variables -# MUPARSER_ROOT ... Path to the installation library +# muparser_ROOT ... Path to the installation library # # It sets the following variables: -# MUPARSER_FOUND ... true if muParser is found on the system -# MUPARSER_LIBRARIES ... full path to muParser library -# MUPARSER_INCLUDES ... muParser include directory +# muparser_FOUND ... true if muparser is found on the system +# muparser_LIBRARIES ... full path to muparser library +# muparser_INCLUDES ... muparser include directory # # It defines the following targets: # muparser::muparser ... muparser library to link against # -find_library(MUPARSER_LIBRARY - NAMES muparser muparserd - HINTS ${MUPARSER_ROOT} ENV MUPARSER_ROOT +find_library(muparser_LIBRARY + NAMES muparser muparserd + HINTS ${muparser_ROOT} ENV muparser_ROOT ) -find_path(MUPARSER_INCLUDE_DIR - muParserDef.h - HINTS ${MUPARSER_ROOT} ENV MUPARSER_ROOT +find_path(muparser_INCLUDE_DIR + muParserDef.h + HINTS ${muparser_ROOT} ENV muparser_ROOT ) -mark_as_advanced(MUPARSER_INCLUDE_DIR MUPARSER_LIBRARY) +mark_as_advanced(muparser_INCLUDE_DIR muparser_LIBRARY) -find_package_handle_standard_args(MUPARSER - DEFAULT_MSG - MUPARSER_LIBRARY MUPARSER_INCLUDE_DIR) +include(FindPackageHandleStandardArgs) +find_package_handle_standard_args(muparser + DEFAULT_MSG + muparser_LIBRARY muparser_INCLUDE_DIR) -if (MUPARSER_FOUND) - set(MUPARSER_LIBRARIES ${MUPARSER_LIBRARY} ) - set(MUPARSER_INCLUDES ${MUPARSER_INCLUDE_DIR} ) +if (muparser_FOUND) + set(muparser_LIBRARIES ${muparser_LIBRARY} ) + set(muparser_INCLUDES ${muparser_INCLUDE_DIR} ) # add the target add_library(muparser::muparser SHARED IMPORTED) set_target_properties(muparser::muparser - PROPERTIES IMPORTED_LOCATION ${MUPARSER_LIBRARIES} - INTERACE_INCLUDE_DIRECTORIES ${MUPARSER_INCLUDES} + PROPERTIES IMPORTED_LOCATION ${muparser_LIBRARIES} + INTERACE_INCLUDE_DIRECTORIES ${muparser_INCLUDES} ) endif() diff --git a/doc/default_files/CMakeLists.txt b/doc/default_files/CMakeLists.txt index 44fa529704471b4434c41b54cb779a3c816695bd..0f9cd64d9f9f3a0ac09b49cbe369e457b85d74a7 100644 --- a/doc/default_files/CMakeLists.txt +++ b/doc/default_files/CMakeLists.txt @@ -119,7 +119,7 @@ endfunction() # copy BC & parameter files file(COPY 2d_infiltr.bcdat 3d_infiltr.bcdat 2d_solute.bcdat 3d_solute.bcdat - param.yml + richards_param.yml transport_param.yml DESTINATION .) # Random field generator diff --git a/doc/default_files/common-parameters.xml b/doc/default_files/common-parameters.xml index 2bdb54d8b62ddefae87e81c072ea2e94fe306ab8..151cbfd2cac13537dc232720f5820b0332bf6ee5 100644 --- a/doc/default_files/common-parameters.xml +++ b/doc/default_files/common-parameters.xml @@ -32,10 +32,15 @@ adding an empty line, make text **bold** or ``monospaced``. - Mode of simulation in DORiE. + Sets the simulation mode of DORiE. + (``richards``) mode calculates the matric head with a + DG scheme and produce water fluxes with Raviart Thomas reconstruction. + (``richards+transport``) mode calculates (``richards``) mode and use + the generated water flux and saturation at each step to calculate the + solute transport model for unsaturated media with a DG scheme. + richards richards, richards+transport - richards, richards+transport diff --git a/doc/default_files/param.yml b/doc/default_files/param.yml deleted file mode 100644 index b95a2c843e3dbd341874661297d8655098fc54e4..0000000000000000000000000000000000000000 --- a/doc/default_files/param.yml +++ /dev/null @@ -1,25 +0,0 @@ -volumes: - sand: - index: 0 - type: MvG - parameters: - alpha: -2.3 - n: 4.17 - k0: 2.2E-5 - theta_r: 0.03 - theta_s: 0.31 - tau: -1.1 - - silt: - index: 1 - type: MvG - parameters: - alpha: -0.7 - n: 1.3 - k0: 1.0E-5 - theta_r: 0.01 - theta_s: 0.41 - tau: 0.0 - -scaling: - type: None diff --git a/doc/default_files/richards-parameters.xml b/doc/default_files/richards-parameters.xml index 7a4834ffe0f27c59ccf6e60af1988b41ba34a5b0..43a0a8ef3e739f1eb565bd8c8804e9155459cac6 100644 --- a/doc/default_files/richards-parameters.xml +++ b/doc/default_files/richards-parameters.xml @@ -227,6 +227,7 @@ adding an empty line, make text **bold** or ``monospaced``. YAML file containing the parameter definitions. path + richards_param.yml diff --git a/doc/default_files/richards_param.yml b/doc/default_files/richards_param.yml new file mode 100644 index 0000000000000000000000000000000000000000..580a49bcf00824b7ed9f1eeffe26a1a48161cffe --- /dev/null +++ b/doc/default_files/richards_param.yml @@ -0,0 +1,27 @@ +volumes: + sand: + index: 0 + richards: + type: MvG + parameters: + alpha: -2.3 + n: 4.17 + k0: 2.2E-5 + theta_r: 0.03 + theta_s: 0.31 + tau: -1.1 + + silt: + index: 1 + richards: + type: MvG + parameters: + alpha: -0.7 + n: 1.3 + k0: 1.0E-5 + theta_r: 0.01 + theta_s: 0.41 + tau: 0.0 + +scaling: + type: None diff --git a/doc/default_files/transport-parameters.xml b/doc/default_files/transport-parameters.xml index 7a7c573125ba0840df70de2010869c375666c702..35d95f2fb775339c1a7f4b3e1d532c9c273972f7 100644 --- a/doc/default_files/transport-parameters.xml +++ b/doc/default_files/transport-parameters.xml @@ -67,6 +67,16 @@ adding an empty line, make text **bold** or ``monospaced``. endOfTransportStep, endOfRichardsStep, none + + Defines whether VTK files should include the hydrodynamic + dispersion tensor. Tensors are written in 3D and have 9 componentents + independently of the world dimension. This can be easily be visualizated + in Paraview with the ``Tensor Glyph`` filter. + + true, false + false + + Plot VTK files with virtually refined grids. VTK only supports bilinear triangulations and displays higher-order solutions @@ -109,6 +119,13 @@ adding an empty line, make text **bold** or ``monospaced``. true, false false + + + Type of the input value for the dirichlet condition. + + soluteConcentration, totalSolute + soluteConcentration + @@ -222,11 +239,12 @@ adding an empty line, make text **bold** or ``monospaced``. - - - Effective hydrodynamic dispersion - float > 0 - 2E-9 + + + YAML file containing the parameter definitions. + + path + transport_param.yml @@ -253,8 +271,7 @@ adding an empty line, make text **bold** or ``monospaced``. - Order of the finite element method used. Values > 1 are not - thoroughly tested. + Order of the finite element method used. 0 0 diff --git a/doc/default_files/transport_param.yml b/doc/default_files/transport_param.yml new file mode 100644 index 0000000000000000000000000000000000000000..08821c65334eca0be177bfe9f34aaa4f371da108 --- /dev/null +++ b/doc/default_files/transport_param.yml @@ -0,0 +1,22 @@ +volumes: + sand: + index: 0 + transport: + type: Deff_MQ1 + Dhm_iso + parameters: + char_length: 1.5E-11 + mol_diff: 2.E-9 + phi: 0.34 + lambda_t: 0.0025 + lambda_l: 0.025 + + silt: + index: 1 + transport: + type: Deff_MQ1 + Dhm_iso + parameters: + char_length: 1.5E-11 + mol_diff: 2.E-9 + phi: 0.34 + lambda_t: 0.0025 + lambda_l: 0.025 \ No newline at end of file diff --git a/doc/introduction/data-io.rst b/doc/introduction/data-io.rst index f57fe3900e02e08fb6dd7e1e0d8739605fc434fc..48ff078d3f1f9f16130e2d38342f93c9e4cfeaf5 100644 --- a/doc/introduction/data-io.rst +++ b/doc/introduction/data-io.rst @@ -64,11 +64,16 @@ Output Files .. object:: Transport - ``solute``: Solute concentration in water phase - :math:`c_w \, [\mathrm{kg}/\mathrm{m}^3]` + :math:`c_w \, [\mathrm{kg}\mathrm{m}^{-3}]` - ``solute_total``: Total solute concentration - :math:`c_t \, [\mathrm{kg}/\mathrm{m}^3]` + :math:`c_t \, [\mathrm{kg}\mathrm{m}^{-3}]` + - ``micro_peclet``: Microscopic peclet number + :math:`\mathsf{pe_\mu} \, [-]` - ``flux_RT``: Reconstructed solute flux - :math:`\mathbf{j}_{s, \mathrm{rc}} \, [\mathrm{kg}/\mathrm{m}^{-2}\mathrm{s}^{-1}]` + :math:`\mathbf{j}_{s, \mathrm{rc}} \, [\mathrm{kg}\,\mathrm{m}^{-2}\mathrm{s}^{-1}]` + (if enabled!) + - ``eff_hd_dispersion``: Hydrodynamic dispersion tensor + :math:`\mathsf{D}_\mathsf{hd} \, [\mathrm{m}^2\mathrm{s}^{-1}]` (if enabled!) * VTK Parallel Collection Files (``.pvtu``): Merging multiple VTK files in case diff --git a/doc/man-cookbook.rst b/doc/man-cookbook.rst index 8063eeb9f6e71c3df5e99dbbbade7def345f703f..e7aa0b59badab143628dff11d2d2dcb5bd1b86e3 100644 --- a/doc/man-cookbook.rst +++ b/doc/man-cookbook.rst @@ -7,7 +7,7 @@ Cook Book Introduction ============ -This part of the documentation is intended for first-time DORiE users. It explaines the basic usage of the program, how to execute a simulation and how to analyze its results. It will also showcase more complex features like Adaptive Grid Refinement. +This part of the documentation is intended for first-time DORiE users. It explains the basic usage of the program, how to execute a simulation and how to analyze its results. It will also showcase more complex features like Adaptive Grid Refinement. Prerequisites ------------- @@ -57,7 +57,7 @@ Let's begin with the ``output``. We want the program to give us at least some ou .. code-block:: none [output] - verbose = 1 + logLevel = info outputPath = ./sand fileName = sand diff --git a/doc/manual/parameter-file.rst b/doc/manual/parameter-file.rst index b609fc24cd2df5e6e3fdab7d2496249f2cd01370..7153dd7237f2d62da93fbcf2ab9128723e80b24a 100644 --- a/doc/manual/parameter-file.rst +++ b/doc/manual/parameter-file.rst @@ -4,11 +4,16 @@ Soil Parameterization Specifying the domain properties requires input of the :doc:`soil architecture ` (in relation to the grid) and of the soil parameters. The latter incorporate -the *subscale physics*, the representation of the soil hydraulic properties -below the REV scale. DORiE requires several input files for retrieving this -information, depending on the type of grid used for the computation. +the *subscale physics*, the representation of the soil hydraulic and solute +properties below the REV scale. DORiE requires several input files for +retrieving this information, depending on the type of grid used for the +computation. -A YAML_ parameter file is always required. It specifies a set of soil +.. contents:: + :depth: 3 + :local: + +A YAML_ parameter file is always required. It specifies a set of parameterizations, along with a medium index. This index is used to reference the medium specified in the architecture files, or with the "global" indices given in the config file. @@ -24,13 +29,18 @@ YAML Parameter File ------------------- This file is used to specify the parameterization for each medium inside the simulated domain. It follows a simple hierarchical syntax. The file path and -name must be specified via the ``parameters.file`` key of the +name must be specified via the ``.parameters.file`` key of the :doc:`config file `. The top-level mapping must contain the key ``volumes``. This key contains a sequence of arbitrary parameterization names. These, in turn, contain the -parameterization ``type``, the medium ``index``, and the actual ``parameters``. -The medium index must be an **integer**. +medium ``index``, and the model type (``richards`` or ``transport``). +The medium index must be an **integer**. Each model key contains the +parameterization ``type``, and the actual ``parameters``. + +.. note:: Parameterization data for model ``transport`` is only required if + the model is actually enabled in the :doc:`config file settings + `, via ``simulation.mode = richards+transport``. Heterogeneities throughout the entire domain are set via the top level key ``scaling``. It must contain a supported scaling ``type``, which may default @@ -38,6 +48,10 @@ to ``none``. In case an actual scaling is to be applied, the section must contain the key ``data``, which specifies the datasets required for this type of scaling. +.. warning:: + Transport parameters like `porosity` or `characteristic length` currently + are not affected by small scale heterogoeneities. + Every ``scaling_section`` has the same keys: The H5 filepath is given by ``file``, and the internal dataset path by ``dataset``. You can then choose an ``interpolation`` method for the dataset, which requires you to insert values @@ -51,11 +65,19 @@ shapes. volumes: : - type: index: - parameters: - : - # more parameters ... + richards: + type: + parameters: + : + # more parameters ... + + # if transport mode is enabled ... + transport: + type: + parameters: + : + # more parameters ... # more volumes ... @@ -74,9 +96,9 @@ shapes. You find a documentation of the parameterizations implemented in DORiE along with the parameters they require below. -The parameterization ``type`` must match the static ``type`` member of one of -them. Likewise, the parameters must match the names of parameters -these objects use. +The parameterization ``type`` must match a static parameterization name or a +valid combination of them. Likewise, the parameters must match the names of +parameters these objects use. .. _man-parameter_scaling: @@ -101,8 +123,8 @@ This section lists the available types for parameterizations, scalings, and interpolators in a understandable format. You can also dig through the code documentation below, or the code itself, to retrieve this information. -Parameterizations -^^^^^^^^^^^^^^^^^ +Richards Parameterizations +~~~~~~~~~~~~~~~~~~~~~~~~~~ As the *soil porosity* :math:`\phi \, [-]` and the *residual water content* :math:`\theta_s \, [-]` vary between soil types, it is convenient to define the @@ -111,13 +133,15 @@ As the *soil porosity* :math:`\phi \, [-]` and the *residual water content* .. math:: \Theta = \frac{\theta_w - \theta_r}{\phi - \theta_r}, - \quad \Theta \in \left[ 0, 1 \right] + \quad \Theta \in \left[ 0, 1 \right], +where :math:`\theta_w \, [-]` is the volumetric soil water content. One typically assumes that the entire pore space can be filled up with water, hence we set the *saturated water content* :math:`\theta_s \, [-]` equal to the porosity, :math:`\theta_s = \phi`. -.. object:: Mualem–van Genuchten parameterization +Mualem–van Genuchten +"""""""""""""""""""" Implements the following parameterization functions: @@ -136,8 +160,8 @@ porosity, :math:`\theta_s = \phi`. Parameters: * ``theta_r``: Residual water content :math:`\theta_r` * ``theta_s``: Saturated water content / Porosity :math:`\theta_s` - * ``k0``: Saturated conductivity :math:`K_0 \, [\mathrm{ms}^{-1}]` - * ``alpha``: Air-entry value :math:`\alpha \, [\mathrm{m}^{-1}]` + * ``k0``: Saturated conductivity :math:`K_0 \, [\text{ms}^{-1}]` + * ``alpha``: Air-entry value :math:`\alpha \, [\text{m}^{-1}]` * ``n``: Retention curve shape factor :math:`n` * ``tau``: Skew factor :math:`\tau` @@ -145,6 +169,7 @@ porosity, :math:`\theta_s = \phi`. .. code-block:: yaml + richards: type: MvG parameters: theta_r: @@ -154,8 +179,323 @@ porosity, :math:`\theta_s = \phi`. n: tau: +Transport Parameterizations +~~~~~~~~~~~~~~~~~~~~~~~~~~~ + +Regardless of the parameterization, the transport simulation always computes +the microscopic péclet number, for which it requires the characteristic pore +length, the molecular diffusion, and the fluid velocity. The latter is directly +provided by the richards simulation while the other two have to be specified +for each volume: + +Permanent parameters: + * ``mol_diff``: Molecular diffusion + :math:`D_m \, [\text{m}^2\text{s}^{-1}]` + * ``char_length``: Characteristic pore length :math:`\ell \, [\text{m}]` + +.. note:: + We support two options for specifying tensors through the YAML syntax. + You may either specify every entry of the tensor with a dedicated key, like + + .. code-block:: yaml + + _xx: + _xy: + # ... + + or give an entire sequence that will be mapped to the respective entries, + + .. code-block:: yaml + + : [, # , ... + ] + + The sequence is interpreted as a flattened tensor with row-major layout. + +.. warning:: + You must omit any component containing the spatial direction ``z`` in a + 2D setup. + +The hydrodynamic dispersion tensor :math:`\mathsf{D}_\text{hd} \, +[\text{m}^2\text{s}^{-1}]` is the main parameter to provide in the +transport simulation. Below you will find several parameterization for this. + +Constant +"""""""" + + In this case, the hydrodynamic dispersion tensor is inserted directly + component-by-compoment. + + .. note:: + From a physical point of view, the hydrodynamic tensor must be + symmetric, but the user input is not verified by DORiE in this regard. + + .. math:: + + \mathsf{D}_\text{hd} = \text{const}. + + * ``type: Dhd_const`` + + Parameters: + * ``hydrodynamic_disp_``: (i, j)-th component of the + hydrodynamic dispersion tensor, + :math:`\left( \mathsf{D}_\text{hd} \right)_{ij} \, [\text{m}^2\text{s}^{-1}]`, + **or** + * ``hydrodynamic_disp``: Flattened hydrodynamic dispersion tensor + :math:`\mathsf{D}_\text{hd} \, [\text{m}^2\text{s}^{-1}]`. + + YAML template: + + .. code-block:: yaml + + transport: + type: Dhd_const + parameters: + mol_diff: + char_length: + + # sequence notation, or... + hydrodynamic_disp: [ ] + + # component notation + hydrodynamic_disp_xx: + hydrodynamic_disp_xy: + hydrodynamic_disp_xz: + hydrodynamic_disp_yx: + hydrodynamic_disp_yy: + hydrodynamic_disp_yz: + hydrodynamic_disp_zx: + hydrodynamic_disp_zy: + hydrodynamic_disp_zz: + +Power Law +""""""""" + Implements the following parameterization function: + + .. math:: + + D_\text{hd} = \gamma D_m \operatorname{pe}^\alpha. + + * ``type: Dhd_pl`` + + Parameters: + * ``gamma``: Scale the power law :math:`\gamma \, [-]` + * ``alpha``: Exponent of the power law :math:`\alpha \, [-]` + * ``mol_diff``: Molecular diffusion + :math:`D_m \, [\text{m}^2\text{s}^{-1}]` + + The Peclét number :math:`\operatorname{pe}` is specified in the + :doc:`config file `. + + YAML template: + + .. code-block:: yaml + + transport: + type: Dhd_pl + parameters: + mol_diff: + char_length: + alpha: + gamma: + +Superposition +""""""""""""" + + The hydrodynamic dispersion tensor is said to be the superposition of + several diffusion/dispersion processes. In DORiE this possible by summing + several valid parameterizations types. Currently we provide + parameterizations for the *effective diffusion* + :math:`D_\text{eff}` and for the *effective hydromechanic tensor* + :math:`\mathsf{D}_\text{hm}` identified by the key prefixes ``Deff`` and + ``Dhm`` respectively. + + .. math:: + + \mathsf{D}_\text{hd} = \mathsf{D}_\text{hm} + D_\text{eff}. + + * ``type: + `` + + Each of the types are further parameterized and can be found below. + +Effective Diffusion +^^^^^^^^^^^^^^^^^^^ + +Constant Effective Diffusion +############################ + + In this case, the effective diffusion is inserted directly, + + .. math:: + + D_\text{eff} = \text{const}. + + * ``Deff_type: Deff_const`` + + Parameters: + + * ``eff_diff``: Effective diffusion + :math:`D_\text{eff} \, [\text{m}^2\text{s}^{-1}]` + + YAML template: + + .. code-block:: yaml + + transport: + type: + Deff_const + parameters: + mol_diff: + char_length: + eff_diff: + # parameters ... + +Milligton-Quirk I Effective Diffusion +##################################### + + Implements the following parameterization function: + + .. math:: + + D_\text{eff} = D_m \frac{\theta_w^{7/3}}{\phi^{2/3}}. + + where the volumetric water content :math:`\theta_w \, [-]` + is automatically obtained from the Richards simulation. + + * ``Deff_type: Deff_MQ1`` + + Parameters: + * ``mol_diff``: Molecular diffusion + :math:`D_m \, [\text{m}^2\text{s}^{-1}]` + * ``phi``: Soil porosity :math:`\phi \, [-]` + + YAML template: + + .. code-block:: yaml + + transport: + type: + Deff_MQ1 + parameters: + mol_diff: + char_length: + phi: + # parameters ... + +Milligton-Quirk II Effective Diffusion +###################################### + + Implements the following parameterization function: + + .. math:: + + D_\text{eff} = D_m \frac{\theta_w}{\phi^{2/3}}. + + where the volumetric water content :math:`\theta_w \, [-]` + is automatically obtained from the Richards simulation. + + * ``Deff_type: Deff_MQ2`` + + Parameters: + * ``mol_diff``: Molecular diffusion + :math:`D_m \, [\text{m}^2\text{s}^{-1}]` + * ``phi``: Soil porosity :math:`\phi \, [-]` + + YAML template: + + .. code-block:: yaml + + transport: + type: + Deff_MQ1 + parameters: + mol_diff: + char_length: + phi: + # parameters ... + +Effective Hydromechanic Dispersion +^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ + +Constant Effective Hydromechanic Dispersion Tensor +################################################## + + In this case, the effective hydromechanic dispersion tensor is inserted + directly. + + .. math:: + + \mathsf{D}_\text{hm} = \text{const}. + + * ``Dhm_type: Dhm_const`` + + Parameters: + * ``eff_hydromechanic_disp_``: (i, j)-th component of the + hydromechanic dispersion tensor, + :math:`\left(\mathsf{D}_\text{hm}\right)_{ij} \, [\text{m}^2\text{s}^{-1}]`, + **or** + * ``eff_hydromechanic_disp``: Flattened hydromechanic dispersion tensor + :math:`\mathsf{D}_\text{hm} \, [\text{m}^2\text{s}^{-1}]`. + + YAML template: + + .. code-block:: yaml + + transport: + type: Dhm_const + + parameters: + mol_diff: + char_length: + + # sequence notation, or... + eff_hydromechanic_disp: [ ] + + # component notation + eff_hydromechanic_disp_xx: + eff_hydromechanic_disp_xy: + eff_hydromechanic_disp_xz: + eff_hydromechanic_disp_yx: + eff_hydromechanic_disp_yy: + eff_hydromechanic_disp_yz: + eff_hydromechanic_disp_zx: + eff_hydromechanic_disp_zy: + eff_hydromechanic_disp_zz: + + # parameters ... + +Effective Hydromechanic Dispersion Tensor for Isotropic Media +############################################################# + + Implements the following parameterization function: + + .. math:: + + \left( \mathsf{D}_\text{hm} \right)_{ij} + = (\lambda_l-\lambda_t)\frac{v_i v_j}{\lvert \mathbf{v} \rvert} + + \delta_{ij}\lambda_t \lvert \mathbf{v} \rvert, + + where :math:`\mathbf{v} \, [\text{ms}^{-1}]` is the local fluid velocity + and :math:`\delta_{ij}` is the Kronecker delta. + + * ``Dhm_type: Dhm_iso`` + + Parameters: + * ``lambda_l``: Longitudinal dispersivity :math:`\lambda_l \, [\text{m}^2\text{s}^{-1}]` + * ``lambda_t``: Transverse dispersivity :math:`\lambda_t \, [\text{m}^2\text{s}^{-1}]` + + YAML template: + + .. code-block:: yaml + + transport: + type: Dhm_iso + + parameters: + mol_diff: + char_length: + lambda_l: + lambda_t: + # parameters ... + + Scalings -^^^^^^^^ +~~~~~~~~ Every ``scaling_section`` has the following layout: .. code-block:: yaml @@ -178,35 +518,36 @@ coordinates in the respective spatial dimensions. omitted. Only omitting the respective *values* of both keys will lead to a parser error. -.. object:: Miller scaling +Miller Scaling +"""""""""""""" - Assumes geometric similarity between scaled regions. Applies the following - scaling: + Assumes geometric similarity between scaled regions. Applies the following + scaling: - .. math :: + .. math :: - \Theta &= \Theta (h_m \cdot \xi_M)\\ - K &= K (\Theta) \cdot \xi_M^2 + \Theta &= \Theta (h_m \cdot \xi_M)\\ + K &= K (\Theta) \cdot \xi_M^2 - * ``type: Miller`` + * ``type: Miller`` - Scaling sections: - * ``scale_miller``: Scaling factor :math:`\xi_M` to be applied onto - matric head and conductivity simultaneously. + Scaling sections: + * ``scale_miller``: Scaling factor :math:`\xi_M` to be applied onto + matric head and conductivity simultaneously. - YAML template: + YAML template: - .. code-block:: yaml + .. code-block:: yaml - type: Miller - data: - scale_miller: - # ... + type: Miller + data: + scale_miller: + # ... .. _man-interpolators: Interpolators -^^^^^^^^^^^^^ +------------- .. object:: Nearest neighbor interpolator diff --git a/doc/python/utilities.rst b/doc/python/utilities.rst index 4e75c9c15548be3024a647d87d64a8e01795fb70..d251d143d52a2de25cbe414bea9036ccb85e8e6b 100644 --- a/doc/python/utilities.rst +++ b/doc/python/utilities.rst @@ -53,12 +53,17 @@ Generic VTK Reader ------------------ This is a flexible VTK XML file reader which can evaluate all data sets written -by DORiE. +by DORiE. It is derived from the abstract base class ``Mapping`` and thus +supports accessing the data arrays stored inside like a Python ``dict``. Use +the data array names as keys. .. autoclass:: dorie.utilities.vtktools.vtkreader.VTKReader - :members: - .. automethod:: __init__ +The ``VTKReader`` returns instances of ``VTKDataArray``, which is a wrapper +around the underlying VTK data array. + +.. autoclass:: dorie.utilities.vtktools.vtkreader.VTKDataArray + :members: dorie.utilities.check_path -------------------------- diff --git a/docker/dune-env.dockerfile b/docker/dune-env.dockerfile index 58c09e06357a71363d54e73952f990e518b6de50..b64c690b80d30ff0327ff5259515b2bf93431814 100644 --- a/docker/dune-env.dockerfile +++ b/docker/dune-env.dockerfile @@ -43,7 +43,7 @@ ENV LANG en_US.utf8 WORKDIR /opt/dune RUN git clone https://gitlab.dune-project.org/staging/dune-uggrid.git -b releases/2.6 \ - && git clone https://gitlab.dune-project.org/oklein/dune-randomfield.git -b master \ + && git clone https://gitlab.dune-project.org/oklein/dune-randomfield.git -b releases/2.6 \ && git clone https://gitlab.dune-project.org/core/dune-common.git -b releases/2.6 \ && git clone https://gitlab.dune-project.org/core/dune-geometry.git -b releases/2.6 \ && git clone https://gitlab.dune-project.org/core/dune-grid.git -b releases/2.6 \ diff --git a/dune/dorie-rfg/CMakeLists.txt b/dune/dorie-rfg/CMakeLists.txt index da3308f3fc79159472b9b86f213e4623fca77ec4..800172eb2e4e1f8c5058398c0a1eccb9c8979ee3 100644 --- a/dune/dorie-rfg/CMakeLists.txt +++ b/dune/dorie-rfg/CMakeLists.txt @@ -1,2 +1,2 @@ -add_executable("dorie-rfg" dorie-rfg.cc) -dune_target_link_libraries(dorie-rfg ${DUNE_LIBS}) \ No newline at end of file +add_executable(dorie-rfg dorie-rfg.cc) +dune_target_link_libraries(dorie-rfg ${DUNE_LIBS}) diff --git a/dune/dorie-rfg/dorie-rfg.cc b/dune/dorie-rfg/dorie-rfg.cc index 5c0d701d40f227db3fd9ea2abdf768f410d2b6f9..5c26b3e32566767cb92a6f3b393b2c918383100f 100644 --- a/dune/dorie-rfg/dorie-rfg.cc +++ b/dune/dorie-rfg/dorie-rfg.cc @@ -1,122 +1,101 @@ -// -*- tab-width: 4; indent-tabs-mode: nil -*- -/** \file - -*/ #ifdef HAVE_CONFIG_H -#include "config.h" + #include "config.h" #endif -// common includes -#include -#include -#include -#include -#include -#include -#include #include #include -#include -#include -#include +#include +#include -// DUNE includes -// Do not treat DUNE warnings as errors -#pragma GCC diagnostic push -#pragma GCC diagnostic warning "-Wall" -#include +#include #include -#include -#include -#include -#include -#pragma GCC diagnostic pop +#include +#include -// dorie-rfg includes #include -/// Traits for the Random Field -template +/// Set up dummy traits required by RandomField (usually supplied by grid) +template struct GridTraits { - enum {dim = dimension}; - using RangeField = double; - using Scalar = Dune::FieldVector; - using DomainField = double; - using Domain = Dune::FieldVector; + enum {dim = dimension}; + using RangeField = double; + using Scalar = Dune::FieldVector; + using DomainField = double; + using Domain = Dune::FieldVector; }; int main(int argc, char** argv) { - try{ - //Initialize Mpi - Dune::MPIHelper::instance(argc, argv); - - if (argc!=2) - DUNE_THROW(Dune::IOError,"No parameter file specified!"); - const std::string inifilename = argv[1]; - - #if !(HDF5_PARALLEL) - if (helper.size() > 1) - DUNE_THROW(Dune::Exception,"Parallel HDF5 library is needed to run in parallel"); - #endif - - // Read ini file - Dune::ParameterTree inifile; - Dune::ParameterTreeParser ptreeparser; - ptreeparser.readINITree(inifilename, inifile); - const unsigned int dim = inifile.get("grid.dimensions"); - - // Attempt to create output directory - const std::string outputPath - = inifile.get("general.tempDir"); - mkdir(outputPath.c_str(), S_IRWXU | S_IRWXG | S_IROTH | S_IXOTH); - int result = access(outputPath.c_str(), W_OK); - if (result != 0) - DUNE_THROW(Dune::IOError,"Output folder " << outputPath << " not writable"); - - const std::string outputFile = outputPath + "/field"; - - // standard values - inifile["stochastic.anisotropy"] = "axiparallel"; - - // extract seed - const unsigned int seed = inifile.get("stochastic.seed"); - - // Create RFG objects - switch(dim){ - case 2: - { - using Traits = GridTraits<2>; - Dune::RandomField::RandomField field(inifile); - field.generate(seed); - field.writeToFile(outputFile); - } - break; - case 3: - { - using Traits = GridTraits<3>; - Dune::RandomField::RandomField field(inifile); - field.generate(seed); - field.writeToFile(outputFile); - } - break; - default: - DUNE_THROW(Dune::NotImplemented,"Only 2 and 3-dimensional fields are supported"); - } - - return 0; - } - - catch (Dune::Exception &e) - { - std::cerr << "Dune reported error: " << e << std::endl; - return 1; - } - catch (...) - { - std::cerr << "Unknown exception thrown!" << std::endl; - throw; - return 1; - } + try{ + //Initialize Mpi + Dune::MPIHelper::instance(argc, argv); + + if (argc!=2) + DUNE_THROW(Dune::IOError,"No parameter file specified!"); + const std::string inifilename = argv[1]; + + // Read ini file + Dune::ParameterTree inifile; + Dune::ParameterTreeParser ptreeparser; + ptreeparser.readINITree(inifilename, inifile); + const unsigned int dim = inifile.get("grid.dimensions"); + + // Attempt to create output directory + const std::string outputPath + = inifile.get("general.tempDir"); + const int status = mkdir(outputPath.c_str(), + S_IRWXU | S_IRWXG | S_IROTH | S_IXOTH); + // allow failure because directory exists + if (status != 0 && errno != EEXIST) + DUNE_THROW(Dune::IOError, + "Output folder " << outputPath << " not writable. " + "Error by system: " << std::strerror(errno)); + + // fix output filename to ensure the file is found by Python frontend + const std::string outputFile = outputPath + "/field"; + + // standard values + inifile["stochastic.anisotropy"] = "axiparallel"; + + // extract seed + const unsigned int seed = inifile.get("stochastic.seed"); + + // Create RFG objects + if (dim == 2) { + using Traits = GridTraits<2>; + Dune::RandomField::RandomField field(inifile); + field.generate(seed); + field.writeToFile(outputFile); + } + else if (dim == 3) { + using Traits = GridTraits<3>; + Dune::RandomField::RandomField field(inifile); + field.generate(seed); + field.writeToFile(outputFile); + } + else { + DUNE_THROW(Dune::NotImplemented, + "Only 2 and 3-dimensional fields are supported"); + } + + return 0; + } + + catch (Dune::Exception &e) + { + std::cerr << "Dune reported error: " << e << std::endl; + return 1; + } + catch (std::exception& e) + { + std::cerr << "Exception occurred: " << e.what() << std::endl; + return 1; + } + catch (...) + { + std::cerr << "Unknown exception thrown!" << std::endl; + throw; + return 1; + } } diff --git a/dune/dorie/CMakeLists.txt b/dune/dorie/CMakeLists.txt index aa8e18c2961036c3c8df6b358c21209d7fbbfc00..f840084fb7c1b64a84b896e1ee4e82f9cc7fc405 100644 --- a/dune/dorie/CMakeLists.txt +++ b/dune/dorie/CMakeLists.txt @@ -4,16 +4,19 @@ if(dune-testtools_FOUND) add_subdirectory("test") endif() -add_executable("richards" richards.cc) -dune_target_link_libraries(richards dorie-richards ${DUNE_LIBS} spdlog) +add_executable(richards richards.cc) +dune_target_link_libraries(richards ${DUNE_LIBS}) +target_link_libraries(richards + dorie-richards spdlog muparser::muparser hdf5 yaml-cpp) -add_executable("transport" transport.cc) -dune_target_link_libraries(transport dorie-richards ${DUNE_LIBS} spdlog) -dune_target_link_libraries(transport dorie-transport ${DUNE_LIBS} spdlog) +add_executable(transport transport.cc) +dune_target_link_libraries(transport ${DUNE_LIBS}) +target_link_libraries(transport + dorie-richards dorie-transport spdlog muparser::muparser hdf5 yaml-cpp) add_custom_target("dorie" DEPENDS richards transport) # enable setting operator scheme from config file if(EXPERIMENTAL_DG_FEATURES) target_compile_definitions("dorie" PUBLIC -DEXPERIMENTAL_DG_FEATURES) -endif() \ No newline at end of file +endif() diff --git a/dune/dorie/common/boundary_condition.hh b/dune/dorie/common/boundary_condition.hh index 25977206dd5a2e8f921a6bb3e430e704ecfd27a5..44cf12efc63e63e634356ab3cbd32a6bf2c61a29 100644 --- a/dune/dorie/common/boundary_condition.hh +++ b/dune/dorie/common/boundary_condition.hh @@ -816,6 +816,8 @@ namespace Dorie{ std::cout << "outflow : "; else if(bcData[i][j].type[k]==BoundaryCondition::Dirichlet) std::cout << "dirichlet : "; + else if(bcData[i][j].type[k]==BoundaryCondition::Outflow) + std::cout << "outflow : "; } index++; } diff --git a/dune/dorie/common/grid_function_writer.hh b/dune/dorie/common/grid_function_writer.hh index 0ef0103b0e4fbb6d8e15f6d08200bea783ea95c7..3573db57b7e2b7caf4fc72f8f68703df3406bf29 100644 --- a/dune/dorie/common/grid_function_writer.hh +++ b/dune/dorie/common/grid_function_writer.hh @@ -16,6 +16,72 @@ namespace Dune{ namespace Dorie{ +/** + * @brief Converts PDELab grid function to match a dune-function. + * @note This adapter differs from the one provided in PDELab in that this + * one has the current interface of dune-grid to write vtk files. + * The main difference using this adapter is that it allows to write + * more than 3 components. + * + * @tparam GF PDELab grid function type. + */ +template +class VTKGridFunctionAdapter +{ + using GridFunction = typename std::decay::type; + using GridView = typename GF::Traits::GridViewType; + using Entity = typename GridView::template Codim<0>::Entity; + using Range = typename GF::Traits::RangeType; + using Domain = typename GF::Traits::DomainType; +public: + + /** + * @brief Constructs the object + * + * @param[in] gf pointer to a constant grid function + */ + VTKGridFunctionAdapter(std::shared_ptr gf) + : _gf(gf) + , _entity(nullptr) + {} + + /** + * @brief Binds an entity to the grid function + * + * @param[in] e Entity + */ + void bind(const Entity& e) + { + _entity = &e; + } + + /** + * @brief Unbinds the previously binded entity + */ + void unbind() + { + _entity = nullptr; + } + + /** + * @brief Evaluates the grid function at a given coordinate + * + * @param[in] x Local coordinate with respect to the last binded entity + */ + Range operator()(const Domain& x) const + { + Range y; + assert(_entity); + _gf->evaluate(*_entity,x,y); + return y; + } + +private: + const Entity* _entity; + std::shared_ptr _gf; +}; + + /*-------------------------------------------------------------------------*//** * @brief Class for grid fucntion vtk sequence writer. This class works * exacly the same way as a VTKSequence writer but it receives @@ -39,13 +105,20 @@ public: * @tparam GF Type of the grid function */ template - void addVertexData(const std::shared_ptr &gf_ptr, std::string name) + void addVertexData(std::shared_ptr gf_ptr, std::string name) { static_assert(std::is_same::value, "GridFunctionVTKSequenceWriter only works with only one GridView type."); - auto vtk_gf_ptr = std::make_shared >(gf_ptr, name); - Dune::VTKSequenceWriter::addVertexData(vtk_gf_ptr); + Dune::Dorie::VTKGridFunctionAdapter vtk_gf(gf_ptr); + + const int vtk_ncomps = GF::Traits::dimRange; + const int dim = GF::Traits::dimDomain; + auto vtk_type = (vtk_ncomps == dim) ? VTK::FieldInfo::Type::vector + : VTK::FieldInfo::Type::scalar; + auto vtk_info = VTK::FieldInfo(name,vtk_type,vtk_ncomps); + + this->vtkWriter()->addVertexData(vtk_gf,vtk_info); } /*-----------------------------------------------------------------------*//** @@ -62,14 +135,21 @@ public: static_assert(std::is_same::value, "GridFunctionVTKSequenceWriter only works with only one GridView type."); - auto vtk_gf_ptr = std::make_shared >(gf_ptr, name); - Dune::VTKSequenceWriter::addCellData(vtk_gf_ptr); + Dune::Dorie::VTKGridFunctionAdapter vtk_gf(gf_ptr); + + const int vtk_ncomps = GF::Traits::dimRange; + const int dim = GF::Traits::dimDomain; + auto vtk_type = (vtk_ncomps == dim) ? VTK::FieldInfo::Type::vector + : VTK::FieldInfo::Type::scalar; + auto vtk_info = VTK::FieldInfo(name,vtk_type,vtk_ncomps); + + this->vtkWriter()->addCellData(vtk_gf,vtk_info); } /*-----------------------------------------------------------------------*//** - * @brief Clear all the attached VTKFunctions. (Necessary for each - * iteration if VTKFunctions are managed internaly with pointers - * instead of references) + * @brief Clear all the attached VTKFunctions. + * @details This is necessary for each iteration if VTKFunctions are + * managed internaly with pointers instead of references */ void clear() { diff --git a/dune/dorie/common/initial_condition/analytic.hh b/dune/dorie/common/initial_condition/analytic.hh index c14578e98069d49e964c3ce2131da4ce361f1cfc..7df8f3bf9ac80ca90f574688aaa09b132808d3ad 100644 --- a/dune/dorie/common/initial_condition/analytic.hh +++ b/dune/dorie/common/initial_condition/analytic.hh @@ -20,7 +20,7 @@ namespace Dorie { * @brief Class for analytic initial conditions. * * @ingroup InitialConditions - * @author Santiago Ospina + * @author Santiago Ospina De Los Ríos * @date 2019 * * @tparam T BaseTraits diff --git a/dune/dorie/common/initial_condition/factory.hh b/dune/dorie/common/initial_condition/factory.hh index 5d3ed2e3fdaf433b9a76d7e65dd06a6a1bdd5c00..edc10765156c688c107f662fbd4f9322c827b0ad 100644 --- a/dune/dorie/common/initial_condition/factory.hh +++ b/dune/dorie/common/initial_condition/factory.hh @@ -17,7 +17,7 @@ namespace Dorie{ * @brief Abstract Base Class for initial conditions. * * @ingroup InitialConditions - * @author Santiago Ospina + * @author Santiago Ospina De Los Ríos * @date 2019 * * @tparam T BaseTraits diff --git a/dune/dorie/common/initial_condition/h5.hh b/dune/dorie/common/initial_condition/h5.hh index cbf811352cb6784762d33068b61c620d5fc1a18c..10a15942e8370388ddeb75d63575e40ac26c3dbe 100644 --- a/dune/dorie/common/initial_condition/h5.hh +++ b/dune/dorie/common/initial_condition/h5.hh @@ -18,7 +18,7 @@ namespace Dorie{ * @brief Class for initial conditions from hdf5 data files. * * @ingroup InitialConditions - * @author Santiago Ospina + * @author Santiago Ospina De Los Ríos * @date 2019 * * @tparam T BaseTraits diff --git a/dune/dorie/common/initial_condition/initial_condition.hh b/dune/dorie/common/initial_condition/initial_condition.hh index bb46a9f6683488eadb42a7dd08917c25b3c64027..29a3957c72c6bb6f3aa2a3cee0bab3cb72596ccc 100644 --- a/dune/dorie/common/initial_condition/initial_condition.hh +++ b/dune/dorie/common/initial_condition/initial_condition.hh @@ -13,7 +13,7 @@ namespace Dorie{ * @brief Abstract Base Class for initial conditions. * * @ingroup InitialConditions - * @author Santiago Ospina + * @author Santiago Ospina De Los Ríos * @date 2019 * * @tparam T BaseTraits diff --git a/dune/dorie/common/parameterization_factory.hh b/dune/dorie/common/parameterization_factory.hh new file mode 100644 index 0000000000000000000000000000000000000000..3017e9c2c50a66542141506922acc8d9c988c0ac --- /dev/null +++ b/dune/dorie/common/parameterization_factory.hh @@ -0,0 +1,120 @@ +#ifndef DUNE_DORIE_PARAMETERIZATION_READER_HH +#define DUNE_DORIE_PARAMETERIZATION_READER_HH + +#include +#include +#include +#include + +#include + +#include + +#include +#include + +namespace Dune{ +namespace Dorie{ + +/*-------------------------------------------------------------------------*//** + * @brief Interface and default reader for parameterization factories + * + * @tparam Param The parameter type to build + */ +template +struct ParameterizationFactory +{ + virtual ~ParameterizationFactory() = default; + + /** + * @brief Parameterization selector + * @details Given a type_node in yaml, it selects a -polymorphic- + * parameter. + * + * @param[in] type_node YAML node describing the type of parameterization to + * choose from. + * @param[in] name Name that describes the parametrized material. + * + * @return Shared pointer to the selected parameter. + */ + virtual std::shared_ptr selector( + const YAML::Node& type_node, + const std::string& name) const = 0; + + + /** + * @brief Parameterization reader + * @details It reads parameterizations and its parameters. For each volume, + * it assings a index to a parameterization (choosen by the + * parameterization selector) and later filled with the paramaters + * in the yaml file. + * + * @param[in] param_file YAML node with information about type and + * parameters + * @param[in] param_name Name of the parameterization. + * @param[in] log Logger + * + * @return (Unsorted) Map with parameterizations for each material. + */ + std::unordered_map> + reader( const YAML::Node& param_file, + const std::string& param_name, + const std::shared_ptr log) const + { + // mapping: index on grid to parameterization object (will be returned) + std::unordered_map> ret; + + // set to check for duplicate indices + std::set param_index_set; + + // iterate over all volume nodes + for (auto&& node : param_file["volumes"]) + { + const auto volume = node.first.as(); + const YAML::Node param_info = node.second; + + // read relevant data + const auto index = param_info["index"].as(); + if (param_index_set.count(index) != 0) { + log->error("Index '{}' in volume '{}' of parameterization '{}'" + "already exists.", + index, volume, param_name); + DUNE_THROW(IOError, "Parameter file indices must be unique"); + } + param_index_set.emplace(index); + + // get parameterization object + const YAML::Node& type_node = param_info[param_name]; + auto parameterization = this->selector(type_node, volume); + auto parameter_values = parameterization->parameters(); + + // fill parameterization with parameter entries + log->trace(" Reading parameters. Volume: {}, Index: {}", + volume, index); + + auto yaml_parameters = yaml_sequence_to_map(type_node["parameters"]); + + for (auto&& [name, value] : parameter_values) + { + try { + value = yaml_parameters[name].template as(); + } + catch (...) { + log->error(" Parameter '{}' is missing in '{}' parameterization," + " volume {}, index {}.", + name, param_name, volume, index); + DUNE_THROW(IOError, "Missing parameters!"); + } + } + + ret.emplace(index, parameterization); + } + + return ret; + } +}; + +} +} + +#endif // DUNE_DORIE_PARAMETERIZATION_READER_HH \ No newline at end of file diff --git a/dune/dorie/common/setup.hh b/dune/dorie/common/setup.hh index 5e68fb06c4bb3dcac84fbe2ba6259949fa3bb080..4f5ed21ed5b8c8bece886bec1876acf8117540e3 100644 --- a/dune/dorie/common/setup.hh +++ b/dune/dorie/common/setup.hh @@ -63,7 +63,7 @@ namespace Setup const std::string inifilename = argv[1]; // Read ini file - log->info("Reading parameter file: {}", inifilename); + log->info("Reading configuration file: {}", inifilename); Dune::ParameterTree inifile; Dune::ParameterTreeParser ptreeparser; ptreeparser.readINITree(inifilename, inifile); diff --git a/dune/dorie/common/typedefs.hh b/dune/dorie/common/typedefs.hh index 9f57c9ffdd09f6d4932aff9c477c2ad05a06f382..57407b17944ef20894113e5d0771765d39d5057c 100644 --- a/dune/dorie/common/typedefs.hh +++ b/dune/dorie/common/typedefs.hh @@ -13,7 +13,7 @@ namespace Dune{ enum Type { Neumann, //!< Fixed flux at boundary Dirichlet, //!< Fixed matric head at boundary - Outflow, //!< Fixed flux at boundary + Outflow, //!< Fixed outflow boundary condition Other //!< UNUSED }; diff --git a/dune/dorie/common/util.hh b/dune/dorie/common/util.hh index bbce2333c7e77dc77bff5f9a9e9e72bbe2cff728..50acb29f91e5ea53f9f92bee6a544e882d85e234 100644 --- a/dune/dorie/common/util.hh +++ b/dune/dorie/common/util.hh @@ -195,16 +195,6 @@ public: } }; -/// A wrapper for begin and end time stamps of a time interval -/** \tparam TF Arithmetic data type for time stamps - */ -template -struct TimeInterval { - TF begin; - TF end; -}; - - /// Traits struct defining basic types based on grid an geometry template struct BaseTraits @@ -214,7 +204,6 @@ struct BaseTraits using TF = double; using TimeField = TF; - using TimeInterval = Dune::Dorie::TimeInterval; using RF = double; using RangeField = RF; diff --git a/dune/dorie/common/utility.hh b/dune/dorie/common/utility.hh index fbd2c4c159eac7cc3ae9cb2c7a084114cf9ed67b..7ca6275343a87081f4b4191904f163fd2ffbcc0e 100644 --- a/dune/dorie/common/utility.hh +++ b/dune/dorie/common/utility.hh @@ -8,6 +8,9 @@ #include #include +#include + +#include #include #include @@ -113,6 +116,61 @@ inline bool is_none (const std::string& in) return false; } +/*-------------------------------------------------------------------------*//** + * @brief Transform a Yaml node with sequences to map + * @details This function check all the nodes within 'in' and transform every + * sequence into a map. In particular it transform vectors and + * tensors into single key-value map pair by adding a suffix + * indicating the component. + * + * @param[in] in Map Yaml node. + * + * @return Map Yaml node with only key-double pairs. + */ +inline YAML::Node yaml_sequence_to_map(const YAML::Node& in) +{ + YAML::Node out; + for (auto&& node : in) + { + const auto key = node.first.as(); + const YAML::Node& value = node.second; + + if (value.IsScalar()) { + out[key] = value; + } else if (value.IsSequence()) { + if (value.size() == 1) { // scalar + out[key] = value[0].as(); + } else if (value.size() == 2) { // assume vector in 2D + out[key+"_x"] = value[0].as(); + out[key+"_y"] = value[1].as(); + } else if (value.size() == 3) { // assume vector in 3D + out[key+"_x"] = value[0].as(); + out[key+"_y"] = value[1].as(); + out[key+"_z"] = value[2].as(); + } else if (value.size() == 4) { // assume tensor in 2D + out[key+"_xx"] = value[0].as(); + out[key+"_xy"] = value[1].as(); + out[key+"_yx"] = value[2].as(); + out[key+"_yy"] = value[3].as(); + } else if (value.size() == 9) { // assume tensor in 3D + out[key+"_xx"] = value[0].as(); + out[key+"_xy"] = value[1].as(); + out[key+"_xz"] = value[2].as(); + out[key+"_yx"] = value[3].as(); + out[key+"_yy"] = value[4].as(); + out[key+"_yz"] = value[5].as(); + out[key+"_zx"] = value[6].as(); + out[key+"_zy"] = value[7].as(); + out[key+"_zz"] = value[8].as(); + } + } else { + DUNE_THROW(Dune::IOError, + "Yaml sequence '" << key << "' is not convertible to map pair!"); + } + } + return out; +} + } // namespace Dorie } // namespace Dune diff --git a/dune/dorie/model/richards/adapters/conductivity.hh b/dune/dorie/model/richards/adapters/conductivity.hh index 1f0e06ee891d6e7f2cbdae55a30ac12a7d4d40e9..b0a57d65c2daa19d90dcb42cb0867e44f586678b 100644 --- a/dune/dorie/model/richards/adapters/conductivity.hh +++ b/dune/dorie/model/richards/adapters/conductivity.hh @@ -9,7 +9,8 @@ namespace Dune{ /*---------------------------------------------------------------------*//** * @brief Converts a parametrization P that contains conductivity * information into a grid function (in PDELab sense) - * + * @ingroup RichardsParam + * * @tparam T BasicTraits with information about the function domain * @tparam P Parametrization class with saturated conductivity * information diff --git a/dune/dorie/model/richards/adapters/saturation.hh b/dune/dorie/model/richards/adapters/saturation.hh index f4bc8080f4782def4dab3075c8b775d7c60b3c73..e14805c15f084de81862b2d26185ae1344c4229f 100755 --- a/dune/dorie/model/richards/adapters/saturation.hh +++ b/dune/dorie/model/richards/adapters/saturation.hh @@ -10,7 +10,8 @@ namespace Dune{ * @brief Converts a parametrization P that contains saturation * information and a matric head grid function GF into a grid * function (in PDELab sense) of saturation. - * + * @ingroup RichardsParam + * * @tparam T BasicTraits with information about the function domain * @tparam P Parametrization class with saturation information * @tparam GF Grid function of the matric head diff --git a/dune/dorie/model/richards/adapters/water_content.hh b/dune/dorie/model/richards/adapters/water_content.hh index c2977b71cdd890530b26e701d94d6d6d5e9856e8..4032f7af9a3ca66af29b79df323950eb25f0475a 100644 --- a/dune/dorie/model/richards/adapters/water_content.hh +++ b/dune/dorie/model/richards/adapters/water_content.hh @@ -10,7 +10,8 @@ namespace Dune{ * @brief Converts a parametrization P that contains conductivity * information and a matric head grid function GF into a grid * function (in PDELab sense) of water content. - * + * @ingroup RichardsParam + * * @tparam T BasicTraits with information about the function domain * @tparam P Parametrization class with conductivity information * @tparam GF Grid function of the matric head diff --git a/dune/dorie/model/richards/adapters/water_flux.hh b/dune/dorie/model/richards/adapters/water_flux.hh index 5459d33d1f4f113c3960b08e4e4a1da00bd34768..13a12e2f215d3666bc41f9e7c7ac0b46d254997e 100644 --- a/dune/dorie/model/richards/adapters/water_flux.hh +++ b/dune/dorie/model/richards/adapters/water_flux.hh @@ -11,7 +11,8 @@ namespace Dune{ * @brief Converts an grid function space GFS and some particular * coefficients for such space X to match the water flux * equation with a GridFunctionInterface. - * + * @ingroup RichardsParam + * * @tparam GFS Grid Function Space * @tparam X Coefficients for the GFS * @tparam P Parameters diff --git a/dune/dorie/model/richards/flow_parameters.hh b/dune/dorie/model/richards/flow_parameters.hh index ca3e134992b2a58c3b4a1fc9dc25af4fe5cbc7f8..3aa9ef401f9c6d7a8e2c4343a5b3b11beeca73da 100644 --- a/dune/dorie/model/richards/flow_parameters.hh +++ b/dune/dorie/model/richards/flow_parameters.hh @@ -14,8 +14,7 @@ #include #include -#include -#include +#include #include namespace Dune { @@ -64,21 +63,23 @@ private: /// Index type used by the element mapper using Index = typename Mapper::Index; /// Base class of the parameterization - using Parameterization = Dorie::RichardsParameterization; + using ParameterizationType = Dorie::Parameterization::Richards; + /// Parameterization factory + using ParameterizationFactory = Dorie::Parameterization::RichardsFactory; /// Struct for storing scaling factors - using Scaling = Dorie::ScalingFactors; + using Scaling = Dorie::Parameterization::ScalingFactors; /// Base class for scaling adapters - using ScaleAdapter = Dorie::ScalingAdapter; + using ScaleAdapter = Dorie::Parameterization::ScalingAdapter; /// Storage for parameters and skaling factors using ParameterStorage = std::unordered_map, Scaling> + std::tuple, Scaling> >; /// Need this gruesome typedef because 'map::value_type' has const key - using Cache = std::tuple, Scaling>; + using Cache = std::tuple, Scaling>; - using ConstCache = std::tuple, Scaling>; + using ConstCache = std::tuple, Scaling>; /// Configuration file tree const Dune::ParameterTree& _config; /// Grid view of the coarsest grid configuration (level 0) @@ -99,7 +100,7 @@ private: /// Check if an entity has been cached void verify_cache () const { - if (not std::get>(_cache)) { + if (not std::get>(_cache)) { _log->error("Parameterization cache is empty. Call 'bind' before " "accessing the cache or the parameterization"); DUNE_THROW(Dune::InvalidStateException, @@ -126,7 +127,7 @@ public: for(const auto& [index, tuple] : flow_param._param) { const auto& [p, sk] = tuple; // make a hard copy of parameterization - std::shared_ptr> _p = p->clone(); + std::shared_ptr> _p = p->clone(); this->_param.emplace(index,std::make_tuple(_p,sk)); } } @@ -197,9 +198,9 @@ public: { verify_cache(); const auto& par = - std::get>(_cache); + std::get>(_cache); const auto cond_f = par->conductivity_f(); - using Saturation = typename Parameterization::Saturation; + using Saturation = typename ParameterizationType::Saturation; const auto& xi_cond = std::get(_cache).scale_cond; @@ -217,9 +218,9 @@ public: { verify_cache(); const auto& par = - std::get>(_cache); + std::get>(_cache); const auto sat_f = par->saturation_f(); - using MatricHead = typename Parameterization::MatricHead; + using MatricHead = typename ParameterizationType::MatricHead; const auto& scale_head = std::get(_cache).scale_head; @@ -237,13 +238,13 @@ public: { verify_cache(); const auto& par = - std::get>(_cache); + std::get>(_cache); // get water content function and apply the scaling const auto& scale_por = std::get(_cache).scale_por; const auto wc_f = par->water_content_f(scale_por); - using Saturation = typename Parameterization::Saturation; + using Saturation = typename ParameterizationType::Saturation; return [wc_f](const RF saturation) { return wc_f(Saturation{saturation}).value; }; @@ -258,12 +259,12 @@ private: const auto param_file_name = _config.get("parameters.file"); _log->info("Loading parameter data file: {}", param_file_name); YAML::Node param_file = YAML::LoadFile(param_file_name); - // create the parameterization data - const auto parameterization_map = read_parameter_file(param_file); + ParameterizationFactory factory; + const auto parameterization_map = factory.reader(param_file,"richards",_log); // build global scaling - using SAF = Dorie::ScalingAdapterFactory; + using SAF = Dorie::Parameterization::ScalingAdapterFactory; auto scale_cfg = param_file["scaling"]; auto scaling_adapter = SAF::create(scale_cfg["type"].as(), scale_cfg["data"], @@ -275,95 +276,14 @@ private: // evaluate element index and position const auto index = _mapper.index(elem); const auto pos = elem.geometry().center(); - // read values and insert const auto p = parameterization_map.at(element_index_map.at(index)); const auto scale = scaling_adapter->evaluate(pos); _param.emplace(index,std::make_tuple(p, scale)); } - // check that mapper can map to parameterization assert(_param.size() == _mapper.size()); } - - /// Read the YAML parameter file and insert information into a map. - /** \param param_file The top YAML node to read from (file root) - */ - std::map> read_parameter_file ( - const YAML::Node& param_file - ) - { - // mapping: index on grid to parameterization object (will be returned) - std::map> ret; - - // set to check for duplicate indices - std::set param_index_set; - - // iterate over all volume nodes - for (auto&& node : param_file["volumes"]) - { - const auto name = node.first.as(); - const YAML::Node param_info = node.second; - - // read relevant data - const auto type = param_info["type"].as(); - const auto index = param_info["index"].as(); - if (param_index_set.count(index) != 0) { - _log->error("Index '{}' of parameterization '{}'" - "already exists.", - index, name); - DUNE_THROW(IOError, "Parameter file indices must be unique"); - } - param_index_set.emplace(index); - - // get parameterization object - auto parameterization = parameterization_selector(type, name); - auto parameter_values = parameterization->parameters(); - - // fill parameterization with parameter entries - _log->trace(" Reading parameters. Volume: {}, Type: {}, Index: {}", - name, type, index); - const YAML::Node yaml_parameters = param_info["parameters"]; - for (auto&& yaml_param : yaml_parameters) - { - const auto p_name = yaml_param.first.as(); - const auto p_value = yaml_param.second.as(); - - // try to insert the value into the parameter map - try { - parameter_values.at(p_name) = p_value; - } - catch (std::out_of_range& e) { - _log->error(" Parameter '{}' not part of " - "parameterization '{}'", - p_name, type); - DUNE_THROW(IOError, "Error inserting parameter values"); - } - } - - ret.emplace(index, parameterization); - } - - return ret; - } - - /// Return a default-initialized parameterization object - /** \param type The type indicating the parameterization implementation - * \param name The name of the parameterization object (layer type) - */ - std::shared_ptr parameterization_selector ( - const std::string type, - const std::string name) const - { - if (type == MualemVanGenuchtenParameterization::type) { - return std::make_shared>(name); - } - else { - _log->error("Parameterization '{}' has unknown type '{}'", - name, type); - DUNE_THROW(IOError, "Unknown parameterization type"); - } - } }; } // namespace Dune diff --git a/dune/dorie/model/richards/modules.dox b/dune/dorie/model/richards/modules.dox index 52dd574bdd80c38f11ff004b511c9ec6d4ec379c..19f2a67ddc69face5003abd94bac4207657ee5d6 100644 --- a/dune/dorie/model/richards/modules.dox +++ b/dune/dorie/model/richards/modules.dox @@ -51,7 +51,7 @@ in the @ref LocalOperators. ## Scaling The parameterizations implemented with the -Dune::Dorie::RichardsParameterization interface are envisaged to be constant +Dune::Dorie::Parameterization::Richards interface are envisaged to be constant throughout a single medium layer. However, one typically finds small-scale variations at the scale of interest. To account for these heterogeneities, several scaling schemes have been developed. @@ -59,7 +59,8 @@ several scaling schemes have been developed. We implement scaling by storing a set of scaling factors for every grid cell on the coarsest level. Unlike the parameterization objects, these do not map to a certain medium but only the cell itself (and its child cells). The scaling -factors are stored inside instances of Dune::Dorie::ScalingFactors. +factors are stored inside instances of +Dune::Dorie::Parameterization::ScalingFactors. For evaluating scaling factor data, the Dune::Dorie::ScalingAdapter interface is used, which internally uses the Interpolator interface. These structures are diff --git a/dune/dorie/model/richards/parameterization/factory.hh b/dune/dorie/model/richards/parameterization/factory.hh new file mode 100644 index 0000000000000000000000000000000000000000..1cb27d42a8fa44295bef41dcb726d9129f577077 --- /dev/null +++ b/dune/dorie/model/richards/parameterization/factory.hh @@ -0,0 +1,50 @@ +#ifndef DUNE_DORIE_PARAM_RICHARDS_FACTORY_HH +#define DUNE_DORIE_PARAM_RICHARDS_FACTORY_HH + +#include + +#include +#include +#include +#include + +namespace Dune { +namespace Dorie { +namespace Parameterization { + +/*-------------------------------------------------------------------------*//** + * @brief Factory of Transport parameters + * + * @tparam Traits The base traits + */ +template +struct RichardsFactory + : public ParameterizationFactory> +{ + /// Return a default-initialized parameterization object + /** \param type The type indicating the parameterization implementation + * \param name The name of the parameterization object (layer type) + */ + std::shared_ptr> + selector( + const YAML::Node& type_node, + const std::string& name) const override + { + auto log = Dorie::get_logger(log_richards); + const auto type = type_node["type"].as(); + if (type == Parameterization::MualemVanGenuchten::type) { + return std::make_shared>(name); + } + else { + log->error("Richards parameterization '{}' has unknown type '{}'", + name, type); + DUNE_THROW(IOError, "Unknown Richards parameterization type"); + } + } +}; + +} // namespace Parameterization +} // namespace Dorie +} // namespace Dune + +#endif // DUNE_DORIE_PARAM_RICHARDS_FACTORY_HH \ No newline at end of file diff --git a/dune/dorie/model/richards/parameterization/interface.hh b/dune/dorie/model/richards/parameterization/interface.hh index 1e4af1cc9bc3b9967fe25bbe5118a87453032b2f..ed8b8d9271dfc766a223d13fe5a546540f01914a 100644 --- a/dune/dorie/model/richards/parameterization/interface.hh +++ b/dune/dorie/model/richards/parameterization/interface.hh @@ -1,5 +1,5 @@ -#ifndef DUNE_DORIE_PARAM_RICHARDS_HH -#define DUNE_DORIE_PARAM_RICHARDS_HH +#ifndef DUNE_DORIE_PARAM_RICHARDS_INTERFACE_HH +#define DUNE_DORIE_PARAM_RICHARDS_INTERFACE_HH #include #include @@ -7,6 +7,7 @@ namespace Dune { namespace Dorie { +namespace Parameterization { /// Interface for parameterizations of Richards equation. /** Defines the basic parameters and the relation between water content @@ -23,7 +24,7 @@ namespace Dorie { * \ingroup RichardsParam */ template -class RichardsParameterization +class Richards { private: using RF = typename Traits::RF; @@ -91,7 +92,7 @@ public: /** \param name The name associated with this soil layer * \param parameters Tuple of parameters to use in this parameterization */ - RichardsParameterization (const std::string name) : + Richards (const std::string name) : _name(name) { } @@ -100,7 +101,7 @@ public: * \param parameters Tuple of parameters to use in this parameterization */ template - RichardsParameterization (const std::string name, + Richards (const std::string name, const std::tuple parameters) : _theta_r(std::get(parameters)), _theta_s(std::get(parameters)), @@ -109,7 +110,7 @@ public: { } /// Default constructor (virtual). - virtual ~RichardsParameterization () = default; + virtual ~Richards () = default; /// Return the name of this parameterization instance. const std::string& get_name() const { return _name; } @@ -142,27 +143,28 @@ public: virtual std::function conductivity_f () const = 0; - /// Return a map referecing all parameters by their names. + /// Return a multimap referecing all parameters by their names. /** \return Map. Key: Name of parameter (string). - * Value: Value of parameter (double&) + * Value(s): Value of parameter (double&) */ - virtual std::map parameters () = 0; + virtual std::multimap parameters () = 0; /// Return a map referecing all parameters by their names. /** \return Map. Key: Name of parameter (string). * Value: Value of parameter (const double&) */ - virtual std::map parameters () const = 0; + virtual std::multimap parameters () const = 0; /// Return a clone of this object /** \return a unique pointer with a copy of this object. */ - virtual std::unique_ptr> clone () const = 0; + virtual std::unique_ptr> clone () const = 0; }; -} // namespace Dune +} // namespace Parameterization } // namespace Dorie +} // namespace Dune -#endif // DUNE_DORIE_PARAM_RICHARDS_HH +#endif // DUNE_DORIE_PARAM_RICHARDS_INTERFACE_HH \ No newline at end of file diff --git a/dune/dorie/model/richards/parameterization/mualem_van_genuchten.hh b/dune/dorie/model/richards/parameterization/mualem_van_genuchten.hh index d08bc3cd7adc579b50c450273c304e5279b7b771..11c77ad529ee1f26061db8dfe8f88fe8476d26df 100644 --- a/dune/dorie/model/richards/parameterization/mualem_van_genuchten.hh +++ b/dune/dorie/model/richards/parameterization/mualem_van_genuchten.hh @@ -9,6 +9,7 @@ namespace Dune { namespace Dorie { +namespace Parameterization { /// Representation of the Mualem–van Genuchten parameterization /** Implements the following relations: @@ -23,12 +24,12 @@ namespace Dorie { * \date 2018 */ template -class MualemVanGenuchtenParameterization : - public RichardsParameterization +class MualemVanGenuchten : + public Richards { private: using RF = typename Traits::RF; - using Base = RichardsParameterization; + using Base = Richards; public: @@ -73,8 +74,8 @@ public: /// Construct with default-initialized parameters /** \param name The name associated with this soil layer */ - MualemVanGenuchtenParameterization (const std::string name) : - RichardsParameterization(name) + MualemVanGenuchten (const std::string name) : + Richards(name) { } /// Construct from a tuple of parameters @@ -82,17 +83,17 @@ public: * \param parameters Tuple of parameters to use in this parameterization */ template - MualemVanGenuchtenParameterization ( + MualemVanGenuchten ( const std::string name, const std::tuple parameters) : - RichardsParameterization(name, parameters), + Richards(name, parameters), _alpha(std::get(parameters)), _tau(std::get(parameters)), _n(std::get(parameters)) { } /// Add default destructor to clarify override - ~MualemVanGenuchtenParameterization () override = default; + ~MualemVanGenuchten () override = default; /// Return a scaled conductivity function /** Saturation -> Conductivity @@ -130,7 +131,7 @@ public: /// Return a map of parameter names and values for manipulation /** \return Map: Parameter name, parameter value in this object */ - std::map parameters () override + std::multimap parameters () override { return { {ResidualWaterContent::name, this->_theta_r.value}, @@ -145,7 +146,7 @@ public: /// Return a map of parameter names and values for manipulation /** \return Map: Parameter name, parameter value in this object */ - std::map parameters () const override + std::multimap parameters () const override { return { {ResidualWaterContent::name, this->_theta_r.value}, @@ -157,14 +158,18 @@ public: }; } - std::unique_ptr> clone () const override + /// Clone the plymorphic class + /** \return unique pointer to the cloned object + */ + std::unique_ptr> clone () const override { - using ThisType = MualemVanGenuchtenParameterization; + using ThisType = MualemVanGenuchten; return std::make_unique(*this); } }; -} // namespace Dune +} // namespace Parameterization } // namespace Dorie +} // namespace Dune #endif // DUNE_DORIE_PARAM_MVG_HH diff --git a/dune/dorie/model/richards/parameterization/scaling.hh b/dune/dorie/model/richards/parameterization/scaling.hh index d2cc4879e6c69116c065c599bda64b7c36d819e2..d65563a5937eec863b24811fde78441b9a81357d 100644 --- a/dune/dorie/model/richards/parameterization/scaling.hh +++ b/dune/dorie/model/richards/parameterization/scaling.hh @@ -16,6 +16,7 @@ namespace Dune { namespace Dorie { +namespace Parameterization { /// Define scaling factors used for any scaling type. /** The variables are initialized to their default values @@ -252,7 +253,7 @@ public: } }; - +} // namespace Parameterization } // namespace Dorie } // namespace Dune diff --git a/dune/dorie/model/richards/richards.cc b/dune/dorie/model/richards/richards.cc index e4730fe03ef7d7646842bdc1a4ea8f37beee3893..13bf00945cf8625de0e18288dfda0c1279aa8061 100644 --- a/dune/dorie/model/richards/richards.cc +++ b/dune/dorie/model/richards/richards.cc @@ -245,6 +245,9 @@ void RichardsSimulation::step() template void RichardsSimulation::mark_grid() { + if (adaptivity_policy() != AdaptivityPolicy::WaterFlux) + DUNE_THROW(Dune::NotImplemented,"Not known adaptivity policy for transport solver"); + adaptivity->mark_grid(*grid, gv, *gfs, *fparam, *fboundary, time_before+dt_before, *u); } @@ -276,34 +279,34 @@ void RichardsSimulation::write_data () const { // create adapters auto head = get_matric_head(); - auto wflux = std::make_shared(gfs, u, fparam); - auto cond = std::make_shared(gv, fparam); - auto wc = std::make_shared(head, gv, fparam); - auto sat = std::make_shared(head, gv, fparam); + auto wflux = std::make_shared(gfs, u, fparam); + auto cond = std::make_shared(gv, fparam); + auto wc = std::make_shared(head, gv, fparam); + auto sat = std::make_shared(head, gv, fparam); if (inifile.get("output.vertexData")) { - vtkwriter->template addVertexData(head,"head"); - vtkwriter->template addVertexData(wflux,"flux"); - vtkwriter->template addVertexData(cond,"K_0"); - vtkwriter->template addVertexData(wc,"theta_w"); - vtkwriter->template addVertexData(sat,"Theta"); + vtkwriter->addVertexData(head,"head"); + vtkwriter->addVertexData(wflux,"flux"); + vtkwriter->addVertexData(cond,"K_0"); + vtkwriter->addVertexData(wc,"theta_w"); + vtkwriter->addVertexData(sat,"Theta"); if constexpr (enable_rt_engine) if (enable_fluxrc) { auto wfluxr = get_water_flux_reconstructed(); auto RT_name = "flux_RT" + std::to_string(flux_order); - vtkwriter->template addVertexData(wfluxr,RT_name); + vtkwriter->addVertexData(wfluxr,RT_name); } } else { - vtkwriter->template addCellData(head,"head"); - vtkwriter->template addCellData(wflux,"flux"); - vtkwriter->template addCellData(cond,"K_0"); - vtkwriter->template addCellData(wc,"theta_w"); - vtkwriter->template addCellData(sat,"Theta"); + vtkwriter->addCellData(head,"head"); + vtkwriter->addCellData(wflux,"flux"); + vtkwriter->addCellData(cond,"K_0"); + vtkwriter->addCellData(wc,"theta_w"); + vtkwriter->addCellData(sat,"Theta"); if constexpr (enable_rt_engine) if (enable_fluxrc) { auto wfluxr = get_water_flux_reconstructed(); auto RT_name = "flux_RT" + std::to_string(flux_order); - vtkwriter->template addCellData(wfluxr,RT_name); + vtkwriter->addCellData(wfluxr,RT_name); } } diff --git a/dune/dorie/model/transport/adapters/effective_hydrodynamic_dispersion.hh b/dune/dorie/model/transport/adapters/effective_hydrodynamic_dispersion.hh new file mode 100644 index 0000000000000000000000000000000000000000..1d408521c8594609ef5af29026928b0db1fdc428 --- /dev/null +++ b/dune/dorie/model/transport/adapters/effective_hydrodynamic_dispersion.hh @@ -0,0 +1,121 @@ +#ifndef DUNE_DORIE_HYDRODYNAMIC_DISPERSION_ADAPTER_HH +#define DUNE_DORIE_HYDRODYNAMIC_DISPERSION_ADAPTER_HH + +#include + +namespace Dune{ +namespace Dorie{ + +/*-------------------------------------------------------------------------*//** + * @brief Converts an interface to match an effective hydrodynamic dispersion. + * @details For a given paramaterization, water flux and water content, this + * class returns the hydrodynamic dispersion tensor (od dimension + * dimRange) in an aliased vector. + * @ingroup TransportParam + * + * @tparam T The base traits + * @tparam Parameter The transport parameter class + * @tparam GFWaterFlux The water flux grid function class + * @tparam GFWaterContent The water content grid function class + * @tparam dimRange The range dimension of the aliased tensor + */ +template +class EffectiveHydrodynamicDispersionAdapter + : public Dune::PDELab::GridFunctionBase< + Dune::PDELab::GridFunctionTraits< + typename T::GridView, + typename T::RangeField, + dimRange*dimRange, + Dune::FieldVector + > + ,EffectiveHydrodynamicDispersionAdapter + > +{ +public: + using Traits = Dune::PDELab::GridFunctionTraits< + typename T::GridView, + typename T::RangeField, + dimRange*dimRange, + Dune::FieldVector + >; + +private: + using RF = typename T::RangeField; + using Domain = typename T::Domain; + +public: + /*-------------------------------------------------------------------*//** + * @brief Constructor for the EffectiveHydrodynamicDispersionAdapter + * + * @param gv GridView + * @param p Parametrization class + * @see RichardsEquationParameter + */ + EffectiveHydrodynamicDispersionAdapter(const typename Traits::GridViewType& gv, + std::shared_ptr param, + std::shared_ptr gf_water_flux, + std::shared_ptr gf_water_content) + : _gv(gv) + , _param(param) + , _gf_water_flux(gf_water_flux) + , _gf_water_content(gf_water_content) + {} + + /*-----------------------------------------------------------------------*//** + * @brief Evaluation of the effective hydrodynamic dispersion for a given + * entity e in an entity position x + * + * @param[in] e Entity of a grid + * @param[in] x Position in local coordinates with respect the entity + * @param y Effective hydrodynamic dispersion at position x + */ + void evaluate ( const typename Traits::ElementType& e, + const typename Traits::DomainType& x, + typename Traits::RangeType& y) const + { + typename GFWaterContent::Traits::RangeType water_content; + typename GFWaterFlux::Traits::RangeType water_flux; + + // evaluate water flux + _gf_water_flux->evaluate(e,x,water_flux); + // evaluate water flux + _gf_water_content->evaluate(e,x,water_content); + // bind to entity + _param->bind(e); + // evaluate hydrodynamic dispersion + auto y_tensor = _param->hydrodynamic_dispersion_f()(water_flux,water_content); + + using Range = typename Traits::RangeType; + using TensorRange = decltype(y_tensor); + + int count = 0; + y = Range(0.); + + for (int i = 0; i < dimRange; ++i) + for (int j = 0; j < dimRange; ++j, ++count) + if (i _param; + const std::shared_ptr _gf_water_flux; + const std::shared_ptr _gf_water_content; +}; + +} // namespace Dune +} // namespace Dorie + +#endif // DUNE_DORIE_HYDRODYNAMIC_DISPERSION_ADAPTER_HH \ No newline at end of file diff --git a/dune/dorie/model/transport/adapters/peclet.hh b/dune/dorie/model/transport/adapters/peclet.hh new file mode 100644 index 0000000000000000000000000000000000000000..a096590006b24e5cc91556a6e8fab38cafbf4001 --- /dev/null +++ b/dune/dorie/model/transport/adapters/peclet.hh @@ -0,0 +1,108 @@ +#ifndef DUNE_DORIE_PECLET_ADAPTER_HH +#define DUNE_DORIE_PECLET_ADAPTER_HH + +#include + +namespace Dune{ +namespace Dorie{ + +/** + * @brief Converts an interface to match a microscopic peclet number. + * @ingroup TransportParam + * + * @tparam T The base traits + * @tparam Parameter The transport parameter class + * @tparam GFWaterFlux The water flux grid function class + * @tparam GFWaterContent The water content grid function class + */ +template +class PecletAdapter + : public Dune::PDELab::GridFunctionBase< + Dune::PDELab::GridFunctionTraits< + typename T::GridView, + typename T::RangeField, + 1, + typename T::Scalar + > + ,PecletAdapter + > +{ +public: + using Traits = Dune::PDELab::GridFunctionTraits< + typename T::GridView, + typename T::RangeField, + 1, + typename T::Scalar + >; + +private: + using RF = typename T::RangeField; + using Domain = typename T::Domain; + +public: + /*-----------------------------------------------------------------------*//** + * @brief Constructor for the PecletAdapter + * + * @param gv GridView + * @param[in] param The parameter + * @param[in] gf_water_flux The gf water flux + * @param[in] gf_water_content The gf water content + * @param p Parametrization class + */ + PecletAdapter(const typename Traits::GridViewType& gv, + std::shared_ptr param, + std::shared_ptr gf_water_flux, + std::shared_ptr gf_water_content) + : _gv(gv) + , _param(param) + , _gf_water_flux(gf_water_flux) + , _gf_water_content(gf_water_content) + {} + + /*-----------------------------------------------------------------------*//** + * @brief Evaluation of the microscopic peclet number for a given entity + * e in an entity position x + * + * @param[in] e Entity of a grid + * @param[in] x Position in local coordinates with respect the entity + * @param y Microscopic peclet number at position x + */ + void inline evaluate ( const typename Traits::ElementType& e, + const typename Traits::DomainType& x, + typename Traits::RangeType& y) const + { + typename GFWaterFlux::Traits::RangeType water_flux; + typename GFWaterContent::Traits::RangeType water_content; + + // evaluate water flux + _gf_water_flux->evaluate(e,x,water_flux); + // evaluate water flux + _gf_water_content->evaluate(e,x,water_content); + // bind to entity + _param->bind(e); + // evaluate peclet from parameterization + y = _param->peclet_f()(water_flux,water_content); + } + + /*-------------------------------------------------------------------*//** + * @brief Function that returns a grid view valid for this grid + * function + * + * @return Grid view + */ + const typename Traits::GridViewType& getGridView() const + { + return _gv; + } + +private: + const typename Traits::GridViewType _gv; + const std::shared_ptr _param; + const std::shared_ptr _gf_water_flux; + const std::shared_ptr _gf_water_content; +}; + +} // namespace Dune +} // namespace Dorie + +#endif // DUNE_DORIE_PECLET_ADAPTER_HH \ No newline at end of file diff --git a/dune/dorie/model/transport/local_operator_FV.hh b/dune/dorie/model/transport/local_operator_FV.hh index ea2b13b6add20a79c883d18b46bed40ed24ad437..c36fa6fd430f87af2d138478cfe5e5ff88bae12f 100644 --- a/dune/dorie/model/transport/local_operator_FV.hh +++ b/dune/dorie/model/transport/local_operator_FV.hh @@ -1,6 +1,6 @@ // -*- tab-width: 4; indent-tabs-mode: nil -*- -#ifndef DUNE_DORIE_TRANSPORT_OPERATOR_HH -#define DUNE_DORIE_TRANSPORT_OPERATOR_HH +#ifndef DUNE_DORIE_TRANSPORT_OPERATOR_FV_HH +#define DUNE_DORIE_TRANSPORT_OPERATOR_FV_HH #include @@ -36,401 +36,427 @@ namespace Operator { * @f} * * @author Santiago Ospina De Los Ríos - * @date 2018 + * @date 2018-2019 * @ingroup LocalOperators * @ingroup TransportModel * - * @todo Use diffusion coefficient from parameterization - * - * @tparam Boundary Type of the class providing boundary conditions - * @tparam GFWaterFlux Type of a grid function which provides the water - * flux - * @tparam GFWaterContent Type of a grid function which provides the water - * content + * @tparam Boundary Type of the class providing boundary + * conditions + * @tparam GridFunctionWaterFlux Type of a grid function which provides + * the water flux + * @tparam GridFunctionWaterContent Type of a grid function which provides + * the water content */ -template +template class TransportFVSpatialOperator - : public Dune::PDELab::NumericalJacobianVolume< - TransportFVSpatialOperator< - Boundary, - GFWaterFlux, - GFWaterContent - > - > - , public Dune::PDELab::NumericalJacobianSkeleton< - TransportFVSpatialOperator< - Boundary, - GFWaterFlux, - GFWaterContent - > - > - , public Dune::PDELab::NumericalJacobianBoundary< - TransportFVSpatialOperator< - Boundary, - GFWaterFlux, - GFWaterContent - > - > - , public Dune::PDELab::FullSkeletonPattern - , public Dune::PDELab::FullVolumePattern - , public Dune::PDELab::LocalOperatorDefaultFlags - , public Dune::PDELab::InstationaryLocalOperatorDefaultMethods + : public Dune::PDELab::NumericalJacobianVolume< + TransportFVSpatialOperator< + Parameter, + Boundary, + GFWaterFlux, + GFWaterContent + > + > + , public Dune::PDELab::NumericalJacobianSkeleton< + TransportFVSpatialOperator< + Parameter, + Boundary, + GFWaterFlux, + GFWaterContent + > + > + , public Dune::PDELab::NumericalJacobianBoundary< + TransportFVSpatialOperator< + Parameter, + Boundary, + GFWaterFlux, + GFWaterContent + > + > + , public Dune::PDELab::FullSkeletonPattern + , public Dune::PDELab::FullVolumePattern + , public Dune::PDELab::LocalOperatorDefaultFlags + , public Dune::PDELab::InstationaryLocalOperatorDefaultMethods { public: - // Pattern assembly flags - enum { doPatternVolume = true }; - enum { doPatternSkeleton = true }; - - // Residual assembly flags - enum { doAlphaBoundary = true }; - enum { doAlphaSkeleton = true }; + // Pattern assembly flags + enum { doPatternVolume = true }; + enum { doPatternSkeleton = true }; + // Residual assembly flags + enum { doAlphaBoundary = true }; + enum { doAlphaSkeleton = true }; private: - static_assert(std::is_same< - typename GFWaterFlux::Traits::GridViewType, - typename GFWaterContent::Traits::GridViewType>::value, - "TransportFVSpatialOperator: GFWaterFlux and" - "GFWaterContent has to use the same grid view."); - - // Define range type for water flux and water content - using WaterFlux = typename GFWaterFlux::Traits::RangeType; - using WaterContent = typename GFWaterContent::Traits::RangeType; + static_assert(std::is_same< + typename GFWaterFlux::Traits::GridViewType, + typename GFWaterContent::Traits::GridViewType>::value, + "TransportFVSpatialOperator: GFWaterFlux and" + "GFWaterContent has to use the same grid view."); -public: + // Define domain field + using DF = typename GFWaterContent::Traits::DomainFieldType; + using RF = typename GFWaterContent::Traits::RangeFieldType; - /** - * @brief Constructor of TransportFVSpatialOperator - * - * @param boundary The boundary terms - * @param[in] diff_coeff The diffusion coefficient - */ - TransportFVSpatialOperator( - std::shared_ptr boundary, - double diff_coeff - ) : _boundary(boundary) - , _time(0.) - , _diff_coeff(diff_coeff) - { } + // Define range type for water flux and water content + using WaterFlux = typename GFWaterFlux::Traits::RangeType; + using WaterContent = typename GFWaterContent::Traits::RangeType; - /*---------------------------------------------------------------------*//** - * @brief Skeleton integral depending on test and ansatz functions. - * Each face is only visited once since this method is symmetric - * - * @param[in] ig The intersection entity of the grid (inside + outside - * entities) - * @param[in] lfsu_i The inside ansatz local function space - * @param[in] x_i The coefficients of the lfsu_i - * @param[in] lfsv_i The inside test local function space - * @param[in] lfsu_o The outside ansatz local function space - * @param[in] x_o The coefficients of the lfsu_o - * @param[in] lfsv_o The outside test local function space - * @param r_i The view of the residual vector w.r.t lfsu_i - * @param r_o The view of the residual vector w.r.t lfsu_o - * - * @tparam IG The type for ig - * @tparam LFSU The type for lfsu_i and lfsu_o - * @tparam X The type for x_i and x_o - * @tparam LFSV The type for lfsv_i and lfsv_o - * @tparam R The type for r_i and r_o - */ - template - void alpha_skeleton(const IG& ig, - const LFSU& lfsu_i, const X& x_i, const LFSV& lfsv_i, - const LFSU& lfsu_o, const X& x_o, const LFSV& lfsv_o, - R& r_i, R& r_o) const - { - // Get entity dim - // const int dim = IG::Entity::dimension; - - // Get local basis traits from local function space - using LocalBasisTraitsU = typename LFSU::Traits::FiniteElementType:: - Traits::LocalBasisType::Traits; - - // Get range, range field and gradient for trial space - using RangeU = typename LocalBasisTraitsU::RangeType; - // using RangeFieldU = typename LocalBasisTraitsU::RangeFieldType; - // using GradientU = Dune::FieldVector; - - // Assert that polynomial degree is always 0 - assert(lfsu_i.finiteElement().localBasis().order()==0); - assert(lfsu_o.finiteElement().localBasis().order()==0); - assert(lfsv_i.finiteElement().localBasis().order()==0); - assert(lfsv_o.finiteElement().localBasis().order()==0); - - // Define entities for each frame of reference - const auto& entity_f = ig.intersection(); - const auto& entity_i = ig.inside(); - const auto& entity_o = ig.outside(); - - // Get geometries - auto geo_f = entity_f.geometry(); - auto geo_i = entity_i.geometry(); - auto geo_o = entity_o.geometry(); - - // Get geometry of intersection in local coord. of neighbor entities - auto geo_in_i = entity_f.geometryInInside(); - - // Get volume of entity - const auto volume_f = geo_f.volume(); - - // Get normal - auto normal_f = entity_f.centerUnitOuterNormal(); - - // Entity centers in references elements - const auto& ref_el_i = referenceElement(geo_i); - const auto& ref_el_o = referenceElement(geo_o); - const auto& center_position_i = ref_el_i.position(0,0); - const auto& center_position_o = ref_el_o.position(0,0); - - // Get center position of the face w.r.t inside entity - auto face_position_i = geo_in_i.center(); - - // Evaluate water content - WaterContent water_content_i,water_content_o; - _gf_water_content->evaluate(entity_i,center_position_i,water_content_i); - _gf_water_content->evaluate(entity_o,center_position_o,water_content_o); - - // Eval flux field assume H(div) velocity field => may choose any side - WaterFlux water_flux; - _gf_water_flux->evaluate(entity_i, face_position_i, water_flux); - auto water_flux_n = water_flux*normal_f; - // auto water_vel_i = water_flux_n*water_content_i; - // auto water_vel_o = water_flux_n*water_content_o; - - // Compute the effective hydrodynamic dispersion - - // bind parameterization and retrieve functions - // _param->bind(entity_i); - // const auto hydr_disp_func_i = _param->hydrodynamic_dispersion_f(); - // auto hydr_disp_i = hydr_disp_func_i(water_vel_i,water_content_i); - - // _param->bind(entity_o); - // const auto hydr_disp_func_o = _param->hydrodynamic_dispersion_f(); - // auto hydr_disp_o = hydr_disp_func_o(water_vel_o,water_content_o); - - // Evaluate diffusion coeff. from both sides - // GradientU diff_coeff_i, diff_coeff_o; - // hydr_disp_i.mv(normal_f,diff_coeff_i); - // hydr_disp_o.mv(normal_f,diff_coeff_o); - - auto diff_coeff_i = _diff_coeff * water_content_i; // FIXME - auto diff_coeff_o = _diff_coeff * water_content_o; // FIXME - - // Take harmonic average - auto diff_coeff_f = 2.0/( 1.0/(diff_coeff_i + 1E-30) - + 1.0/(diff_coeff_o + 1E-30) ); - - // Inside/outside solute value - RangeU u_i = x_i(lfsu_i,0); - RangeU u_o = x_o(lfsu_o,0); - - // Upwinding - const auto& u_upwind = (water_flux_n>=0) ? u_i : u_o; - - // Entity centers in global coordinates - auto center_position_i_g = geo_i.global(center_position_i); - auto center_position_o_g = geo_o.global(center_position_o); - - // Distance between the two entity centers - center_position_i_g -= center_position_o_g; - auto distance = center_position_i_g.two_norm(); - - // Finite difference of u between the two entities - RangeU dudn = (u_o-u_i)/distance; - - // Solute flux in normal direction w.r.t the intersection - auto soulte_flux_n = u_upwind*water_flux_n - diff_coeff_f*dudn; - - // Symmetric contribution to residual on inside element - r_i.accumulate(lfsv_i, 0, soulte_flux_n*volume_f ); - r_o.accumulate(lfsv_o, 0, -soulte_flux_n*volume_f ); - } + // Extract world dimension + enum {dim = GFWaterContent::Traits::dimDomain}; + using Tensor = Dune::FieldMatrix; +public: - /*---------------------------------------------------------------------*//** - * @brief Boundary integral depending on test and ansatz functions. - * - * @param[in] ig The intersection entity of the grid (inside + outside - * entities) - * @param[in] lfsu_i The inside ansatz local function space - * @param[in] x_i The coefficients of the lfsu_i - * @param[in] lfsv_i The inside test local function space - * @param r_i The view of the residual vector w.r.t lfsu_i - * - * @tparam IG The type for ig - * @tparam LFSU The type for lfsu_i - * @tparam X The type for x_i - * @tparam LFSV The type for lfsv_i - * @tparam R The type for r_i - */ - template - void alpha_boundary (const IG& ig, const LFSU& lfsu_i, const X& x_i, - const LFSV& lfsv_i, R& r_i) const + /** + * @brief Constructor of TransportFVSpatialOperator + * + * @param[in] param The parameters + * @param[in] boundary The boundary terms + * @param[in] gf_water_flux The grid function of water flux + * @param[in] gf_water_content The grid function of water content + * @param[in] dirichlet_mode The dirichlet mode + */ + TransportFVSpatialOperator( + std::shared_ptr param, + std::shared_ptr boundary, + DirichletMode::Type dirichlet_mode = DirichletMode::SoluteConcentration + ) : _param(param) + , _boundary(boundary) + , _dirichlet_mode(dirichlet_mode) + { + assert(_gf_water_flux); + assert(_gf_water_content); + } + + /*---------------------------------------------------------------------*//** + * @brief Skeleton integral depending on test and ansatz functions. + * Each face is only visited once since this method is symmetric + * + * @param[in] ig The intersection entity of the grid (inside + outside + * entities) + * @param[in] lfsu_i The inside ansatz local function space + * @param[in] x_i The coefficients of the lfsu_i + * @param[in] lfsv_i The inside test local function space + * @param[in] lfsu_o The outside ansatz local function space + * @param[in] x_o The coefficients of the lfsu_o + * @param[in] lfsv_o The outside test local function space + * @param r_i The view of the residual vector w.r.t lfsu_i + * @param r_o The view of the residual vector w.r.t lfsu_o + * + * @tparam IG The type for ig + * @tparam LFSU The type for lfsu_i and lfsu_o + * @tparam X The type for x_i and x_o + * @tparam LFSV The type for lfsv_i and lfsv_o + * @tparam R The type for r_i and r_o + */ + template + void alpha_skeleton(const IG& ig, + const LFSU& lfsu_i, const X& x_i, const LFSV& lfsv_i, + const LFSU& lfsu_o, const X& x_o, const LFSV& lfsv_o, + R& r_i, R& r_o) const + { + // Get entity dim + const int dim = IG::Entity::dimension; + + // Get local basis traits from local function space + using LocalBasisTraitsU = typename LFSU::Traits::FiniteElementType:: + Traits::LocalBasisType::Traits; + + // Get range, range field and gradient for trial space + using RangeU = typename LocalBasisTraitsU::RangeType; + using RangeFieldU = typename LocalBasisTraitsU::RangeFieldType; + using GradientU = Dune::FieldVector; + + // Assert that polynomial degree is always 0 + assert(lfsu_i.finiteElement().localBasis().order()==0); + assert(lfsu_o.finiteElement().localBasis().order()==0); + assert(lfsv_i.finiteElement().localBasis().order()==0); + assert(lfsv_o.finiteElement().localBasis().order()==0); + + // Define entities for each frame of reference + const auto& entity_f = ig.intersection(); + const auto& entity_i = ig.inside(); + const auto& entity_o = ig.outside(); + + // Get geometries + auto geo_f = entity_f.geometry(); + auto geo_i = entity_i.geometry(); + auto geo_o = entity_o.geometry(); + + // Get geometry of intersection in local coord. of neighbor entities + auto geo_in_i = entity_f.geometryInInside(); + + // Get volume of entities + const auto volume_f = geo_f.volume(); + + // Get normal + auto normal_f = entity_f.centerUnitOuterNormal(); + + // Entity centers in references elements + const auto& ref_el_i = referenceElement(geo_i); + const auto& ref_el_o = referenceElement(geo_o); + const auto& center_position_i = ref_el_i.position(0,0); + const auto& center_position_o = ref_el_o.position(0,0); + + // Get center position of the face w.r.t inside entity + auto face_position_i = geo_in_i.center(); + + // Evaluate water content + WaterContent water_content_i,water_content_o; + _gf_water_content->evaluate(entity_i,center_position_i,water_content_i); + _gf_water_content->evaluate(entity_o,center_position_o,water_content_o); + + // Eval flux field assume H(div) velocity field => may choose any side + WaterFlux water_flux; + _gf_water_flux->evaluate(entity_i, face_position_i, water_flux); + auto water_flux_n = water_flux*normal_f; + + // Compute the effective hydrodynamic dispersion + // bind parameterization and retrieve functions + _param->bind(entity_i); + const auto hydr_disp_func_i = _param->hydrodynamic_dispersion_f(); + auto hydr_disp_i = hydr_disp_func_i(water_flux,water_content_i); + + _param->bind(entity_o); + const auto hydr_disp_func_o = _param->hydrodynamic_dispersion_f(); + auto hydr_disp_o = hydr_disp_func_o(water_flux,water_content_o); + + // Evaluate diffusion coeff. from both sides + GradientU diff_coeff_i, diff_coeff_o; + hydr_disp_i.mv(normal_f,diff_coeff_i); + hydr_disp_o.mv(normal_f,diff_coeff_o); + + diff_coeff_i *= water_content_i; + diff_coeff_o *= water_content_o; + + // Take harmonic average + auto diff_coeff_f = 2.0/( 1.0/(diff_coeff_i*normal_f + 1E-30) + + 1.0/(diff_coeff_o*normal_f + 1E-30) ); + + // Inside/outside solute value + RangeU u_i = x_i(lfsu_i,0); + RangeU u_o = x_o(lfsu_o,0); + + // Upwinding + const auto& u_upwind = (water_flux_n>=0) ? u_i : u_o; + + // Entity centers in global coordinates + auto center_position_i_g = geo_i.global(center_position_i); + auto center_position_o_g = geo_o.global(center_position_o); + + // Distance between the two entity centers + center_position_i_g -= center_position_o_g; + auto distance = center_position_i_g.two_norm(); + + // Finite difference of u between the two entities + RangeU dudn = (u_o-u_i)/distance; + + // Solute flux in normal direction w.r.t the intersection + auto soulte_flux_n = u_upwind*water_flux_n - diff_coeff_f*dudn; + + // Symmetric contribution to residual on inside element + r_i.accumulate(lfsv_i, 0, soulte_flux_n*volume_f ); + r_o.accumulate(lfsv_o, 0, -soulte_flux_n*volume_f ); + } + + + /*---------------------------------------------------------------------*//** + * @brief Boundary integral depending on test and ansatz functions. + * + * @param[in] ig The intersection entity of the grid (inside + outside + * entities) + * @param[in] lfsu_i The inside ansatz local function space + * @param[in] x_i The coefficients of the lfsu_i + * @param[in] lfsv_i The inside test local function space + * @param r_i The view of the residual vector w.r.t lfsu_i + * + * @tparam IG The type for ig + * @tparam LFSU The type for lfsu_i + * @tparam X The type for x_i + * @tparam LFSV The type for lfsv_i + * @tparam R The type for r_i + */ + template + void alpha_boundary (const IG& ig, const LFSU& lfsu_i, const X& x_i, + const LFSV& lfsv_i, R& r_i) const + { + // Get entity dim + const int dim = IG::Entity::dimension; + + // Get local basis traits from local function space + using LocalBasisTraitsU = typename LFSU::Traits::FiniteElementType:: + Traits::LocalBasisType::Traits; + + // Get range, range field and gradient for trial space + using RangeU = typename LocalBasisTraitsU::RangeType; + using RangeFieldU = typename LocalBasisTraitsU::RangeFieldType; + using GradientU = Dune::FieldVector; + + // Assert that polynomial degree is always 0 + assert(lfsu_i.finiteElement().localBasis().order()==0); + assert(lfsv_i.finiteElement().localBasis().order()==0); + + // Get entities + const auto& entity_f = ig.intersection(); + const auto& entity_i = ig.inside(); + + // Get geometries + auto geo_f = entity_f.geometry(); + const auto& ref_el_f = referenceElement(geo_f); + const auto& center_position_f = ref_el_f.position(0,0); + + // Get normal + auto normal_f = entity_f.centerUnitOuterNormal(); + + // Get boundary condition type + auto bc = _boundary->bc(entity_f,center_position_f,_time); + + // Get geometry of intersection in local coord. of neighbor entity + auto geo_in_i = ig.geometryInInside(); + auto face_position_i = geo_in_i.center(); + + // Evaluate convective term (i.e water flux) + WaterFlux water_flux; + _gf_water_flux->evaluate(entity_i, face_position_i, water_flux); + auto water_flux_n = water_flux*normal_f; + + // bind parameterization and retrieve functions + _param->bind(entity_i); + const auto hydr_disp_func_i = _param->hydrodynamic_dispersion_f(); + + // Face volume for integration + auto volume_f = geo_f.volume(); + + if (BoundaryCondition::isOutflow(bc)) { - // Get entity dim - // const int dim = IG::Entity::dimension; - - // Get local basis traits from local function space - using LocalBasisTraitsU = typename LFSU::Traits::FiniteElementType:: - Traits::LocalBasisType::Traits; - - // Get range, range field and gradient for trial space - using RangeU = typename LocalBasisTraitsU::RangeType; - // using RangeFieldU = typename LocalBasisTraitsU::RangeFieldType; - // using GradientU = Dune::FieldVector; - - // Assert that polynomial degree is always 0 - assert(lfsu_i.finiteElement().localBasis().order()==0); - assert(lfsv_i.finiteElement().localBasis().order()==0); - - // Get entities - const auto& entity_f = ig.intersection(); - const auto& entity_i = ig.inside(); - - // Get geometries - auto geo_f = entity_f.geometry(); - const auto& ref_el_f = referenceElement(geo_f); - const auto& center_position_f = ref_el_f.position(0,0); - - // Get normal - auto normal_f = entity_f.centerUnitOuterNormal(); - - // Get boundary condition type - auto bc = _boundary->bc(entity_f,center_position_f,_time); - - // Get geometry of intersection in local coord. of neighbor entity - auto geo_in_i = ig.geometryInInside(); - auto face_position_i = geo_in_i.center(); - - // Evaluate convective term (i.e water flux) - WaterFlux water_flux; - _gf_water_flux->evaluate(entity_i, face_position_i, water_flux); - auto water_flux_n = water_flux*normal_f; - // auto water_vel_i = water_flux_n*water_content_i; - - // bind parameterization and retrieve functions - // _param->bind(entity_i); - // const auto hydr_disp_func_i = _param->hydrodynamic_dispersion_f(); - - // Face volume for integration - auto volume_f = geo_f.volume(); - - if (BoundaryCondition::isOutflow(bc)) - { - // Check that water is going out - if (Dune::FloatCmp::lt(water_flux_n,0.)) - return; - - // Inside unknown value - RangeU u = x_i(lfsu_i,0); - - // Eval flux boundary condition - auto o = _boundary->o(entity_f,center_position_f,_time); - - // convection term - auto conv_num_flux = (u * water_flux_n + o) * volume_f; - - // integrate o - r_i.accumulate(lfsu_i,0,conv_num_flux); - } - else if (BoundaryCondition::isDirichlet(bc)) - { - // Only allow dirichlet values for influx cases - if (Dune::FloatCmp::gt(water_flux_n,0.)) - return; + // Check that water is going out + if (Dune::FloatCmp::lt(water_flux_n,0.)) + return; - // Get center of entity in global coordinates - auto geo_i = entity_i.geometry(); - auto center_position_i_g = geo_i.center(); + // Inside unknown value + RangeU u = x_i(lfsu_i,0); - // Face center in global coordinates - auto center_position_f_g = geo_f.center(); + // Eval flux boundary condition + auto o = _boundary->o(entity_f,center_position_f,_time); - // Compute distance of these two points - center_position_i_g -= center_position_f_g; - auto distance = center_position_i_g.two_norm(); + // convection term + auto conv_num_flux = (u * water_flux_n + o) * volume_f; - // Get geometry of intersection in local coord. of neighbor entity - auto geo_in_i = ig.geometryInInside(); - auto face_position_i = geo_in_i.center(); - - // Evaluate water content - WaterContent water_content_i; - _gf_water_content->evaluate(entity_i, face_position_i, water_content_i); - - // Evaluate Dirichlet condition - auto g = _boundary->g(entity_f,center_position_f,_time); - - // Compute the effective hydrodynamic dispersion - // auto hydr_disp_i = hydr_disp_func_i(water_vel_i,water_content_i); - - // Calculate diffusion coefficient from inside - // GradientU diff_coeff_i; - // hydr_disp_i.mv(normal_f,diff_coeff_i); - - auto diff_coeff_i = _diff_coeff; // FIXME - - // Inside unknown value - RangeU u = x_i(lfsu_i,0); - - // Solute flux in normal direction w.r.t the intersection - auto soulte_flux_n = (g*water_flux_n) - -(diff_coeff_i*water_content_i*(g-u)/distance); - - // Contribution to residual from Dirichlet boundary - r_i.accumulate(lfsv_i, 0, soulte_flux_n*volume_f); - - } - else if (BoundaryCondition::isNeumann(bc)) - { - // Solute flux in normal direction w.r.t the intersection - auto soulte_flux_n = _boundary->j(entity_f,center_position_f,_time); - - r_i.accumulate(lfsv_i, 0, soulte_flux_n*volume_f ); - } + // integrate o + r_i.accumulate(lfsu_i,0,conv_num_flux); } - - /*---------------------------------------------------------------------*//** - * @brief Sets the time. - * - * @param[in] t time of type RangeField - * - * @tparam RangeField type of the range field - */ - template - void setTime (RangeField t) + else if (BoundaryCondition::isDirichlet(bc)) { - _time = t; - _gf_water_content->setTime(t); - _gf_water_flux->setTime(t); - } + // Only allow dirichlet values for influx cases + if (Dune::FloatCmp::gt(water_flux_n,0.)) + return; + + // Get center of entity in global coordinates + auto geo_i = entity_i.geometry(); + auto center_position_i_g = geo_i.center(); + + // Face center in global coordinates + auto center_position_f_g = geo_f.center(); + + // Compute distance of these two points + center_position_i_g -= center_position_f_g; + auto distance = center_position_i_g.two_norm(); + + // Get geometry of intersection in local coord. of neighbor entity + auto geo_in_i = ig.geometryInInside(); + auto face_position_i = geo_in_i.center(); + + // Evaluate water content + WaterContent water_content_i; + _gf_water_content->evaluate(entity_i, face_position_i, water_content_i); + + // Evaluate Dirichlet condition + auto g = _boundary->g(entity_f,center_position_f,_time); + + // Convert input to total solute if needed + if (_dirichlet_mode == DirichletMode::TotalSolute) + { + if (Dune::FloatCmp::gt(water_content_i,WaterContent{0.})) + g /=water_content_i; + else + g = 0.; + } + + // Compute the effective hydrodynamic dispersion + auto hydr_disp_i = hydr_disp_func_i(water_flux,water_content_i); + + // Calculate diffusion coefficient from inside + GradientU diff_coeff_i; + hydr_disp_i.mv(normal_f,diff_coeff_i); + + // Inside unknown value + RangeU u = x_i(lfsu_i,0); + + // Solute flux in normal direction w.r.t the intersection + auto soulte_flux_n = (g*water_flux_n) + -(diff_coeff_i*water_content_i*(g-u)/distance); + + // Contribution to residual from Dirichlet boundary + r_i.accumulate(lfsv_i, 0, soulte_flux_n*volume_f); - /*---------------------------------------------------------------------*//** - * @brief Sets the the water content grid function. - * - * @param[in] gf_water_content The water content grid function - */ - void set_water_content(std::shared_ptr gf_water_content) - { - _gf_water_content = gf_water_content; } - - /*---------------------------------------------------------------------*//** - * @brief Sets the the water flux grid function. - * - * @param[in] gf_water_content The water flux grid function - */ - void set_water_flux(std::shared_ptr gf_water_flux) + else if (BoundaryCondition::isNeumann(bc)) { - _gf_water_flux = gf_water_flux; + // Solute flux in normal direction w.r.t the intersection + auto soulte_flux_n = _boundary->j(entity_f,center_position_f,_time); + + r_i.accumulate(lfsv_i, 0, soulte_flux_n*volume_f ); } + } + + /*---------------------------------------------------------------------*//** + * @brief Sets the time. + * + * @param[in] t time of type RangeField + * + * @tparam RangeField type of the range field + */ + template + void setTime (RangeField t) + { + _time = t; + _gf_water_content->setTime(t); + _gf_water_flux->setTime(t); + } + + double getTime () const + { + return _time; + } + + /*---------------------------------------------------------------------*//** + * @brief Sets the the water content grid function. + * + * @param[in] gf_water_content The water content grid function + */ + void set_water_content(std::shared_ptr gf_water_content) + { + _gf_water_content = gf_water_content; + } + + /*---------------------------------------------------------------------*//** + * @brief Sets the the water flux grid function. + * + * @param[in] gf_water_content The water flux grid function + */ + void set_water_flux(std::shared_ptr gf_water_flux) + { + _gf_water_flux = gf_water_flux; + } private: - const std::shared_ptr _boundary; - std::shared_ptr _gf_water_flux; - std::shared_ptr _gf_water_content; - double _time; - double _diff_coeff; + + const std::shared_ptr _param; + const std::shared_ptr _boundary; + std::shared_ptr _gf_water_flux; + std::shared_ptr _gf_water_content; + const DirichletMode::Type _dirichlet_mode; + double _time; }; /** @@ -451,13 +477,10 @@ private: * @f} * * @author Santiago Ospina De Los Ríos - * @date 2018 + * @date 2018-2019 * @ingroup LocalOperators * @ingroup TransportModel * - * @bug The water content is not being updated by the set time to the - * right state of the grid function. - * * @tparam GFWaterContent Type of a grid function which provides the water * content */ @@ -474,11 +497,11 @@ template { public: - // Pattern assembly flags - enum { doPatternVolume = true }; + // Pattern assembly flags + enum { doPatternVolume = true }; - // Eesidual assembly flags - enum { doAlphaVolume = true }; + // Eesidual assembly flags + enum { doAlphaVolume = true }; public: @@ -486,75 +509,75 @@ public: : _time(0.) { } - /*---------------------------------------------------------------------*//** - * @brief Volume integral depending on test and ansatz functions - * - * @param[in] eg THe entity of the grid - * @param[in] lfsu The ansatz local function space - * @param[in] x The coefficients of the lfsu - * @param[in] lfsv The test local function space - * @param r The view of the residual vector w.r.t lfsu - * - * @tparam EG The type for eg - * @tparam LFSU The type for lfsu - * @tparam X The type for x - * @tparam LFSV The type for lfsv - * @tparam R The type for r - */ - template - void alpha_volume (const EG& eg, const LFSU& lfsu, const X& x, - const LFSV& lfsv, R& r) const - { - // Get local basis traits from local function space - using LocalBasisTraitsU = typename LFSU::Traits::FiniteElementType:: - Traits::LocalBasisType::Traits; - // Get range for trial space - using RangeU = typename LocalBasisTraitsU::RangeType; - - RangeU u = x(lfsu,0); - - // entity geometry - auto geo = eg.geometry(); - - // inside cell center - const auto& ref_el = referenceElement(geo); - const auto& center_position = ref_el.position(0,0); - - typename GFWaterContent::Traits::RangeType water_content; - - _gf_water_content->evaluate(eg.entity(),center_position,water_content); - - // update residual - r.accumulate(lfsv ,0 , water_content*u*geo.volume()); - } - - /*---------------------------------------------------------------------*//** - * @brief Sets the time. - * - * @param[in] t time of type RangeField - * - * @tparam RangeField type of the range field - */ - template - void setTime (RF t) - { - _time = t; - _gf_water_content->setTime(t); - } - - /*---------------------------------------------------------------------*//** - * @brief Sets the the water content grid function. - * - * @param[in] gf_water_content The water content grid function - */ - void set_water_content(std::shared_ptr gf_water_content) - { - _gf_water_content = gf_water_content; - } + /*---------------------------------------------------------------------*//** + * @brief Volume integral depending on test and ansatz functions + * + * @param[in] eg THe entity of the grid + * @param[in] lfsu The ansatz local function space + * @param[in] x The coefficients of the lfsu + * @param[in] lfsv The test local function space + * @param r The view of the residual vector w.r.t lfsu + * + * @tparam EG The type for eg + * @tparam LFSU The type for lfsu + * @tparam X The type for x + * @tparam LFSV The type for lfsv + * @tparam R The type for r + */ + template + void alpha_volume (const EG& eg, const LFSU& lfsu, const X& x, + const LFSV& lfsv, R& r) const + { + // Get local basis traits from local function space + using LocalBasisTraitsU = typename LFSU::Traits::FiniteElementType:: + Traits::LocalBasisType::Traits; + // Get range for trial space + using RangeU = typename LocalBasisTraitsU::RangeType; + + RangeU u = x(lfsu,0); + + // entity geometry + auto geo = eg.geometry(); + + // inside cell center + const auto& ref_el = referenceElement(geo); + const auto& center_position = ref_el.position(0,0); + + typename GFWaterContent::Traits::RangeType water_content; + + _gf_water_content->evaluate(eg.entity(),center_position,water_content); + + // update residual + r.accumulate(lfsv ,0 , water_content*u*geo.volume()); + } + + /*---------------------------------------------------------------------*//** + * @brief Sets the time. + * + * @param[in] t time of type RangeField + * + * @tparam RangeField type of the range field + */ + template + void setTime (RF t) + { + _time = t; + _gf_water_content->setTime(t); + } + + /*---------------------------------------------------------------------*//** + * @brief Sets the the water content grid function. + * + * @param[in] gf_water_content The water content grid function + */ + void set_water_content(std::shared_ptr gf_water_content) + { + _gf_water_content = gf_water_content; + } private: - std::shared_ptr _gf_water_content; - double _time; + std::shared_ptr _gf_water_content; + double _time; }; } // namespace Operator @@ -562,4 +585,4 @@ private: } // namespace Dune -#endif // DUNE_DORIE_TRANSPORT_OPERATOR_HH +#endif // DUNE_DORIE_TRANSPORT_OPERATOR_FV_HH diff --git a/dune/dorie/model/transport/modules.dox b/dune/dorie/model/transport/modules.dox index abead292be65f6adc0bd88dcc57f483c3fafa908..4b49d02cf6bc9f802e379e15ea5a44c2d85a5823 100644 --- a/dune/dorie/model/transport/modules.dox +++ b/dune/dorie/model/transport/modules.dox @@ -29,10 +29,96 @@ \Gamma_D @f} - It implemented with a finite volume scheme and is written so that could + It's implemented with a finite volume scheme and is written so that could be run independently of the @ref RichardsModel model (e.g. with an stationary case), or fully coupled with a transient case with fluxes provided by the Dune::Dorie::RichardsSimulation class. @} +@defgroup TransportParam Parameterization +@{ +@ingroup TransportModel + +@brief Representing dispersion properties of the solute in the porous medium + +## Overview + +The parameterization condenses the subscale physics from below the +REV scale and represents the effective dispersion process in the far-field +regimme of the transport process. The main task of the parameterization in DORiE +is to estimate the hydrodynamic dispersion tensor. However, it also estimates +the microscopic péclet number. Several options (and their combinations) are +available. + +### Parameterization + +The microscopic peclet number is directly implemented in the base class of the +parameterization Dune::Dorie::Parameterization::Transport::peclet_f(). + +The hydrodynamic dispersion tensor is implemented as a polyphormic object in +Dune::Dorie::Parameterization::Transport. It can take the form of three +different parameterizations: + + 1. **Constant hydrodynamic dispersion tensor** + As the name reveals, it simply requires a constant tensor + Dune::Dorie::Parameterization::ConstHydrodynamicDispersion. + + 2. **Power Law** + In this case, the diagonal of the tensor depends on the microscopic peclet + number, the molecular diffusion, and two user-dependent parameters + Dune::Dorie::Parameterization::PowerLawDispersion. + + 3. **Superposition** + The superposition case assumes that the linear combination of different + dispersion processes lead to the desired hydrodynamic dispersion tensor. + Since the concept is not tied to any particular parameterization, the + individual processes are implemented as well as polyphormic objects. + In this case, there are implemented the two most fundamental processes at + pore scale + Dune::Dorie::Parameterization::HydrodynamicDispersionSuperposition: + + a. *Effective diffusion coefficient* + Dune::Dorie::Parameterization::EffectiveDiffusion + + + Constant coefficient + Dune::Dorie::Parameterization::ConstEffectiveDiffusion + + Milligtong Quirk I + Dune::Dorie::Parameterization::MillingtonQuirk1 + + Milligtong Quirk II + Dune::Dorie::Parameterization::MillingtonQuirk2 + + b. *Effective hydromechanic dispersion tensor* + Dune::Dorie::Parameterization::EffectiveHydromechanicDispersion + + + Constant effective hydromechanic dispersion tensor + Dune::Dorie::Parameterization::ConstEffectiveHydromechanicDispersion + + Hydromechanic dispersion for isotropic media + Dune::Dorie::Parameterization::IsotropicEffectiveHydromechanicDispersion + +Storing this data, providing grid-to-data mappings, and access to this data +are tasks of the Dune::Dorie::TransportParameters interface, that is directly +used in the @ref LocalOperators. + +Special care has to be taken when modifying repeated parameters (e.g. molecular +diffusion) or when storing shared pointers in polyhpromic objects. If not +treated carefully, some parameters may be not reachable anymore with the +parameters() methods or be doubly created when using the method clone(). + +@todo Add scaling factors to porosity values + +### Adapters + +The adapter for the hydrodynamic dispersion tensor is rather specific: + ++ Visualization of glyph tensors in Paraview is only avaiable for 9 component + data (3D tensor). ++ The VTK writer of dune-grid can only write data that is indexable once. H + Hence, the 3x3 tensor has to be flattened. ++ The usual Grid Function to VTK Function adapter in dune-pdelab does not allow + to write more than 3 components of data. Hence we write our own Grid Function + to dune-function adapter (Dune::Dorie::VTKGridFunctionAdapter) which is + natively accepted by the current vtk writers. + +@} + **/ diff --git a/dune/dorie/model/transport/parameterization/factory.hh b/dune/dorie/model/transport/parameterization/factory.hh new file mode 100644 index 0000000000000000000000000000000000000000..47134486335b534024beb237179eb92a8b6f3739 --- /dev/null +++ b/dune/dorie/model/transport/parameterization/factory.hh @@ -0,0 +1,166 @@ +#ifndef DUNE_DORIE_PARAM_TRANSPORT_FACTORY_HH +#define DUNE_DORIE_PARAM_TRANSPORT_FACTORY_HH + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include + +namespace Dune { +namespace Dorie { +namespace Parameterization { + +/*-------------------------------------------------------------------------*//** + * @brief Factory of Transport parameters + * + * @tparam Traits The base traits + */ +template +struct TransportFactory + : public ParameterizationFactory> +{ + + /// Return a default-initialized parameterization object + /** \param type The type indicating the parameterization implementation + * \param name The name of the parameterization object (layer type) + */ + std::shared_ptr> + selector( + const YAML::Node& type_node, + const std::string& name) const override + { + auto log = Dorie::get_logger(log_transport); + + std::shared_ptr> hd_dips; + + uint type; + + try { + mu::Parser parser; + parser.DefineConst(ConstHydrodynamicDispersionParam::type, 1); + parser.DefineConst(PowerLawDispersion::type, 2); + + parser.DefineConst(ConstEffectiveDiffusionParam::type, 10); + parser.DefineConst(MillingtonQuirk1::type, 20); + parser.DefineConst(MillingtonQuirk2::type, 30); + + parser.DefineConst(ConstEffectiveHydromechanicDispersionParam::type, 100); + parser.DefineConst(IsotropicEffectiveHydromechanicDispersion::type, 200); + + parser.SetExpr(type_node["type"].as()); + type = parser.Eval(); + + } catch (mu::Parser::exception_type& e) { + log->error("Evaluating transport parameterization type failed:"); + log->error(" Parsed expression: {}",e.GetExpr()); + log->error(" Token: {}",e.GetToken()); + log->error(" Error position: {}",e.GetPos()); + log->error(" Error code: {}",e.GetCode()); + log->error(" Error message: {}",e.GetMsg()); + DUNE_THROW(IOError, "Error evaluating transport parameterization type"); + } + + // get digit x from right to left of a uint. + // uint ....... xxxxxx + // ^^ ^^ + // digit ...... 54..10 + auto get_digit = [](uint number, int x) + {return (uint)(number*std::pow(10,-1*x))%10;}; + + if (type == 1) { + hd_dips = std::make_shared>(name); + } else if (type == 2) { + hd_dips = std::make_shared>(name); + } else if (get_digit(type,0) == 0) { + std::shared_ptr> eff_diff; + + // Select effective diffusion according to digit 1 of 'type' + uint type_1 = get_digit(type,1); + + // Build parametrization for effective diffusion + if (type_1 == 1) { + eff_diff = std::make_shared>(name); + } else if (type_1 == 2) { + eff_diff = std::make_shared>(name); + } else if (type_1 == 3) { + eff_diff = std::make_shared>(name); + } else { + log->error("Effective Diffusion parameterization '{}' " + "has unknown type '{}'", + name, type_node["type"].as()); + DUNE_THROW(IOError, "Unknown Transport parameterization type"); + } + + // Select type for effective hydromechanic dispersion + std::shared_ptr> eff_hm_disp; + + // Select effective hydromechanic dispersion according to digit 2 of 'type' + uint type_2 = get_digit(type,2); + + // Build parametrization for effective hydromechanic dispersion + if (type_2 == 1) { + eff_hm_disp = std::make_shared>(name); + } else if (type_2 == 2) { + + // Select type for longitudinal and transverse dispersivity + std::shared_ptr> lambda_l; + std::shared_ptr> lambda_t; + + // Build parametrization for longitudinal dispersivity + // if (...) { + lambda_l = std::make_shared>(name); + // } else { + // log->error("Longitudinal Dispersivity parameterization '{}' " + // "has unknown type '{}'", + // name, long_disp_type); + // DUNE_THROW(IOError, "Unknown Transport parameterization type"); + // } + + // Build parametrization for transverse dispersivity + // if (...) { + lambda_t = std::make_shared>(name); + // } else { + // log->error("Transverse Dispersivity parameterization '{}' " + // "has unknown type '{}'", + // name, trans_disp_type); + // DUNE_THROW(IOError, "Unknown Transport parameterization type"); + // } + + eff_hm_disp = std::make_shared>(name,lambda_l,lambda_t); + } else { + log->error("Effective Hydromechanic Dispersion parameterization '{}' " + "has unknown type '{}'", + name, type_node["type"].as()); + DUNE_THROW(IOError, "Unknown Transport parameterization type"); + } + + hd_dips = std::make_shared>(name,eff_diff,eff_hm_disp); + } else { + log->error("Transport parameterization '{}' has unknown type '{}'", + name, type_node["type"].as()); + DUNE_THROW(IOError, "Unknown Transport parameterization type"); + } + return hd_dips; + } +}; + +} // namespace Parameterization +} // namespace Dune +} // namespace Dorie + +#endif // DUNE_DORIE_PARAM_TRANSPORT_FACTORY_HH \ No newline at end of file diff --git a/dune/dorie/model/transport/parameterization/hydrodynamic_dispersion/const.hh b/dune/dorie/model/transport/parameterization/hydrodynamic_dispersion/const.hh new file mode 100644 index 0000000000000000000000000000000000000000..9416e5f0790db990aaf31075d257e003823ad03f --- /dev/null +++ b/dune/dorie/model/transport/parameterization/hydrodynamic_dispersion/const.hh @@ -0,0 +1,181 @@ +#ifndef DUNE_DORIE_PARAM_CONST_HYDRODYNAMIC_DISPERSION_HH +#define DUNE_DORIE_PARAM_CONST_HYDRODYNAMIC_DISPERSION_HH + +#include +#include +#include + +#include + +namespace Dune { +namespace Dorie { +namespace Parameterization { + +/*-------------------------------------------------------------------------*//** + * @brief Class for constant hydrodynamic dispersion. + * + * @ingroup TransportParam + * @author Santiago Ospina De Los Ríos + * @date 2019 + * + * @tparam Traits the base traits + */ +template +class ConstHydrodynamicDispersionParam : + public Transport +{ +private: + using Tensor = typename Traits::Tensor; + using Base = Transport; + +public: + using HydrodynamicDispersion = typename Base::HydrodynamicDispersion; + using WaterFlux = typename Base::WaterFlux; + using WaterContent = typename Base::WaterContent; + + /// Parameter defining the constant hydrodynamic dispersion tensor + struct ConstHydrodynamicDispersion + { + Tensor value; + inline static const std::string name = "hydrodynamic_disp"; + }; + + /// The name of this parameterization type (This is not the instance name) + static inline std::string type = "Dhd_const"; + + ConstHydrodynamicDispersion _const_hydrodynamic_dispersion; + +public: + + /// Construct with default-initialized parameters + /** \param name The name associated with this soil layer + */ + ConstHydrodynamicDispersionParam (const std::string name) : + Transport(name) + { } + + /// Construct from a tuple of parameters + /** \param name The name associated with this soil layer + * \param parameters Tuple of parameters to use in this parameterization + */ + template + ConstHydrodynamicDispersionParam ( + const std::string name, + const std::tuple parameters) : + Base(name, parameters), + _const_hydrodynamic_dispersion(std::get(parameters)) + { } + + /// Add default destructor to clarify override + ~ConstHydrodynamicDispersionParam () override = default; + + /// Return the hydrodynamic dispersion + /** {WaterFlux,WaterContent} -> HydrodynamicDispersion + */ + std::function + hydrodynamic_dispersion_f () const override + { + return [this](const WaterFlux, const WaterContent) { + return HydrodynamicDispersion{_const_hydrodynamic_dispersion.value}; + }; + } + + /// Return a map of parameter names and values for manipulation + /** \return Map: Parameter name, parameter value in this object + */ + std::multimap parameters () override + { + using Map = std::multimap; + + Map base_param = Transport::parameters(); + Map new_param; + Map compsite_map; + + auto _const_type_name = ConstHydrodynamicDispersion::name; + + if constexpr (Traits::dim == 2) + new_param = Map({ + {_const_type_name+"_xx", _const_hydrodynamic_dispersion.value[0][0]}, + {_const_type_name+"_xy", _const_hydrodynamic_dispersion.value[1][0]}, + {_const_type_name+"_yx", _const_hydrodynamic_dispersion.value[0][1]}, + {_const_type_name+"_yy", _const_hydrodynamic_dispersion.value[1][1]} + }); + else if constexpr (Traits::dim == 3) + new_param = Map({ + {_const_type_name+"_xx", _const_hydrodynamic_dispersion.value[0][0]}, + {_const_type_name+"_xy", _const_hydrodynamic_dispersion.value[1][0]}, + {_const_type_name+"_xz", _const_hydrodynamic_dispersion.value[2][0]}, + {_const_type_name+"_yx", _const_hydrodynamic_dispersion.value[0][1]}, + {_const_type_name+"_yy", _const_hydrodynamic_dispersion.value[1][1]}, + {_const_type_name+"_yz", _const_hydrodynamic_dispersion.value[2][1]}, + {_const_type_name+"_zx", _const_hydrodynamic_dispersion.value[0][2]}, + {_const_type_name+"_xy", _const_hydrodynamic_dispersion.value[1][2]}, + {_const_type_name+"_zz", _const_hydrodynamic_dispersion.value[2][2]} + }); + else + static_assert(true, + "ConstHydrodynamicDispersion not available for this dimension"); + + compsite_map.insert(base_param.begin(),base_param.end()); + compsite_map.insert(new_param.begin(),new_param.end()); + + return compsite_map; + } + + /// Return a map of parameter names and values for manipulation + /** \return Map: Parameter name, parameter value in this object + */ + std::multimap parameters () const override + { + using Map = std::multimap; + + Map base_param = Transport::parameters(); + Map new_param; + Map compsite_map; + + auto _const_type_name = ConstHydrodynamicDispersion::name; + + if constexpr (Traits::dim == 2) + new_param = Map({ + {_const_type_name+"_xx", _const_hydrodynamic_dispersion.value[0][0]}, + {_const_type_name+"_xy", _const_hydrodynamic_dispersion.value[1][0]}, + {_const_type_name+"_yx", _const_hydrodynamic_dispersion.value[0][1]}, + {_const_type_name+"_yy", _const_hydrodynamic_dispersion.value[1][1]} + }); + else if constexpr (Traits::dim == 3) + new_param = Map({ + {_const_type_name+"_xx", _const_hydrodynamic_dispersion.value[0][0]}, + {_const_type_name+"_xy", _const_hydrodynamic_dispersion.value[1][0]}, + {_const_type_name+"_xz", _const_hydrodynamic_dispersion.value[2][0]}, + {_const_type_name+"_yx", _const_hydrodynamic_dispersion.value[0][1]}, + {_const_type_name+"_yy", _const_hydrodynamic_dispersion.value[1][1]}, + {_const_type_name+"_yz", _const_hydrodynamic_dispersion.value[2][1]}, + {_const_type_name+"_zx", _const_hydrodynamic_dispersion.value[0][2]}, + {_const_type_name+"_xy", _const_hydrodynamic_dispersion.value[1][2]}, + {_const_type_name+"_zz", _const_hydrodynamic_dispersion.value[2][2]} + }); + else + static_assert(Dune::AlwaysFalse::value, + "ConstHydrodynamicDispersion not available for this dimension"); + + compsite_map.insert(base_param.begin(),base_param.end()); + compsite_map.insert(new_param.begin(),new_param.end()); + + return compsite_map; + } + + /// Clone the plymorphic class + /** \return unique pointer to the cloned object + */ + std::unique_ptr> clone () const override + { + using ThisType = ConstHydrodynamicDispersionParam; + return std::make_unique(*this); + } +}; + +} // namespace Parameterization +} // namespace Dune +} // namespace Dorie + +#endif // DUNE_DORIE_PARAM_CONST_HYDRODYNAMIC_DISPERSION_HH \ No newline at end of file diff --git a/dune/dorie/model/transport/parameterization/hydrodynamic_dispersion/eff_diffusion/const.hh b/dune/dorie/model/transport/parameterization/hydrodynamic_dispersion/eff_diffusion/const.hh new file mode 100644 index 0000000000000000000000000000000000000000..10e161ce00b2605ddf983de172f07b3563d50391 --- /dev/null +++ b/dune/dorie/model/transport/parameterization/hydrodynamic_dispersion/eff_diffusion/const.hh @@ -0,0 +1,118 @@ +#ifndef DUNE_DORIE_PARAM_CONST_EFFECTIVE_DIFFUSION_HH +#define DUNE_DORIE_PARAM_CONST_EFFECTIVE_DIFFUSION_HH + +#include +#include +#include + +#include + +namespace Dune { +namespace Dorie { +namespace Parameterization { + + +/*-------------------------------------------------------------------------*//** + * @brief Class for constant effective diffusion. + * + * @ingroup TransportParam + * @author Santiago Ospina De Los Ríos + * @date 2019 + * + * @tparam Traits the base traits + */ +template +class ConstEffectiveDiffusionParam + : public EffectiveDiffusionParam +{ +private: + using RangeField = typename Traits::RF; + using Base = EffectiveDiffusionParam; + +public: + using EffectiveDiffusion = typename Base::EffectiveDiffusion; + using WaterFlux = typename Base::WaterFlux; + using WaterContent = typename Base::WaterContent; + + /// Parameter defining the constant effective diffusion + struct ConstEffectiveDiffusion + { + RangeField value; + inline static const std::string name = "eff_diff"; + }; + + /// The name of this parameterization type (This is not the instance name) + static inline std::string type = "Deff_const"; + + ConstEffectiveDiffusion _const_eff_diff; + +public: + + /// Construct with default-initialized parameters + /** \param name The name associated with this soil layer + */ + ConstEffectiveDiffusionParam (const std::string name) : + EffectiveDiffusionParam(name) + { } + + /// Construct from a tuple of parameters + /** \param name The name associated with this soil layer + * \param parameters Tuple of parameters to use in this parameterization + */ + template + ConstEffectiveDiffusionParam ( + const std::string name, + const std::tuple parameters) : + EffectiveDiffusionParam(name, parameters), + _const_eff_diff(std::get(parameters)) + { } + + /// Add default destructor to clarify override + ~ConstEffectiveDiffusionParam () override = default; + +public: + + /// Return the effective diffusion + /** {WaterFlux,WaterContent} -> EffectiveDiffusion + */ + std::function + effective_diffusion_f () const override + { + return [this](const WaterFlux, const WaterContent) { + return EffectiveDiffusion{_const_eff_diff.value}; + }; + } + + /// Return a map of parameter names and values for manipulation + /** \return Map: Parameter name, parameter value in this object + */ + std::multimap parameters () override + { + const std::string _const_eff_diff_name = ConstEffectiveDiffusion::name; + return {{_const_eff_diff_name, _const_eff_diff.value}}; + } + + /// Return a map of parameter names and values for manipulation + /** \return Map: Parameter name, parameter value in this object + */ + std::multimap parameters () const override + { + const std::string _const_eff_diff_name = ConstEffectiveDiffusion::name; + return {{_const_eff_diff_name, _const_eff_diff.value}}; + } + + /// Clone the plymorphic class + /** \return unique pointer to the cloned object + */ + std::unique_ptr> clone () const override + { + using ThisType = ConstEffectiveDiffusionParam; + return std::make_unique(*this); + } +}; + +} // namespace Parameterization +} // namespace Dune +} // namespace Dorie + +#endif // DUNE_DORIE_PARAM_CONST_EFFECTIVE_DIFFUSION_HH \ No newline at end of file diff --git a/dune/dorie/model/transport/parameterization/hydrodynamic_dispersion/eff_diffusion/interface.hh b/dune/dorie/model/transport/parameterization/hydrodynamic_dispersion/eff_diffusion/interface.hh new file mode 100644 index 0000000000000000000000000000000000000000..f76e0b885010b5a041d441af034dcf228c3f41a2 --- /dev/null +++ b/dune/dorie/model/transport/parameterization/hydrodynamic_dispersion/eff_diffusion/interface.hh @@ -0,0 +1,99 @@ +#ifndef DUNE_DORIE_PARAM_EFFECTIVE_DIFFUSION_INTERFACE_HH +#define DUNE_DORIE_PARAM_EFFECTIVE_DIFFUSION_INTERFACE_HH + +#include +#include +#include + +#include + +namespace Dune { +namespace Dorie { +namespace Parameterization { + +/*-------------------------------------------------------------------------*//** + * @brief Class for effective diffusion. + * + * @ingroup TransportParam + * @author Santiago Ospina De Los Ríos + * @date 2019 + * + * @tparam Traits the base traits + */ +template +class EffectiveDiffusionParam +{ +private: + using RangeField = typename Traits::RF; + using BaseP = Transport; + +public: + + /// Type of the effective diffusion + struct EffectiveDiffusion + { + RangeField value; + }; + + using WaterFlux = typename BaseP::WaterFlux; + using WaterContent = typename BaseP::WaterContent; + + //! The name of this parameterization instance, associated with the layer. + const std::string _name; + +public: + /// Construct with default-initialized parameters + /** \param name The name associated with this soil layer + */ + EffectiveDiffusionParam (const std::string name) + : _name(name) + { } + + /// Construct from a tuple of parameters + /** \param name The name associated with this soil layer + * \param parameters Tuple of parameters to use in this parameterization + */ + template + EffectiveDiffusionParam ( + const std::string name, + const std::tuple parameters) + : _name(name) + { } + + /// Default constructor (virtual). + virtual ~EffectiveDiffusionParam () = default; + + /// Return the name of this parameterization instance. + const std::string& get_name() const { return _name; } + + /// Return a bound version of the effective diffusion + /** \return Function: {Water Flux, Water Content} -> Effective Diff. + */ + virtual std::function< + EffectiveDiffusion(const WaterFlux, const WaterContent)> + effective_diffusion_f () const = 0; + + /// Return a map referecing all parameters by their names. + /** \return Map. Key: Name of parameter (string). + * Value: Value of parameter (double&) + */ + virtual std::multimap parameters () = 0; + + /// Return a map referecing all parameters by their names. + /** \return Map. Key: Name of parameter (string). + * Value: Value of parameter (const double&) + */ + virtual std::multimap parameters () const = 0; + + /// Clone the plymorphic class + /** \return unique pointer to the cloned object + */ + virtual std::unique_ptr> clone () const = 0; + +}; + +} // namespace Parameterization +} // namespace Dune +} // namespace Dorie + +#endif // DUNE_DORIE_PARAM_EFFECTIVE_DIFFUSION_INTERFACE_HH \ No newline at end of file diff --git a/dune/dorie/model/transport/parameterization/hydrodynamic_dispersion/eff_diffusion/millington_quirk_1.hh b/dune/dorie/model/transport/parameterization/hydrodynamic_dispersion/eff_diffusion/millington_quirk_1.hh new file mode 100644 index 0000000000000000000000000000000000000000..b17189ad61f702875166c8575ac972cc51a5f21c --- /dev/null +++ b/dune/dorie/model/transport/parameterization/hydrodynamic_dispersion/eff_diffusion/millington_quirk_1.hh @@ -0,0 +1,134 @@ +#ifndef DUNE_DORIE_PARAM_MILLIGTON_QUIRK_1_HH +#define DUNE_DORIE_PARAM_MILLIGTON_QUIRK_1_HH + +#include +#include +#include + +#include + +namespace Dune { +namespace Dorie { +namespace Parameterization { + +/*-------------------------------------------------------------------------*//** + * @brief Class for effective diffusion with Milington Quirk I. + * @details Implements the relation @f$ \mathsf{D}^\mathsf{eff} + * =\frac{\theta_w^{7/3}\mathsf{D}_\mathsf{m}}{\phi^2} @f$ + * where @f$ \theta_w @f$ is the water content, + * @f$ \mathsf{D}_\mathsf{m} @f$ the molecular diffusion, and + * @f$ \phi @f$ the porosity of the medium. + * + * @ingroup TransportParam + * @author Santiago Ospina De Los Ríos + * @date 2019 + * + * @tparam Traits the base traits + */ +template +class MillingtonQuirk1 + : public EffectiveDiffusionParam +{ +private: + using RangeField = typename Traits::RF; + using Base = EffectiveDiffusionParam; + +public: + using EffectiveDiffusion = typename Base::EffectiveDiffusion; + using WaterFlux = typename Base::WaterFlux; + using WaterContent = typename Base::WaterContent; + + /// Parameter defining molecular diffusion + struct ModelcularDiffusionType + { + RangeField value; + inline static const std::string name = "mol_diff"; + }; + + /// Parameter defining porosity + struct PorosityType + { + RangeField value; + inline static const std::string name = "phi"; + }; + + /// The name of this parameterization type (This is not the instance name) + static inline std::string type = "Deff_MQ1"; + + ModelcularDiffusionType _mol_diff; + PorosityType _porosity; + +public: + + /// Construct with default-initialized parameters + /** \param name The name associated with this soil layer + */ + MillingtonQuirk1 (const std::string name) : + Base(name) + { } + + /// Construct from a tuple of parameters + /** \param name The name associated with this soil layer + * \param parameters Tuple of parameters to use in this parameterization + */ + template + MillingtonQuirk1 ( + const std::string name, + const std::tuple parameters) : + Base(name, parameters), + _mol_diff(std::get(parameters)), + _porosity(std::get(parameters)) + { } + + /// Add default destructor to clarify override + ~MillingtonQuirk1 () override = default; + + /// Return the effective diffusion + /** {WaterFlux,WaterContent} -> EffectiveDiffusion + */ + std::function + effective_diffusion_f () const override + { + return [this](const WaterFlux water_flux, const WaterContent water_content) { + auto val = _mol_diff.value*std::pow(water_content.value,7./3.)/std::pow(_porosity.value,2./3.); + return EffectiveDiffusion{val}; + }; + } + + /// Return a map of parameter names and values for manipulation + /** \return Map: Parameter name, parameter value in this object + */ + std::multimap parameters () override + { + return { + {ModelcularDiffusionType::name, _mol_diff.value}, + {PorosityType::name, _porosity.value} + }; + } + + /// Return a map of parameter names and values for manipulation + /** \return Map: Parameter name, parameter value in this object + */ + std::multimap parameters () const override + { + return { + {ModelcularDiffusionType::name, _mol_diff.value}, + {PorosityType::name, _porosity.value} + }; + } + + /// Return a map of parameter names and values for manipulation + /** \return Map: Parameter name, parameter value in this object + */ + std::unique_ptr> clone () const override + { + using ThisType = MillingtonQuirk1; + return std::make_unique(*this); + } +}; + +} // namespace Parameterization +} // namespace Dune +} // namespace Dorie + +#endif // DUNE_DORIE_PARAM_MILLIGTON_QUIRK_1_HH \ No newline at end of file diff --git a/dune/dorie/model/transport/parameterization/hydrodynamic_dispersion/eff_diffusion/millington_quirk_2.hh b/dune/dorie/model/transport/parameterization/hydrodynamic_dispersion/eff_diffusion/millington_quirk_2.hh new file mode 100644 index 0000000000000000000000000000000000000000..8a822ca4279d42545681085e8bb3a9a67a5d9e78 --- /dev/null +++ b/dune/dorie/model/transport/parameterization/hydrodynamic_dispersion/eff_diffusion/millington_quirk_2.hh @@ -0,0 +1,134 @@ +#ifndef DUNE_DORIE_PARAM_MILLIGTON_QUIRK_2_HH +#define DUNE_DORIE_PARAM_MILLIGTON_QUIRK_2_HH + +#include +#include +#include + +#include + +namespace Dune { +namespace Dorie { +namespace Parameterization { + +/*-------------------------------------------------------------------------*//** + * @brief Class for effective diffusion with Milington Quirk II. + * @details Implements the relation @f$ \mathsf{D}^\mathsf{eff} + * =\frac{\theta_w\mathsf{D}_\mathsf{m}}{\phi^{2/3}} @f$ + * where @f$ \theta_w @f$ is the water content, + * @f$ \mathsf{D}_\mathsf{m} @f$ the molecular diffusion, and + * @f$ \phi @f$ the porosity of the medium. + * + * @ingroup TransportParam + * @author Santiago Ospina De Los Ríos + * @date 2019 + * + * @tparam Traits the base traits + */ +template +class MillingtonQuirk2 + : public EffectiveDiffusionParam +{ +private: + using RangeField = typename Traits::RF; + using Base = EffectiveDiffusionParam; + +public: + using EffectiveDiffusion = typename Base::EffectiveDiffusion; + using WaterFlux = typename Base::WaterFlux; + using WaterContent = typename Base::WaterContent; + + /// Parameter defining molecular diffusion + struct ModelcularDiffusionType + { + RangeField value; + inline static const std::string name = "mol_diff"; + }; + + /// Parameter defining porosity + struct PorosityType + { + RangeField value; + inline static const std::string name = "phi"; + }; + + /// The name of this parameterization type (This is not the instance name) + static inline std::string type = "Deff_MQ2"; + + ModelcularDiffusionType _mol_diff; + PorosityType _porosity; + +public: + + /// Construct with default-initialized parameters + /** \param name The name associated with this soil layer + */ + MillingtonQuirk2 (const std::string name) : + Base(name) + { } + + /// Construct from a tuple of parameters + /** \param name The name associated with this soil layer + * \param parameters Tuple of parameters to use in this parameterization + */ + template + MillingtonQuirk2 ( + const std::string name, + const std::tuple parameters) : + Base(name, parameters), + _mol_diff(std::get(parameters)), + _porosity(std::get(parameters)) + { } + + /// Add default destructor to clarify override + ~MillingtonQuirk2 () override = default; + + /// Return the effective diffusion + /** {WaterFlux,WaterContent} -> EffectiveDiffusion + */ + std::function + effective_diffusion_f () const override + { + return [this](const WaterFlux water_flux, const WaterContent water_content) { + auto val = _mol_diff.value*water_content.value/std::pow(_porosity.value,2./3.); + return EffectiveDiffusion{val}; + }; + } + + /// Return a map of parameter names and values for manipulation + /** \return Map: Parameter name, parameter value in this object + */ + std::multimap parameters () override + { + return { + {ModelcularDiffusionType::name, _mol_diff.value}, + {PorosityType::name, _porosity.value} + }; + } + + /// Return a map of parameter names and values for manipulation + /** \return Map: Parameter name, parameter value in this object + */ + std::multimap parameters () const override + { + return { + {ModelcularDiffusionType::name, _mol_diff.value}, + {PorosityType::name, _porosity.value} + }; + } + + /// Return a map of parameter names and values for manipulation + /** \return Map: Parameter name, parameter value in this object + */ + std::unique_ptr> clone () const override + { + using ThisType = MillingtonQuirk2; + return std::make_unique(*this); + } +}; + +} // namespace Parameterization +} // namespace Dune +} // namespace Dorie + +#endif // DUNE_DORIE_PARAM_MILLIGTON_QUIRK_2_HH \ No newline at end of file diff --git a/dune/dorie/model/transport/parameterization/hydrodynamic_dispersion/eff_hydromechanic_dispersion/const.hh b/dune/dorie/model/transport/parameterization/hydrodynamic_dispersion/eff_hydromechanic_dispersion/const.hh new file mode 100644 index 0000000000000000000000000000000000000000..25bf7d1cd6615a3b9d809bde78cdb0a85d76cd4b --- /dev/null +++ b/dune/dorie/model/transport/parameterization/hydrodynamic_dispersion/eff_hydromechanic_dispersion/const.hh @@ -0,0 +1,157 @@ +#ifndef DUNE_DORIE_PARAM_CONST_EFF_HM_DISP_HH +#define DUNE_DORIE_PARAM_CONST_EFF_HM_DISP_HH + +#include +#include +#include + +#include + +namespace Dune { +namespace Dorie { +namespace Parameterization { + +/*-------------------------------------------------------------------------*//** + * @brief Class for constant effective hydromechanic dispersion. + * + * @ingroup TransportParam + * @author Santiago Ospina De Los Ríos + * @date 2019 + * + * @tparam Traits the base traits + */ +template +class ConstEffectiveHydromechanicDispersionParam : + public EffectiveHydromechanicDispersionParam +{ +private: + using Tensor = typename Traits::Tensor; + using Base = EffectiveHydromechanicDispersionParam; + +public: + using EffectiveHydromechanicDispersion = typename Base::EffectiveHydromechanicDispersion; + using WaterFlux = typename Base::WaterFlux; + using WaterContent = typename Base::WaterContent; + + /// Parameter defining the constant hydromechanic dispersion tensor + struct ConstEffectiveHydromechanicDispersion + { + Tensor value; + inline static const std::string name = "eff_hydromechanic_disp"; + }; + + /// The name of this parameterization type (This is not the instance name) + static inline std::string type = "Dhm_const"; + + ConstEffectiveHydromechanicDispersion _const_eff_hydromechanic_dispersion; + +public: + + /// Construct with default-initialized parameters + /** \param name The name associated with this soil layer + */ + ConstEffectiveHydromechanicDispersionParam (const std::string name) : + EffectiveHydromechanicDispersionParam(name) + { } + + /// Construct from a tuple of parameters + /** \param name The name associated with this soil layer + * \param parameters Tuple of parameters to use in this parameterization + */ + template + ConstEffectiveHydromechanicDispersionParam ( + const std::string name, + const std::tuple parameters) : + Base(name, parameters), + _const_eff_hydromechanic_dispersion(std::get(parameters)) + { } + + /// Add default destructor to clarify override + ~ConstEffectiveHydromechanicDispersionParam () override = default; + + /// Return the hydromechanic dispersion + /** {WaterFlux,WaterContent} -> EffectiveHydromechanicDispersion + */ + std::function + effective_hydromechanic_dispersion_f () const override + { + return [this](const WaterFlux, const WaterContent) { + return EffectiveHydromechanicDispersion{_const_eff_hydromechanic_dispersion.value}; + }; + } + + /// Return a map of parameter names and values for manipulation + /** \return Map: Parameter name, parameter value in this object + */ + std::multimap parameters () override + { + const std::string _const_eff_hydromechanic_dispersion_name = ConstEffectiveHydromechanicDispersion::name; + if constexpr (Traits::dim == 2) + return { + {_const_eff_hydromechanic_dispersion_name+"_xx", _const_eff_hydromechanic_dispersion.value[0][0]}, + {_const_eff_hydromechanic_dispersion_name+"_xy", _const_eff_hydromechanic_dispersion.value[0][1]}, + {_const_eff_hydromechanic_dispersion_name+"_yx", _const_eff_hydromechanic_dispersion.value[1][0]}, + {_const_eff_hydromechanic_dispersion_name+"_yy", _const_eff_hydromechanic_dispersion.value[1][1]} + }; + else if constexpr (Traits::dim == 3) + return { + {_const_eff_hydromechanic_dispersion_name+"_xx", _const_eff_hydromechanic_dispersion.value[0][0]}, + {_const_eff_hydromechanic_dispersion_name+"_xy", _const_eff_hydromechanic_dispersion.value[0][1]}, + {_const_eff_hydromechanic_dispersion_name+"_xz", _const_eff_hydromechanic_dispersion.value[0][2]}, + {_const_eff_hydromechanic_dispersion_name+"_yx", _const_eff_hydromechanic_dispersion.value[1][0]}, + {_const_eff_hydromechanic_dispersion_name+"_yy", _const_eff_hydromechanic_dispersion.value[1][1]}, + {_const_eff_hydromechanic_dispersion_name+"_yz", _const_eff_hydromechanic_dispersion.value[1][2]}, + {_const_eff_hydromechanic_dispersion_name+"_zx", _const_eff_hydromechanic_dispersion.value[2][0]}, + {_const_eff_hydromechanic_dispersion_name+"_xy", _const_eff_hydromechanic_dispersion.value[2][1]}, + {_const_eff_hydromechanic_dispersion_name+"_zz", _const_eff_hydromechanic_dispersion.value[2][2]} + }; + else + static_assert(true, + "ConstEffectiveHydromechanicDispersion not available for this dimension"); + } + + /// Return a map of parameter names and values for manipulation + /** \return Map: Parameter name, parameter value in this object + */ + std::multimap parameters () const override + { + const std::string _const_eff_hydromechanic_dispersion_name = ConstEffectiveHydromechanicDispersion::name; + if constexpr (Traits::dim == 2) + return { + {_const_eff_hydromechanic_dispersion_name+"_xx", _const_eff_hydromechanic_dispersion.value[0][0]}, + {_const_eff_hydromechanic_dispersion_name+"_xy", _const_eff_hydromechanic_dispersion.value[0][1]}, + {_const_eff_hydromechanic_dispersion_name+"_yx", _const_eff_hydromechanic_dispersion.value[1][0]}, + {_const_eff_hydromechanic_dispersion_name+"_yy", _const_eff_hydromechanic_dispersion.value[1][1]} + }; + else if constexpr (Traits::dim == 3) + return { + {_const_eff_hydromechanic_dispersion_name+"_xx", _const_eff_hydromechanic_dispersion.value[0][0]}, + {_const_eff_hydromechanic_dispersion_name+"_xy", _const_eff_hydromechanic_dispersion.value[0][1]}, + {_const_eff_hydromechanic_dispersion_name+"_xz", _const_eff_hydromechanic_dispersion.value[0][2]}, + {_const_eff_hydromechanic_dispersion_name+"_yx", _const_eff_hydromechanic_dispersion.value[1][0]}, + {_const_eff_hydromechanic_dispersion_name+"_yy", _const_eff_hydromechanic_dispersion.value[1][1]}, + {_const_eff_hydromechanic_dispersion_name+"_yz", _const_eff_hydromechanic_dispersion.value[1][2]}, + {_const_eff_hydromechanic_dispersion_name+"_zx", _const_eff_hydromechanic_dispersion.value[2][0]}, + {_const_eff_hydromechanic_dispersion_name+"_xy", _const_eff_hydromechanic_dispersion.value[2][1]}, + {_const_eff_hydromechanic_dispersion_name+"_zz", _const_eff_hydromechanic_dispersion.value[2][2]} + }; + else + static_assert(Dune::AlwaysFalse::value, + "ConstEffectiveHydromechanicDispersion not available for this dimension"); + } + + /// Return a map of parameter names and values for manipulation + /** \return Map: Parameter name, parameter value in this object + */ + std::unique_ptr> clone () const override + { + using ThisType = ConstEffectiveHydromechanicDispersionParam; + return std::make_unique(*this); + } +}; + +} // namespace Parameterization +} // namespace Dune +} // namespaie + +#endif // DUNE_DORIE_PARAM_CONST_EFF_HM_DISPERSION_HH \ No newline at end of file diff --git a/dune/dorie/model/transport/parameterization/hydrodynamic_dispersion/eff_hydromechanic_dispersion/interface.hh b/dune/dorie/model/transport/parameterization/hydrodynamic_dispersion/eff_hydromechanic_dispersion/interface.hh new file mode 100644 index 0000000000000000000000000000000000000000..e6c3e0aedafeb82eb139fb4e0b96967441f2ccdb --- /dev/null +++ b/dune/dorie/model/transport/parameterization/hydrodynamic_dispersion/eff_hydromechanic_dispersion/interface.hh @@ -0,0 +1,100 @@ +#ifndef DUNE_DORIE_PARAM_EFFECTIVE_HYDROMECHANIC_DISPERSION_INTERFACE_HH +#define DUNE_DORIE_PARAM_EFFECTIVE_HYDROMECHANIC_DISPERSION_INTERFACE_HH + +#include +#include +#include + +#include + +namespace Dune { +namespace Dorie { +namespace Parameterization { + +/*-------------------------------------------------------------------------*//** + * @brief Class for effective hydromechanic dispersion. + * + * @ingroup TransportParam + * @author Santiago Ospina De Los Ríos + * @date 2019 + * + * @tparam Traits the base traits + */ +template +class EffectiveHydromechanicDispersionParam +{ +private: + using BaseP = Transport; + using RangeField = typename Traits::RF; + using Scalar = typename Traits::Scalar; + using Vector = typename Traits::Vector; + using Tensor = typename Traits::Tensor; + +public: + + /// Type of the hydrodynamic dispersion + struct EffectiveHydromechanicDispersion + { + Tensor value; + }; + + using WaterFlux = typename BaseP::WaterFlux; + using WaterContent = typename BaseP::WaterContent; + + //! The name of this parameterization instance, associated with the layer. + const std::string _name; + + /// Construct with default-initialized parameters + /** \param name The name associated with this soil layer + */ + EffectiveHydromechanicDispersionParam (const std::string name) + : _name(name) + { } + + /// Construct from a tuple of parameters + /** \param name The name associated with this soil layer + * \param parameters Tuple of parameters to use in this parameterization + */ + template + EffectiveHydromechanicDispersionParam ( + const std::string name, + const std::tuple parameters) + : _name(name) + { } + + /// Default constructor (virtual). + virtual ~EffectiveHydromechanicDispersionParam () = default; + + /// Return the name of this parameterization instance. + const std::string& get_name() const { return _name; } + + /// Return a bound version of the hydromechanic dispersion tensor + /** \return Function: {Water Flux, Water Content} -> Hydromechanic Disp. Tensor + */ + virtual std::function< + EffectiveHydromechanicDispersion(const WaterFlux, const WaterContent)> + effective_hydromechanic_dispersion_f () const = 0; + + /// Return a map referecing all parameters by their names. + /** \return Map. Key: Name of parameter (string). + * Value: Value of parameter (double&) + */ + virtual std::multimap parameters () = 0; + + /// Return a map referecing all parameters by their names. + /** \return Map. Key: Name of parameter (string). + * Value: Value of parameter (const double&) + */ + virtual std::multimap parameters () const = 0; + + /// Clone the plymorphic class + /** \return unique pointer to the cloned object + */ + virtual std::unique_ptr> clone () const = 0; +}; + +} // namespace Parameterization +} // namespace Dune +} // namespace Dorie + +#endif // DUNE_DORIE_PARAM_EFFECTIVE_HYDROMECHANIC_DISPERSION_INTERFACE_HH \ No newline at end of file diff --git a/dune/dorie/model/transport/parameterization/hydrodynamic_dispersion/eff_hydromechanic_dispersion/isotropic.hh b/dune/dorie/model/transport/parameterization/hydrodynamic_dispersion/eff_hydromechanic_dispersion/isotropic.hh new file mode 100644 index 0000000000000000000000000000000000000000..1528a2f61a18189aa3caf760059f3b5e7cd59ca1 --- /dev/null +++ b/dune/dorie/model/transport/parameterization/hydrodynamic_dispersion/eff_hydromechanic_dispersion/isotropic.hh @@ -0,0 +1,156 @@ +#ifndef DUNE_DORIE_PARAM_ISOTROPIC_EFFECTIVE_HYDROMECHANIC_DISPERSION_HH +#define DUNE_DORIE_PARAM_ISOTROPIC_EFFECTIVE_HYDROMECHANIC_DISPERSION_HH + +#include +#include +#include + +#include +#include +#include + +namespace Dune { +namespace Dorie { +namespace Parameterization { + +/*-------------------------------------------------------------------------*//** + * @brief Class for isotropic effective hydromechanic dispersion. + * @details Implements an interface for the the relation + * @f$ {\mathsf{D}_\mathsf{hm}^\mathsf{eff}}_{ij} = + * [\lamda_l-\lambda_t] \frac{v_iv_j}{|v|}+\lambda_t|v|\delta_ij + * where @f$ v:=j_w/\theta @f$ is the water velocity, + * @f$ \delta_ij @f$ the kroneker delta, @f$ \lambda_t @f$ the + * transverse dispersivity, and @f$ \lambda_l @f$ the longitudinal + * dispersivity + * + * @ingroup TransportParam + * @author Santiago Ospina De Los Ríos + * @date 2019 + * + * @tparam Traits the base traits + */ +template +class IsotropicEffectiveHydromechanicDispersion + : public EffectiveHydromechanicDispersionParam +{ +private: + using RangeField = typename Traits::RF; + using Tensor = typename Traits::Tensor; + using Base = EffectiveHydromechanicDispersionParam; + using LambdaT = TransverseDispersivityParam; + using LambdaL = LongitudinalDispersivityParam; + +public: + using EffectiveHydromechanicDispersion = typename Base::EffectiveHydromechanicDispersion; + using WaterFlux = typename Base::WaterFlux; + using WaterContent = typename Base::WaterContent; + + /// Parameter defining the longinutinal dispersivity + struct LongitudinalDispersivity + { + RangeField value; + }; + + /// Parameter defining the transverse dispersivity + struct TransverseDispersivity + { + RangeField value; + }; + + /// The name of this parameterization type (This is not the instance name) + static inline std::string type = "Dhm_iso"; + + std::shared_ptr _long_disp; + std::shared_ptr _trans_disp; + +public: + + /// Construct with polymorphic objects + /** \param name The name associated with this soil layer + * \param long_disp The class for the longitudinal dispersivity + * \param trans_disp The class for the transverse dispersivity + */ + template + IsotropicEffectiveHydromechanicDispersion ( + const std::string name, + std::shared_ptr long_disp, + std::shared_ptr trans_disp) : + Base(name), + _long_disp(long_disp), + _trans_disp(trans_disp) + { } + + /// Add default destructor to clarify override + ~IsotropicEffectiveHydromechanicDispersion () override = default; + + /// Return the effective hydromechanic dispersion + /** {WaterFlux,WaterContent} -> EffectiveHydromechanicDispersion + */ + std::function + effective_hydromechanic_dispersion_f () const override + { + return [this](const WaterFlux water_flux, const WaterContent water_content) { + const auto lambda_l_f = _long_disp->longitudinal_dispersivity_f(); + const LongitudinalDispersivity lambda_l = lambda_l_f(water_flux,water_content); + const auto lambda_t_f = _trans_disp->transverse_dispersivity_f(); + const TransverseDispersivity lambda_t = lambda_t_f(water_flux,water_content); + Tensor D; + auto difference = (lambda_l.value-lambda_t.value); + auto v_abs = water_flux.value.two_norm(); + for (int i = 0; i < Traits::dim; ++i) + for (int j = 0; j < Traits::dim; ++j) { + D[i][j] = 0.; + if (Dune::FloatCmp::gt(v_abs,0.)) { + D[i][j] += difference*water_flux.value[i]*water_flux.value[j]/(v_abs*water_content.value); // !!! + if (i==j) + D[i][j] += lambda_t.value * v_abs/water_content.value; + } + } + return EffectiveHydromechanicDispersion{D}; + }; + } + + /// Return a map of parameter names and values for manipulation + /** \return Map: Parameter name, parameter value in this object + */ + virtual std::multimap parameters () override + { + auto lambda_l_param = _long_disp->parameters(); + auto lambda_t_param = _trans_disp->parameters(); + + std::multimap compsite_map(lambda_l_param); + compsite_map.insert(lambda_t_param.begin(),lambda_t_param.end()); + + return compsite_map; + } + + /// Return a map of parameter names and values for manipulation + /** \return Map: Parameter name, parameter value in this object + */ + virtual std::multimap parameters () const override + { + auto lambda_l_param = _long_disp->parameters(); + auto lambda_t_param = _trans_disp->parameters(); + + std::multimap compsite_map; + compsite_map.insert(lambda_l_param.begin(),lambda_l_param.end()); + compsite_map.insert(lambda_t_param.begin(),lambda_t_param.end()); + + return compsite_map; + } + + /// Return a map of parameter names and values for manipulation + /** \return Map: Parameter name, parameter value in this object + */ + std::unique_ptr> clone () const override + { + using ThisType = IsotropicEffectiveHydromechanicDispersion; + return std::make_unique(this->get_name(),_long_disp->clone(),_trans_disp->clone()); + } +}; + +} // namespace Parameterization +} // namespace Dune +} // namespace Dorie + +#endif // DUNE_DORIE_PARAM_ISOTROPIC_EFFECTIVE_HYDROMECHANIC_DISPERSION_HH \ No newline at end of file diff --git a/dune/dorie/model/transport/parameterization/hydrodynamic_dispersion/eff_hydromechanic_dispersion/longitudinal/const.hh b/dune/dorie/model/transport/parameterization/hydrodynamic_dispersion/eff_hydromechanic_dispersion/longitudinal/const.hh new file mode 100644 index 0000000000000000000000000000000000000000..4c4d001c4571da59bd0a27b9b44d578dddf18143 --- /dev/null +++ b/dune/dorie/model/transport/parameterization/hydrodynamic_dispersion/eff_hydromechanic_dispersion/longitudinal/const.hh @@ -0,0 +1,117 @@ +#ifndef DUNE_DORIE_PARAM_CONST_LONGITUDINAL_DISPERSIVITY_INTERFACE_HH +#define DUNE_DORIE_PARAM_CONST_LONGITUDINAL_DISPERSIVITY_INTERFACE_HH + +#include +#include +#include + +#include + +namespace Dune { +namespace Dorie { +namespace Parameterization { + +/*-------------------------------------------------------------------------*//** + * @brief Class for constant longitudinal dispersivity. + * + * @ingroup TransportParam + * @author Santiago Ospina De Los Ríos + * @date 2019 + * + * @tparam Traits the base traits + */ +template +class ConstLongitudinalDispersivityParam + : public LongitudinalDispersivityParam +{ +private: + using Base = LongitudinalDispersivityParam; + using RangeField = typename Traits::RF; + +public: + using WaterFlux = typename Base::WaterFlux; + using WaterContent = typename Base::WaterContent; + using LongitudinalDispersivity = typename Base::LongitudinalDispersivity; + + /// Parameter defining the constant longitudinal dispersivity + struct ConstLongitudinalDispersivity + { + RangeField value; + inline static const std::string name = "lambda_l"; + }; + + /// The name of this parameterization type (This is not the instance name) + static inline std::string type = "const"; + + ConstLongitudinalDispersivity _lambda_l; + +public: + + /// Construct with default-initialized parameters + /** \param name The name associated with this soil layer + */ + ConstLongitudinalDispersivityParam (const std::string name) + : Base(name) + { } + + /// Construct from a tuple of parameters + /** \param name The name associated with this soil layer + * \param parameters Tuple of parameters to use in this parameterization + */ + template + ConstLongitudinalDispersivityParam ( + const std::string name, + const std::tuple parameters) : + Base(name, parameters), + _lambda_l(std::get(parameters)) + { } + + /// Add default destructor to clarify override + ~ConstLongitudinalDispersivityParam () override = default; + + /// Return the longitudinal dispersivity + /** {WaterFlux,WaterContent} -> LongitudinalDispersivity + */ + std::function< + LongitudinalDispersivity(const WaterFlux, const WaterContent)> + longitudinal_dispersivity_f () const override + { + return [this](const WaterFlux, const WaterContent) { + return LongitudinalDispersivity{_lambda_l.value}; + }; + } + + /// Return a map of parameter names and values for manipulation + /** \return Map: Parameter name, parameter value in this object + */ + std::multimap parameters () override + { + return {{ConstLongitudinalDispersivity::name,_lambda_l.value}}; + } + + /// Return a map of parameter names and values for manipulation + /** \return Map: Parameter name, parameter value in this object + */ + std::multimap parameters () const override + { + return {{ConstLongitudinalDispersivity::name,_lambda_l.value}}; + } + + /// Return a map of parameter names and values for manipulation + /** \return Map: Parameter name, parameter value in this object + */ + std::unique_ptr> clone () const override + { + using ThisType = ConstLongitudinalDispersivityParam; + return std::make_unique(*this); + } + +private: + const std::string _name; +}; + +} // namespace Parameterizatoion +} // namespace Dune +} // namespace Dorie + +#endif // DUNE_DORIE_PARAM_CONST_LONGITUDINAL_DISPERSIVITY_INTERFACE_HH \ No newline at end of file diff --git a/dune/dorie/model/transport/parameterization/hydrodynamic_dispersion/eff_hydromechanic_dispersion/longitudinal/interface.hh b/dune/dorie/model/transport/parameterization/hydrodynamic_dispersion/eff_hydromechanic_dispersion/longitudinal/interface.hh new file mode 100644 index 0000000000000000000000000000000000000000..bdeac140470566517195a946ef5dfb37ab2cc5d3 --- /dev/null +++ b/dune/dorie/model/transport/parameterization/hydrodynamic_dispersion/eff_hydromechanic_dispersion/longitudinal/interface.hh @@ -0,0 +1,92 @@ +#ifndef DUNE_DORIE_PARAM_LONGITUDINAL_DISPERSIVITY_INTERFACE_HH +#define DUNE_DORIE_PARAM_LONGITUDINAL_DISPERSIVITY_INTERFACE_HH + +#include +#include +#include + +namespace Dune { +namespace Dorie { +namespace Parameterization { + +template +class IsotropicEffectiveHydromechanicDispersion; + +/*-------------------------------------------------------------------------*//** + * @brief Class for longitudinal dispersivity. + * + * @ingroup TransportParam + * @author Santiago Ospina De Los Ríos + * @date 2019 + * + * @tparam Traits the base traits + */ +template +class LongitudinalDispersivityParam +{ +private: + using BaseP = IsotropicEffectiveHydromechanicDispersion; + +public: + + using WaterFlux = typename BaseP::WaterFlux; + using WaterContent = typename BaseP::WaterContent; + using LongitudinalDispersivity = typename BaseP::LongitudinalDispersivity; + + //! The name of this parameterization instance, associated with the layer. + const std::string _name; + + /// Construct with default-initialized parameters + /** \param name The name associated with this soil layer + */ + LongitudinalDispersivityParam (const std::string name) + : _name(name) + { } + + /// Construct from a tuple of parameters + /** \param name The name associated with this soil layer + * \param parameters Tuple of parameters to use in this parameterization + */ + template + LongitudinalDispersivityParam ( + const std::string name, + const std::tuple parameters) + : _name(name) + { } + + /// Default constructor (virtual). + virtual ~LongitudinalDispersivityParam () = default; + + /// Return the name of this parameterization instance. + const std::string& get_name() const { return _name; } + + /// Return a bound version of the longitudinal dispersivity + /** \return Function: {Water Flux, Water Content} -> Longitudinal Dispersivity + */ + virtual std::function< + LongitudinalDispersivity(const WaterFlux, const WaterContent)> + longitudinal_dispersivity_f () const = 0; + + /// Return a map referecing all parameters by their names. + /** \return Map. Key: Name of parameter (string). + * Value: Value of parameter (double&) + */ + virtual std::multimap parameters () = 0; + + /// Return a map referecing all parameters by their names. + /** \return Map. Key: Name of parameter (string). + * Value: Value of parameter (const double&) + */ + virtual std::multimap parameters () const = 0; + + /// Clone the plymorphic class + /** \return unique pointer to the cloned object + */ + virtual std::unique_ptr> clone () const = 0; +}; + +} // namespace Parameterization +} // namespace Dune +} // namespace Dorie + +#endif // DUNE_DORIE_PARAM_LONGITUDINAL_DISPERSIVITY_INTERFACE_HH \ No newline at end of file diff --git a/dune/dorie/model/transport/parameterization/hydrodynamic_dispersion/eff_hydromechanic_dispersion/transverse/const.hh b/dune/dorie/model/transport/parameterization/hydrodynamic_dispersion/eff_hydromechanic_dispersion/transverse/const.hh new file mode 100644 index 0000000000000000000000000000000000000000..db2d673c0f5e0bd7cc3b59d834e8c51beac5351f --- /dev/null +++ b/dune/dorie/model/transport/parameterization/hydrodynamic_dispersion/eff_hydromechanic_dispersion/transverse/const.hh @@ -0,0 +1,117 @@ +#ifndef DUNE_DORIE_PARAM_CONST_TRANSVERSE_DISPERSIVITY_INTERFACE_HH +#define DUNE_DORIE_PARAM_CONST_TRANSVERSE_DISPERSIVITY_INTERFACE_HH + +#include +#include +#include + +#include + +namespace Dune { +namespace Dorie { +namespace Parameterization { + +/*-------------------------------------------------------------------------*//** + * @brief Class for constant transverse dispersivity. + * + * @ingroup TransportParam + * @author Santiago Ospina De Los Ríos + * @date 2019 + * + * @tparam Traits the base traits + */ +template +class ConstTransverseDispersivityParam + : public TransverseDispersivityParam +{ +private: + using Base = TransverseDispersivityParam; + using RangeField = typename Traits::RF; + +public: + using WaterFlux = typename Base::WaterFlux; + using WaterContent = typename Base::WaterContent; + using TransverseDispersivity = typename Base::TransverseDispersivity; + + /// Parameter defining the constant transverse dispersivity + struct ConstTransverseDispersivity + { + RangeField value; + inline static const std::string name = "lambda_t"; + }; + + /// The name of this parameterization type (This is not the instance name) + static inline std::string type = "const"; + + ConstTransverseDispersivity _lambda_t; + +public: + + /// Construct with default-initialized parameters + /** \param name The name associated with this soil layer + */ + ConstTransverseDispersivityParam (const std::string name) + : Base(name) + { } + + /// Construct from a tuple of parameters + /** \param name The name associated with this soil layer + * \param parameters Tuple of parameters to use in this parameterization + */ + template + ConstTransverseDispersivityParam ( + const std::string name, + const std::tuple parameters) : + Base(name, parameters), + _lambda_t(std::get(parameters)) + { } + + /// Add default destructor to clarify override + ~ConstTransverseDispersivityParam () override = default; + + /// Return the transverse dispersivity + /** {WaterFlux,WaterContent} -> TransverseDispersivity + */ + std::function< + TransverseDispersivity(const WaterFlux, const WaterContent)> + transverse_dispersivity_f () const override + { + return [this](const WaterFlux, const WaterContent) { + return TransverseDispersivity{_lambda_t.value}; + }; + } + + /// Return a map of parameter names and values for manipulation + /** \return Map: Parameter name, parameter value in this object + */ + std::multimap parameters () override + { + return {{ConstTransverseDispersivity::name,_lambda_t.value}}; + } + + /// Return a map of parameter names and values for manipulation + /** \return Map: Parameter name, parameter value in this object + */ + std::multimap parameters () const override + { + return {{ConstTransverseDispersivity::name,_lambda_t.value}}; + } + + /// Return a map of parameter names and values for manipulation + /** \return Map: Parameter name, parameter value in this object + */ + std::unique_ptr> clone () const override + { + using ThisType = ConstTransverseDispersivityParam; + return std::make_unique(*this); + } + +private: + const std::string _name; +}; + +} // namespace Parameterization +} // namespace Dune +} // namespace Dorie + +#endif // DUNE_DORIE_PARAM_CONST_TRANSVERSE_DISPERSIVITY_INTERFACE_HH \ No newline at end of file diff --git a/dune/dorie/model/transport/parameterization/hydrodynamic_dispersion/eff_hydromechanic_dispersion/transverse/interface.hh b/dune/dorie/model/transport/parameterization/hydrodynamic_dispersion/eff_hydromechanic_dispersion/transverse/interface.hh new file mode 100644 index 0000000000000000000000000000000000000000..bc88bb15b941bce4b372671ac35632c8c0367729 --- /dev/null +++ b/dune/dorie/model/transport/parameterization/hydrodynamic_dispersion/eff_hydromechanic_dispersion/transverse/interface.hh @@ -0,0 +1,92 @@ +#ifndef DUNE_DORIE_PARAM_TRANSVERSE_DISPERSIVITY_INTERFACE_HH +#define DUNE_DORIE_PARAM_TRANSVERSE_DISPERSIVITY_INTERFACE_HH + +#include +#include +#include + +namespace Dune { +namespace Dorie { +namespace Parameterization { + +template +class IsotropicEffectiveHydromechanicDispersion; + +/*-------------------------------------------------------------------------*//** + * @brief Class for transverse dispersivity. + * + * @ingroup TransportParam + * @author Santiago Ospina De Los Ríos + * @date 2019 + * + * @tparam Traits the base traits + */ +template +class TransverseDispersivityParam +{ +private: + using BaseP = IsotropicEffectiveHydromechanicDispersion; + +public: + + using WaterFlux = typename BaseP::WaterFlux; + using WaterContent = typename BaseP::WaterContent; + using TransverseDispersivity = typename BaseP::TransverseDispersivity; + + //! The name of this parameterization instance, associated with the layer. + const std::string _name; + + /// Construct with default-initialized parameters + /** \param name The name associated with this soil layer + */ + TransverseDispersivityParam (const std::string name) + : _name(name) + { } + + /// Construct from a tuple of parameters + /** \param name The name associated with this soil layer + * \param parameters Tuple of parameters to use in this parameterization + */ + template + TransverseDispersivityParam ( + const std::string name, + const std::tuple parameters) + : _name(name) + { } + + /// Default constructor (virtual). + virtual ~TransverseDispersivityParam () = default; + + /// Return the name of this parameterization instance. + const std::string& get_name() const { return _name; } + + /// Return a bound version of the transverse dispersivity + /** \return Function: {Water Flux, Water Content} -> Transverse Dispersivity + */ + virtual std::function< + TransverseDispersivity(const WaterFlux, const WaterContent)> + transverse_dispersivity_f () const = 0; + + /// Return a map referecing all parameters by their names. + /** \return Map. Key: Name of parameter (string). + * Value: Value of parameter (double&) + */ + virtual std::multimap parameters () = 0; + + /// Return a map referecing all parameters by their names. + /** \return Map. Key: Name of parameter (string). + * Value: Value of parameter (const double&) + */ + virtual std::multimap parameters () const = 0; + + /// Return a clone of this object + /** \return a unique pointer with a copy of this object. + */ + virtual std::unique_ptr> clone () const = 0; +}; + +} // namespace Parameterization +} // namespace Dune +} // namespace Dorie + +#endif // DUNE_DORIE_PARAM_TRANSVERSE_DISPERSIVITY_INTERFACE_HH \ No newline at end of file diff --git a/dune/dorie/model/transport/parameterization/hydrodynamic_dispersion/power_law.hh b/dune/dorie/model/transport/parameterization/hydrodynamic_dispersion/power_law.hh new file mode 100644 index 0000000000000000000000000000000000000000..bc0a0ab0c8223e46ada21a0aee0a55ca8f7c2bc8 --- /dev/null +++ b/dune/dorie/model/transport/parameterization/hydrodynamic_dispersion/power_law.hh @@ -0,0 +1,155 @@ +#ifndef DUNE_DORIE_PARAM_POWER_LAW_HD_DISP_HH +#define DUNE_DORIE_PARAM_POWER_LAW_HD_DISP_HH + +#include +#include +#include + +#include + +namespace Dune { +namespace Dorie { +namespace Parameterization { + +/*-------------------------------------------------------------------------*//** + * @brief Class for power law hydrodynamic dispersion. + * @details Implements the relation @f$ \mathsf{D}_\mathsf{hd} + * =I \mathsf{D}_\mathsf{m}\gamma\mathsf{pe}^\alpha @f$ + * where @f$ I @f$ is the identity matrix, + * @f$ \mathsf{D}_\mathsf{m} @f$ the molecular diffusion, and + * @f$ \alpha @f$ and @f$ \gamma @f$ parameters of the power law. + * + * @ingroup TransportParam + * @author Santiago Ospina De Los Ríos + * @date 2019 + * + * @tparam Traits the base traits + */ +template +class PowerLawDispersion : + public Transport +{ +private: + using RangeField = typename Traits::RF; + using Base = Transport; + +public: + + using HydrodynamicDispersion = typename Base::HydrodynamicDispersion; + using WaterFlux = typename Base::WaterFlux; + using WaterContent = typename Base::WaterContent; + using Peclet = typename Transport::Peclet; + + /// Parameter defining the gamma constant + struct Gamma + { + RangeField value; + inline static const std::string name = "gamma"; + }; + + /// Parameter defining the alpha constant + struct Alpha + { + RangeField value; + inline static const std::string name = "alpha"; + }; + + /// The name of this parameterization type (This is not the instance name) + static inline std::string type = "Dhd_pl"; + + Gamma _gamma; + Alpha _alpha; + +public: + + /// Construct with default-initialized parameters + /** \param name The name associated with this soil layer + */ + PowerLawDispersion (const std::string name) : + Base(name) + { } + + /// Construct from a tuple of parameters + /** \param name The name associated with this soil layer + * \param parameters Tuple of parameters to use in this parameterization + */ + template + PowerLawDispersion ( + const std::string name, + const std::tuple parameters) : + Base(name, parameters), + _gamma(std::get(parameters)), + _alpha(std::get(parameters)) + { } + + /// Add default destructor to clarify override + ~PowerLawDispersion () override = default; + + /// Return the hydrodynamic dispersion + /** {WaterFlux,WaterContent} -> HydrodynamicDispersion + */ + std::function + hydrodynamic_dispersion_f () const override + { + return [this](const WaterFlux water_flux, const WaterContent water_content) { + auto diff = this->_mol_diff.value; + auto pe = this->peclet_f()(water_flux,water_content).value; + auto val = diff*_gamma.value*std::pow(pe,_alpha.value); + HydrodynamicDispersion eff_hdr_disp{0.}; + for (int i = 0; i < Traits::dim; ++i) + eff_hdr_disp.value[i][i] = val; + return eff_hdr_disp; + }; + } + + /// Return a map of parameter names and values for manipulation + /** \return Map: Parameter name, parameter value in this object + */ + std::multimap parameters () override + { + using Map = std::multimap; + + Map base_param = Transport::parameters(); + Map new_param({{Alpha::name, _alpha.value}, + {Gamma::name, _gamma.value}}); + Map compsite_map; + + compsite_map.insert(base_param.begin(),base_param.end()); + compsite_map.insert(new_param.begin(),new_param.end()); + + return compsite_map; + } + + /// Return a map of parameter names and values for manipulation + /** \return Map: Parameter name, parameter value in this object + */ + std::multimap parameters () const override + { + using Map = std::multimap; + + Map base_param = Transport::parameters(); + Map new_param({{Alpha::name, _alpha.value}, + {Gamma::name, _gamma.value}}); + Map compsite_map; + + compsite_map.insert(base_param.begin(),base_param.end()); + compsite_map.insert(new_param.begin(),new_param.end()); + + return compsite_map; + } + + /// Return a map of parameter names and values for manipulation + /** \return Map: Parameter name, parameter value in this object + */ + std::unique_ptr> clone () const override + { + using ThisType = PowerLawDispersion; + return std::make_unique(*this); + } +}; + +} // namespace Parameterization +} // namespace Dune +} // namespace Dorie + +#endif // DUNE_DORIE_PARAM_POWER_LAW_HD_DISP_HH \ No newline at end of file diff --git a/dune/dorie/model/transport/parameterization/hydrodynamic_dispersion/superposition.hh b/dune/dorie/model/transport/parameterization/hydrodynamic_dispersion/superposition.hh new file mode 100644 index 0000000000000000000000000000000000000000..b442bf5411d2c5d93348195b4441f18f98434ba0 --- /dev/null +++ b/dune/dorie/model/transport/parameterization/hydrodynamic_dispersion/superposition.hh @@ -0,0 +1,135 @@ +#ifndef DUNE_DORIE_PARAM_SUPERPOSITION_HD_DISP_HH +#define DUNE_DORIE_PARAM_SUPERPOSITION_HD_DISP_HH + +#include +#include +#include + +#include +#include +#include + +namespace Dune { +namespace Dorie { +namespace Parameterization { + +/*-------------------------------------------------------------------------*//** + * @brief Class for superposition hydrodynamic dispersion. + * @details Implements the relation @f$ \mathsf{D}_\mathsf{hd} + * =I \mathsf{D}^\mathsf{eff} + \mathsf{D}_\mathsf{hm}^\mathsf{eff} + * @f$ where @f$ I @f$ is the identity matrix, + * @f$ \mathsf{D}^\mathsf{eff} @f$ the effective diffusion + * coefficient, and @f$ \mathsf{D}_\mathsf{hm}^\mathsf{eff} @f$ the + * hydromechanic dispersion tensor. Both of them implemented by + * polymorphic objects following EffectiveDiffusion and + * EffectiveHydromechanicDispersion interfaces. + * + * @ingroup TransportParam + * @author Santiago Ospina De Los Ríos + * @date 2019 + * + * @tparam Traits the base traits + */ +template +class HydrodynamicDispersionSuperposition : + public Transport +{ +private: + using RangeField = typename Traits::RF; + using Base = Transport; + static constexpr int dim = Traits::dim; + +public: + using HydrodynamicDispersion = typename Base::HydrodynamicDispersion; + using WaterFlux = typename Base::WaterFlux; + using WaterContent = typename Base::WaterContent; + + /// The name of this parameterization type (This is not the instance name) + static inline std::string type = "superposition"; + + std::shared_ptr> _eff_diff; + std::shared_ptr> _eff_hm_disp; + +public: + + /// Construct with polymorphic objects + /** \param name The name associated with this soil layer + * \param eff_diff The class for the effective diffusion coefficient + * \param eff_hm_disp The class for the effective hydromechanic tensor + */ + HydrodynamicDispersionSuperposition (const std::string name, + std::shared_ptr> eff_diff, + std::shared_ptr> eff_hm_disp + ) : Base(name), + _eff_diff(eff_diff), + _eff_hm_disp(eff_hm_disp) + { } + + + /// Add default destructor to clarify override + ~HydrodynamicDispersionSuperposition () override = default; + + /// Return the hydrodynamic dispersion + /** {WaterFlux,WaterContent} -> HydrodynamicDispersion + */ + std::function + hydrodynamic_dispersion_f () const override + { + return [this](const WaterFlux water_flux, const WaterContent water_content) { + auto eff_diff = _eff_diff->effective_diffusion_f()(water_flux,water_content); + auto eff_hm_disp = _eff_hm_disp->effective_hydromechanic_dispersion_f()(water_flux,water_content); + for (int i = 0; i < dim; ++i) + eff_hm_disp.value[i][i]+=eff_diff.value; + return HydrodynamicDispersion{eff_hm_disp.value}; + }; + } + + /// Return a map of parameter names and values for manipulation + /** \return Map: Parameter name, parameter value in this object + */ + std::multimap parameters () override + { + auto base_param = Transport::parameters(); + auto eff_hydr_disp_param = _eff_hm_disp->parameters(); + auto eff_diff_param = _eff_diff->parameters(); + + std::multimap compsite_map; + compsite_map.insert(base_param.begin(),base_param.end()); + compsite_map.insert(eff_hydr_disp_param.begin(),eff_hydr_disp_param.end()); + compsite_map.insert(eff_diff_param.begin(),eff_diff_param.end()); + + return compsite_map; + } + + /// Return a map of parameter names and values for manipulation + /** \return Map: Parameter name, parameter value in this object + */ + std::multimap parameters () const override + { + auto base_param = Transport::parameters(); + auto eff_hydr_disp_param = _eff_hm_disp->parameters(); + auto eff_diff_param = _eff_diff->parameters(); + + std::multimap compsite_map; + compsite_map.insert(base_param.begin(),base_param.end()); + compsite_map.insert(eff_hydr_disp_param.begin(),eff_hydr_disp_param.end()); + compsite_map.insert(eff_diff_param.begin(),eff_diff_param.end()); + + return compsite_map; + } + + /// Return a map of parameter names and values for manipulation + /** \return Map: Parameter name, parameter value in this object + */ + std::unique_ptr> clone () const override + { + using ThisType = HydrodynamicDispersionSuperposition; + return std::make_unique(this->get_name(),_eff_diff->clone(),_eff_hm_disp->clone()); + } +}; + +} // namespace Parameterization +} // namespace Dune +} // namespace Dorie + +#endif // DUNE_DORIE_PARAM_SUPERPOSITION_HD_DISP_HH \ No newline at end of file diff --git a/dune/dorie/model/transport/parameterization/interface.hh b/dune/dorie/model/transport/parameterization/interface.hh new file mode 100644 index 0000000000000000000000000000000000000000..16764f3658d12eea8f561fc7f0d4110bfbe1aa69 --- /dev/null +++ b/dune/dorie/model/transport/parameterization/interface.hh @@ -0,0 +1,148 @@ +#ifndef DUNE_DORIE_PARAM_TRANSPORT_INTERFACE_HH +#define DUNE_DORIE_PARAM_TRANSPORT_INTERFACE_HH + +#include +#include +#include + +namespace Dune { +namespace Dorie { +namespace Parameterization { + +/// Class defining the interface for parameterizations of Transport equation +template +class Transport +{ +private: + using RangeField = typename Traits::RF; + using Vector = typename Traits::Vector; + using Tensor = typename Traits::Tensor; + +public: + // These are standard types for all parameterizations of Transport Equation + /// Type of the hydrodynamic dispersion + struct HydrodynamicDispersion + { + Tensor value; + }; + + /// Type of water content + struct WaterContent + { + RangeField value; + }; + + // Type of water flux + struct WaterFlux + { + Vector value; + }; + + /// Type of microscopic peclet number + struct Peclet + { + RangeField value; + }; + + /// Parameter defining the characteristic length + struct CharacteristicLength + { + CharacteristicLength() {} + RangeField value; + inline static const std::string name = "char_length"; + }; + + /// Parameter defining the molecluar diffusion + struct MolecularDiffusion + { + MolecularDiffusion() {} + RangeField value; + inline static const std::string name = "mol_diff"; + }; + + //! Value of the characteristic length. + CharacteristicLength _char_length; + //! Value of the molecluar diffusion. + MolecularDiffusion _mol_diff; + //! The name of this parameterization instance, associated with the layer. + const std::string _name; + +public: + /// Construct with default-initialized parameters + /** \param name The name associated with this soil layer + */ + Transport (const std::string name) : _name(name) + { } + + /// Construct from a tuple of parameters + /** \param name The name associated with this soil layer + * \param parameters Tuple of parameters to use in this parameterization + */ + template + Transport ( + const std::string name, + const std::tuple parameters) : + _char_length(std::get(parameters)), + _mol_diff(std::get(parameters)) + { } + + /// Default constructor (virtual). + virtual ~Transport () = default; + + /// Return the name of this parameterization instance. + const std::string& get_name() const { return _name; } + + /// Return a bound version of the hydrodynamic dispersion tensor + /** \return Function: {Water Flux, Water Content} -> Hydrodynamic Disp. Tensor + */ + virtual std::function< + HydrodynamicDispersion(const WaterFlux water_flux, const WaterContent water_content)> + hydrodynamic_dispersion_f () const = 0; + + /// Return a bound version of the microscopic peclet function + /** \return Function: {Water Flux, Water Content} -> Peclet + */ + std::function< + Peclet(const WaterFlux, const WaterContent)> + peclet_f () const + { + return [this](const WaterFlux water_flux, const WaterContent water_content){ + return Peclet{_char_length.value*water_flux.value.two_norm()/(_mol_diff.value*water_content.value)}; + }; + } + + /// Return a map referecing all parameters by their names. + /** \return Map. Key: Name of parameter (string). + * Value: Value of parameter (double&) + */ + virtual std::multimap parameters () + { + return { + {_char_length.name, _char_length.value}, + {_mol_diff.name, _mol_diff.value} + }; + } + + /// Return a map referecing all parameters by their names. + /** \return Map. Key: Name of parameter (string). + * Value: Value of parameter (const double&) + */ + virtual std::multimap parameters () const + { + return { + {_char_length.name, _char_length.value}, + {_mol_diff.name, _mol_diff.value} + }; + } + + /// Clone the plymorphic class + /** \return unique pointer to the cloned object + */ + virtual std::unique_ptr> clone () const = 0; +}; + +} // namespace Parameterization +} // namespace Dune +} // namespace Dorie + +#endif // DUNE_DORIE_PARAM_TRANSPORT_INTERFACE_HH \ No newline at end of file diff --git a/dune/dorie/model/transport/solute_boundary.hh b/dune/dorie/model/transport/solute_boundary.hh index 4ee21c33ae48c967c0e30e5d28cf7d9a39b85826..8306ea30770776ac1b000f2723be84a7a0887520 100644 --- a/dune/dorie/model/transport/solute_boundary.hh +++ b/dune/dorie/model/transport/solute_boundary.hh @@ -11,6 +11,17 @@ namespace Dune{ namespace Dorie{ +/** + * @brief Interior penalty type for dg methods + * @ingroup LocalOperators + */ +struct DirichletMode +{ + enum Type { SoluteConcentration, //!< Solute concentration + TotalSolute, //!< Total solute + }; +}; + /*-------------------------------------------------------------------------*//** * @brief Boundary type and condition value queries for solute. * @details This class containts functions that return the type of boundary diff --git a/dune/dorie/model/transport/solute_parameters.hh b/dune/dorie/model/transport/solute_parameters.hh new file mode 100644 index 0000000000000000000000000000000000000000..57611114d524e40e4917bf8be9f5baa664872b9c --- /dev/null +++ b/dune/dorie/model/transport/solute_parameters.hh @@ -0,0 +1,230 @@ +#ifndef DUNE_DORIE_PARAM_TRANSPORT_HH +#define DUNE_DORIE_PARAM_TRANSPORT_HH + +#include +#include +#include +#include +#include +#include +#include + +#include + +#include +#include + +#include + +namespace Dune { +namespace Dorie { + +/// Top-level parameterization interface for the local operator. +/** This class loads and stores the parameters, maps them to grid entities, + * and provides functions to access the parameterization. For doing so, the + * FlowParameters::bind function has to be called first to cache the + * parameters of the given grid entity. Parameters are only stored for entities + * of codim 0. + * + * The parameterization consists of three structures: + * - This top-level class serves as interface and storage. It also casts to + * data types used by the DUNE operators + * - The Dorie::Parameterization::Transport interface that provides strong + * data types for the base parameters and purely virtual functions that + * need to be defined by derived parameterizations + * - The actual parameterization defining all parameters and functions + * + * \author Lukas Riedel + * \author Santiago Ospina De Los Ríos + * \date 2018-2019 + * \ingroup RichardsParam + */ +template +class SoluteParameters +{ +private: + using Grid = typename Traits::Grid; + using LevelGridView = typename Grid::LevelGridView; + using Mapper = typename Dune::MultipleCodimMultipleGeomTypeMapper< + LevelGridView + >; + using Scalar = typename Traits::Scalar; + using Vector = typename Traits::Vector; + using Tensor = typename Traits::Tensor; + + // Define the storage and cache for parameterizations and scaling + /// Index type used by the element mapper + using Index = typename Mapper::Index; + /// Base class of the parameterization + using ParameterizationType = typename Dorie::Parameterization::Transport; + /// Parameterization factory + using ParameterizationFactory = Dorie::Parameterization::TransportFactory; + + /// Storage for parameters + using ParameterStorage = std::unordered_map>; + + /// Need this gruesome typedef because 'map::value_type' has const key + using Cache = std::shared_ptr; + + using ConstCache = std::shared_ptr; + /// Configuration file tree + const Dune::ParameterTree& _config; + /// Grid view of the coarsest grid configuration (level 0) + const LevelGridView _gv; + /// Logger for this instance + const std::shared_ptr _log; + /// Mapper for mapping grid elements to indices + Mapper _mapper; + /// Map for storing parameterization information for each grid entity + ParameterStorage _param; + /// Currently cached parameterization (bound to element) + mutable Cache _cache; + +private: + + /// Check if an entity has been cached + void verify_cache () const + { + if (not _cache) { + _log->error("Parameterization cache is empty. Call 'bind' before " + "accessing the cache or the parameterization"); + DUNE_THROW(Dune::InvalidStateException, + "No parameterization cached"); + } + } + +public: + + /// Copy constructor for solute parameters. + /** Create a level grid view of level 0, create a mapper on it, and store + * the config tree. + */ + SoluteParameters ( + const SoluteParameters& solute_param + ): + _config(solute_param._config), + _gv(solute_param._gv), + _log(solute_param._log), + _mapper(solute_param._mapper) + { + // copy parameterization map + this->_param.clear(); + for(const auto& [index, p] : solute_param._param) { + // make a hard copy of parameterization + std::shared_ptr> _p = p->clone(); + this->_param.emplace(index,_p); + } + } + + /// Create this object and load the data. + /** Create a level grid view of level 0, create a mapper on it, and store + * the config tree. + * \param config Configuration file tree + * \param grid Shared pointer to the grid + * \param element_index_map The mapping from grid entity index to + * parameterization index. + */ + SoluteParameters ( + const Dune::ParameterTree& config, + const std::shared_ptr grid, + const std::vector& element_index_map + ): + _config(config), + _gv(grid->levelGridView(0)), + _log(Dorie::get_logger(log_transport)), + _mapper(_gv, Dune::mcmgElementLayout()) + { + build_parameterization(element_index_map); + } + + /// Return a copy of the current cache with constant parameterization + ConstCache cache() const + { + verify_cache(); + return _cache; + } + + /// Return the current cache + const Cache& cache() + { + verify_cache(); + return _cache; + } + + /// Bind to a grid entity. Required to call parameterization functions! + /** \param entity Grid entity (codim 0) to bind to + */ + template + void bind (Entity entity) const + { + // retrieve index of top father element of chosen entity + while(entity.hasFather()) { + entity = entity.father(); + } + const auto index = _mapper.index(entity); + _cache = _param.find(index)->second; + } + + /// Return the hydrodynamic dispersion function + /** Uses the function of the underlying parameterization. Cast to + * operator-internal Tensor. + * \return Function: {WaterFlux,WaterContent} -> HydrdyDisp + */ + std::function hydrodynamic_dispersion_f () const + { + verify_cache(); + const auto hd_dips = _cache->hydrodynamic_dispersion_f(); + + using WaterFlux = typename ParameterizationType::WaterFlux; + using WaterContent = typename ParameterizationType::WaterContent; + return [hd_dips](const Vector water_flux, const Scalar water_content) { + return hd_dips(WaterFlux{water_flux},WaterContent{water_content}).value; + }; + } + + /// Return the microscopic peclet fuction + /** Uses the function of the underlying parameterization. Cast to + * operator-internal RF. + * \return Function: {WaterFlux,WaterContent} -> MicroPeclet + */ + auto peclet_f () const + { + verify_cache(); + const auto peclet = _cache->peclet_f(); + + using WaterFlux = typename ParameterizationType::WaterFlux; + using WaterContent = typename ParameterizationType::WaterContent; + return [peclet](const Vector water_flux, const Scalar water_content) { + return peclet(WaterFlux{water_flux},WaterContent{water_content}).value; + }; + } + +private: + + /// Build the parameterization from an element mapping and the input file. + void build_parameterization (const std::vector& element_index_map) + { + // Open the YAML file + const auto param_file_name = _config.get("parameters.file"); + _log->info("Loading parameter data file: {}", param_file_name); + YAML::Node param_file = YAML::LoadFile(param_file_name); + + // create the parameterization data + ParameterizationFactory factory; + const auto parameterization_map = factory.reader(param_file,"transport",_log); + + // insert parameterization for each element + for (auto&& elem : elements(_gv)) { + const auto index = _mapper.index(elem); + auto p = parameterization_map.at(element_index_map.at(index)); + _param.emplace(index,p); + } + // check that mapper can map to parameterization + assert(_param.size() == _mapper.size()); + } +}; + +} // namespace Dune +} // namespace Dorie + +#endif // DUNE_DORIE_PARAM_TRANSPORT_HH \ No newline at end of file diff --git a/dune/dorie/model/transport/transport.cc b/dune/dorie/model/transport/transport.cc index 49e48fe404253749d0f8f104448505c7a7f806ff..81710876688457a3fc2f9ebe4d22d6106205c5e1 100644 --- a/dune/dorie/model/transport/transport.cc +++ b/dune/dorie/model/transport/transport.cc @@ -44,8 +44,12 @@ TransportSimulation::TransportSimulation( cc = std::make_unique(); Dune::PDELab::constraints(*gfs,*cc,false); + // --- Create the new parameter class + auto element_map = _grid_creator.element_index_map(); + sparam = std::make_shared(inifile, grid, element_map); + // --- Operator Helper Classes --- - sboundary = std::make_unique(inifile); + sboundary = std::make_shared(inifile); sinitial = SoluteInitialFactory::create(inifile, gv, this->_log); controller = std::make_unique(inifile, @@ -102,8 +106,7 @@ void TransportSimulation::operator_setup() this->_log->debug(" Total number of DOF: {}", dof_count); // --- Local Operators --- - double diff_coeff = inifile.get("parameters.effHydrDips"); - slop = std::make_unique(sboundary, diff_coeff); + slop = std::make_unique(sparam, sboundary); tlop = std::make_unique(); Dune::PDELab::constraints(*this->gfs,*this->cc,false); @@ -219,23 +222,34 @@ void TransportSimulation::write_data () const { if (vtkwriter) { + auto peclet = std::make_shared(gv, sparam, water_flux_gf, water_content_gf); + auto d_hd = std::make_shared(gv, sparam, water_flux_gf, water_content_gf); + if (inifile.get("output.vertexData")) { - vtkwriter->template addVertexData(get_solute(),"solute"); - vtkwriter->template addVertexData(get_total_solute(),"total_solute"); + vtkwriter->addVertexData(get_solute(),"solute"); + vtkwriter->addVertexData(get_total_solute(),"total_solute"); + vtkwriter->addVertexData(peclet,"micro_peclet"); + + if (inifile.get("output.writeDispersionTensor")) + vtkwriter->addVertexData(d_hd,"eff_hd_dispersion"); + if constexpr (enable_rt_engine) if (enable_fluxrc) { auto RT_name = "flux_RT" + std::to_string(flux_order); - vtkwriter->template addVertexData - (get_solute_flux_reconstructed(),RT_name); + vtkwriter->addVertexData(get_solute_flux_reconstructed(),RT_name); } } else { - vtkwriter->template addCellData(get_solute(),"solute"); - vtkwriter->template addCellData(get_total_solute(),"total_solute"); + vtkwriter->addCellData(get_solute(),"solute"); + vtkwriter->addCellData(get_total_solute(),"total_solute"); + vtkwriter->addCellData(peclet,"micro_peclet"); + + if (inifile.get("output.writeDispersionTensor")) + vtkwriter->addCellData(d_hd,"eff_hd_dispersion"); + if constexpr (enable_rt_engine) if (enable_fluxrc) { auto RT_name = "flux_RT" + std::to_string(flux_order); - vtkwriter->template addCellData - (get_solute_flux_reconstructed(),RT_name); + vtkwriter->addCellData(get_solute_flux_reconstructed(),RT_name); } } diff --git a/dune/dorie/model/transport/transport.hh b/dune/dorie/model/transport/transport.hh index 114cae7c3672008561f21b5440e39cede0e6b380..50a2de53283bbaecf85c57185be47d16e93510e3 100644 --- a/dune/dorie/model/transport/transport.hh +++ b/dune/dorie/model/transport/transport.hh @@ -24,6 +24,9 @@ #include #include +#include +#include +#include #include #include #include @@ -81,6 +84,8 @@ struct TransportSimulationTraits : public BaseTraits // -- DORiE Class Definitions -- // /// Wrapper for grid and volume & boundary mappings using GridCreator = Dune::Dorie::GridCreator; + /// Class defining the solute parameterization + using SoluteParameters = Dune::Dorie::SoluteParameters; /// Class defining boundary conditions using SoluteBoundary = Dune::Dorie::SoluteBoundary; /// Class defining the initial condition @@ -89,7 +94,8 @@ struct TransportSimulationTraits : public BaseTraits using SoluteInitialFactory = Dune::Dorie::TransportInitialConditionFactory; /// Local spatial operator - using SLOP = Dune::Dorie::Operator::TransportFVSpatialOperator; /// Local temporal operator @@ -154,7 +160,10 @@ struct TransportSimulationTraits : public BaseTraits /// Water Flux reconstruction static constexpr bool enable_rt_engine = Dune::Dorie::RaviartThomasFluxReconstructionHelper::enable; using GFFluxReconstruction = std::conditional_t,void>; - + /// Peclet + using GFPeclet = Dune::Dorie::PecletAdapter; + /// Effective Hydrodynamic Dipspersion + using GFEffectiveHydrodynamicDispersion = Dune::Dorie::EffectiveHydrodynamicDispersionAdapter; // -- Utility Class Definitions -- // /// VTK Output writer base class @@ -211,6 +220,8 @@ protected: // -- D/*ORiE Class Definitions -- // /// Wrapper for grid and volume & boundary mappings using GridCreator = typename Traits::GridCreator; + /// Class defining the solute parameterization + using SoluteParameters = typename Traits::SoluteParameters; /// Class defining boundary conditions using SoluteBoundary = typename Traits::SoluteBoundary; /// Class defining the initial condition @@ -253,6 +264,10 @@ protected: using GFTotalSolute = typename Traits::GFTotalSolute; /// Water Flux reconstruction using GFFluxReconstruction = typename Traits::GFFluxReconstruction; + /// Peclet + using GFPeclet = typename Traits::GFPeclet; + /// Effective Hydrodynamic Dipspersion + using GFEffectiveHydrodynamicDispersion = typename Traits::GFEffectiveHydrodynamicDispersion; // // -- Utility Class Definitions -- // /// VTK Output writer base class @@ -273,6 +288,7 @@ protected: std::shared_ptr gfs; std::unique_ptr cc; + std::shared_ptr sparam; std::shared_ptr sboundary; std::unique_ptr sinitial; std::unique_ptr controller; diff --git a/dune/dorie/test/CMakeLists.txt b/dune/dorie/test/CMakeLists.txt index d2ec4be70dfcfe779d54d32b09196c452db84e4a..8f0c52ee1fe354d75dd9a9df361174f1b9229b11 100644 --- a/dune/dorie/test/CMakeLists.txt +++ b/dune/dorie/test/CMakeLists.txt @@ -1,19 +1,18 @@ # copy ALL the config files and stuff file(COPY scaling.h5 scaling.yml - test-parameterization-scaled.yml - test-parameterization-unscaled.yml + test-param-richards-scaled.yml + test-param-richards-unscaled.yml + test-param-transport.yml + test-param-factory.yml DESTINATION .) # scaling adapter test dorie_add_unit_test(SOURCES test-simulation-base.cc NAME test-simulation-base) -target_link_libraries(ut-test-simulation-base spdlog) dorie_add_unit_test(SOURCES test-grid-function-container.cc NAME test-grid-function-container) -target_link_libraries(ut-test-grid-function-container spdlog) dorie_add_unit_test(SOURCES test-interpolators.cc NAME test-interpolators) -target_link_libraries(ut-test-interpolators spdlog) dorie_add_unit_test(SOURCES test-scaling-adapters.cc NAME test-scaling-adapters) -target_link_libraries(ut-test-scaling-adapters spdlog) +dorie_add_unit_test(SOURCES test-param-factory.cc NAME test-param-factory) # grid creator test dorie_add_metaini_test( @@ -23,7 +22,6 @@ dorie_add_metaini_test( METAINI test-grid-creation.mini.in CREATED_TARGETS gc_target ) -target_link_libraries(${gc_target} spdlog) configure_file(test-grid-creation-global.mini.in test-grid-creation-global.mini) dune_add_system_test( @@ -44,31 +42,38 @@ dorie_add_metaini_test( METAINI test-initial-condition.mini.in CREATED_TARGETS ic_target ) -target_link_libraries(${ic_target} spdlog) # parameterization and scaling test dorie_add_metaini_test( UNIT_TEST - SOURCE test-parameterization.cc - BASENAME test-parameterization - METAINI test-parameterization.mini.in + SOURCE test-param-richards.cc + BASENAME test-param-richards + METAINI test-param-richards.mini.in CREATED_TARGETS param_target ) target_link_libraries(${param_target} spdlog) dorie_add_metaini_test( UNIT_TEST - SOURCE test-parameterization-scaled.cc - BASENAME test-parameterization-scaled - METAINI test-parameterization-scaled.mini.in + SOURCE test-param-richards-scaled.cc + BASENAME test-param-richards-scaled + METAINI test-param-richards-scaled.mini.in + CREATED_TARGETS param_target +) + +# parameterization and scaling test +dorie_add_metaini_test( + UNIT_TEST + SOURCE test-param-transport.cc + BASENAME test-param-transport + METAINI test-param-transport.mini.in CREATED_TARGETS param_target ) -target_link_libraries(${param_target} spdlog) -add_custom_target(test-parameterization - COMMAND ctest --output-on-failure --tests-regex ^.+test-parameterization.+$ +add_custom_target(test-param + COMMAND ctest --output-on-failure --tests-regex ^.+test-param.+$ ) -add_dependencies(test-parameterization prepare_testing) +add_dependencies(test-param prepare_testing) # run two tests in parallel dune_add_test( @@ -86,4 +91,4 @@ dune_add_test( MPI_RANKS 2 TIMEOUT 30 CMD_ARGS ${CMAKE_CURRENT_BINARY_DIR}/gc_files_0003.ini -) \ No newline at end of file +) diff --git a/dune/dorie/test/test-param-factory.cc b/dune/dorie/test/test-param-factory.cc new file mode 100644 index 0000000000000000000000000000000000000000..066b8c22160c99de62b0ce66f78887fde8cda6e6 --- /dev/null +++ b/dune/dorie/test/test-param-factory.cc @@ -0,0 +1,189 @@ +#ifdef HAVE_CONFIG_H + #include "config.h" +#endif + +#include + +#include +#include +#include + +#include + +#include "test-parameterization.hh" + +/// Return parameter values based on the material index +/** \param index index assigned to the material + * \return Map of parameters and associated values + * \throw runtime_error index is not known + */ +std::multimap get_parameter_map (const int index) +{ + if (index == 0) { // scalar + return std::multimap { + {"val", 1.} + }; + } + else if (index == 1) { // 2Dvector + return std::multimap { + {"val_x", 1.}, + {"val_y", 2.} + }; + } + else if (index == 2) { // 2Dtensor + return std::multimap { + {"val_xx", 1.}, + {"val_xy", 2.}, + {"val_yx", 3.}, + {"val_yy", 4.} + }; + } + else if (index == 3) { // 3Dvector + return std::multimap { + {"val_x", 1.}, + {"val_y", 2.}, + {"val_z", 3.} + }; + } + else if (index == 4) { // 3Dtensor + return std::multimap { + {"val_xx", 1.}, + {"val_xy", 2.}, + {"val_xz", 3.}, + {"val_yx", 4.}, + {"val_yy", 5.}, + {"val_yz", 6.}, + {"val_zx", 7.}, + {"val_zy", 8.}, + {"val_zz", 9.} + }; + } + else { + throw std::runtime_error("Unknown index " + std::to_string(index)); + } +} + +struct Param { + virtual ~Param() = default; + virtual std::multimap parameters() = 0; +}; + +struct ParamScalar : public Param { + std::multimap parameters() {return {{"val",val[0]}};} + Dune::FieldVector val; +}; + +struct Param2DVec : public Param { + std::multimap parameters() { + return { {"val_x",val[0]}, + {"val_y",val[1]} }; + } + Dune::FieldVector val; +}; + +struct Param3DVec : public Param { + std::multimap parameters() { + return { {"val_x",val[0]}, + {"val_y",val[1]}, + {"val_z",val[2]} }; + } + Dune::FieldVector val; +}; + +struct Param2DTensor : public Param { +std::multimap parameters() { + return { {"val_xx",val[0][0]}, + {"val_xy",val[0][1]}, + {"val_yx",val[1][0]}, + {"val_yy",val[1][1]} }; + } + Dune::FieldMatrix val; +}; + +struct Param3DTensor : public Param { +std::multimap parameters() { + return { {"val_xx",val[0][0]}, + {"val_xy",val[0][1]}, + {"val_xz",val[0][2]}, + {"val_yx",val[1][0]}, + {"val_yy",val[1][1]}, + {"val_yz",val[1][2]}, + {"val_zx",val[2][0]}, + {"val_zy",val[2][1]}, + {"val_zz",val[2][2]}, }; + } + Dune::FieldMatrix val; +}; + +// Instaintiate factory for Param class +struct ParamFactory : public Dune::Dorie::ParameterizationFactory +{ + std::shared_ptr selector( + const YAML::Node& type_node, + const std::string& name) const override + { + const auto type = type_node["type"].as(); + if (type == "scalar") + return std::make_shared(); + else if (type == "2Dvector") + return std::make_shared(); + else if (type == "2Dtensor") + return std::make_shared(); + else if (type == "3Dvector") + return std::make_shared(); + else if (type == "3Dtensor") + return std::make_shared(); + else + DUNE_THROW(Dune::NotImplemented, + "Parameterization {" << type << "} is not implemented"); + } +}; + + +template +bool compare_parameter_maps(const ParamMap& param_map, const ParamGetter& param_getter) +{ + bool pass = true; + for (const auto& [index,param] : param_map) + { + pass &= compare_maps(param->parameters(),param_getter(index)); + } + return pass; +} + +int main (int argc, char** argv) +{ + try{ + // initialize MPI if needed + auto& helper = Dune::MPIHelper::instance(argc, argv); + + // build the logger + const auto log = Dune::Dorie::create_base_logger(helper, + spdlog::level::trace); + + const YAML::Node param_file = YAML::LoadFile("test-param-factory.yml"); + + // Create factory and read parameters + ParamFactory factory; + auto param_map = factory.reader(param_file,"param",log); + + // Wrap map getter in a lambda + auto parameter_getter = [](int index){return get_parameter_map(index);}; + + // Compare parameters from yaml with hardcoded values + return !compare_parameter_maps(param_map,parameter_getter); + } + + catch (Dune::Exception &e) { + std::cerr << "Dune reported error: " << e << std::endl; + return 1; + } + catch (std::exception &e) { + std::cerr << "Exception thrown: " << e.what() << std::endl; + return 1; + } + catch (...) { + std::cerr << "Unknown exception!" << std::endl; + return 1; + } +} diff --git a/dune/dorie/test/test-param-factory.yml b/dune/dorie/test/test-param-factory.yml new file mode 100644 index 0000000000000000000000000000000000000000..9ed9379b948d2fc875417114ab7ca9c3fc150b56 --- /dev/null +++ b/dune/dorie/test/test-param-factory.yml @@ -0,0 +1,38 @@ +volumes: + material_0: + index: 0 + param: + type: "scalar" + parameters: + val: 1 + + material_1: + index: 1 + param: + type: "2Dvector" + parameters: + val: [1, 2] + + material_2: + index: 2 + param: + type: "2Dtensor" + parameters: + val: [1, 2, + 3, 4] + + material_3: + index: 3 + param: + type: "3Dvector" + parameters: + val: [1, 2, 3] + + material_4: + index: 4 + param: + type: "3Dtensor" + parameters: + val: [1, 2, 3, + 4, 5, 6, + 7, 8, 9] \ No newline at end of file diff --git a/dune/dorie/test/test-parameterization-scaled.cc b/dune/dorie/test/test-param-richards-scaled.cc similarity index 94% rename from dune/dorie/test/test-parameterization-scaled.cc rename to dune/dorie/test/test-param-richards-scaled.cc index 5c013655da293b86402ffc43799758f01de8b1b4..1978727c9224b21024952c67066de59b799d478e 100644 --- a/dune/dorie/test/test-parameterization-scaled.cc +++ b/dune/dorie/test/test-param-richards-scaled.cc @@ -44,7 +44,7 @@ void test_scaled_parameterization (FlowParameters& param_scaled, param_unscaled.bind(cell); // retrieve actual scaling - using Scale = Dune::Dorie::ScalingFactors::RF>; + using Scale = Dune::Dorie::Parameterization::ScalingFactors::RF>; auto scale = std::get(param_scaled.cache()); // check results of scaled functions @@ -87,13 +87,13 @@ int main (int argc, char** argv) helper, spdlog::level::from_str(log_level)); - inifile["parameters.file"] = "test-parameterization-scaled.yml"; + inifile["parameters.file"] = "test-param-richards-scaled.yml"; Dune::Dorie::FlowParameters> param_scaled(inifile, grid, index_map); // use unscaled params now - inifile["parameters.file"] = "test-parameterization-unscaled.yml"; + inifile["parameters.file"] = "test-param-richards-unscaled.yml"; Dune::Dorie::FlowParameters> param_unscaled(inifile, grid, index_map); diff --git a/dune/dorie/test/test-parameterization-scaled.mini.in b/dune/dorie/test/test-param-richards-scaled.mini.in similarity index 92% rename from dune/dorie/test/test-parameterization-scaled.mini.in rename to dune/dorie/test/test-param-richards-scaled.mini.in index c2211c5197219f4e2b1296eefdbbef1da4670ebe..0942038bcf0951ea4837edb9f2afc5a42e7e14ca 100644 --- a/dune/dorie/test/test-parameterization-scaled.mini.in +++ b/dune/dorie/test/test-param-richards-scaled.mini.in @@ -13,4 +13,4 @@ extensions = 1 1 [grid.mapping] file = "{_asset_path}/maps/cell_ids.h5" -volume = ut/param_test +volume = ut/param_test_2 diff --git a/dune/dorie/test/test-param-richards-scaled.yml b/dune/dorie/test/test-param-richards-scaled.yml new file mode 100644 index 0000000000000000000000000000000000000000..95e5a74a908e14b833b33130e0659d16b585305e --- /dev/null +++ b/dune/dorie/test/test-param-richards-scaled.yml @@ -0,0 +1,34 @@ +volumes: + sand: + index: 0 + richards: + type: MvG + parameters: + alpha: -2.3 + n: 4.17 + k0: 2.2E-5 + theta_r: 0.03 + theta_s: 0.31 + tau: -1.1 + + silt: + index: 1 + richards: + type: MvG + parameters: + alpha: -0.7 + n: 1.3 + k0: 1.0E-5 + theta_r: 0.01 + theta_s: 0.41 + tau: 0.0 + +scaling: + type: Miller + data: + scale_miller: + file: scaling.h5 + dataset: twos + interpolation: nearest + extensions: [1, 1] + offset: [0, 0] diff --git a/dune/dorie/test/test-param-richards-unscaled.yml b/dune/dorie/test/test-param-richards-unscaled.yml new file mode 100644 index 0000000000000000000000000000000000000000..df866c0085619b3f5af52bd75977efc1d95dee0c --- /dev/null +++ b/dune/dorie/test/test-param-richards-unscaled.yml @@ -0,0 +1,27 @@ +volumes: + sand: + index: 0 + richards: + type: MvG + parameters: + alpha: -2.3 + n: 4.17 + k0: 2.2E-5 + theta_r: 0.03 + theta_s: 0.31 + tau: -1.1 + + silt: + index: 1 + richards: + type: MvG + parameters: + alpha: -0.7 + n: 1.3 + k0: 1.0E-5 + theta_r: 0.01 + theta_s: 0.41 + tau: 0.0 + +scaling: + type: none diff --git a/dune/dorie/test/test-parameterization.cc b/dune/dorie/test/test-param-richards.cc similarity index 78% rename from dune/dorie/test/test-parameterization.cc rename to dune/dorie/test/test-param-richards.cc index d918a7f1fc265401da3dd1035bc4235cb4b7ed56..c6dbcf265488d0df5f6e00812a32a21b7d767169 100644 --- a/dune/dorie/test/test-parameterization.cc +++ b/dune/dorie/test/test-param-richards.cc @@ -26,10 +26,10 @@ * \return Map of parameters and associated values * \throw runtime_error Medium name is not known */ -std::map get_parameter_map (const std::string medium) +std::multimap get_parameter_map (const std::string medium) { if (medium == "sand") { - return std::map { + return std::multimap { {"alpha", -2.3}, {"n", 4.17}, {"k0", 2.2e-5}, @@ -39,7 +39,7 @@ std::map get_parameter_map (const std::string medium) }; } else if (medium == "silt") { - return std::map { + return std::multimap { {"alpha", -0.7}, {"n", 1.3}, {"k0", 1.0e-5}, @@ -53,33 +53,6 @@ std::map get_parameter_map (const std::string medium) } } -/// Compare the key-values pairs for two maps -/** This does not check if the set of keys compares equal. - * map_b may have more keys than map_a. - * \return True of the values match - */ -template -bool compare_maps (const MapA& map_a, const MapB& map_b) -{ - for (auto&& [key, value_a]: map_a) { - // try to access the second map - try { - auto value_b = map_b.at(key); - // compare values - if (not Dune::FloatCmp::eq(value_a, value_b)) { - return false; - } - } - catch (std::out_of_range& e){ - return false; - } - catch (...) { - throw; - } - } - return true; -} - /// Test the parameters and parameterization functions /** The standard parameters for Sand and Silt are hardcoded here. * We check the actual values and then the return values of the @@ -101,7 +74,7 @@ void test_new_parameters (const FlowParameters& fparam, { // bind parameterization to cell and fetch parameters fparam.bind(cell); - using RP = Dune::Dorie::RichardsParameterization>; + using RP = Dune::Dorie::Parameterization::Richards>; auto par = std::get>(fparam.cache()); auto parameters = par->parameters(); @@ -116,14 +89,14 @@ void test_new_parameters (const FlowParameters& fparam, assert(FloatCmp::eq(sat_f(0.0), 1.0)); auto cond_f = fparam.conductivity_f(); - assert(FloatCmp::eq(cond_f(1.0), param_sand.at("k0")) - or FloatCmp::eq(cond_f(1.0), param_silt.at("k0"))); + assert(FloatCmp::eq(cond_f(1.0), param_sand.find("k0")->second) + or FloatCmp::eq(cond_f(1.0), param_silt.find("k0")->second)); auto wc_f = fparam.water_content_f(); - assert(FloatCmp::eq(wc_f(1.0), param_sand.at("theta_s")) - or FloatCmp::eq(wc_f(1.0), param_silt.at("theta_s"))); - assert(FloatCmp::eq(wc_f(0.0), param_sand.at("theta_r")) - or FloatCmp::eq(wc_f(0.0), param_silt.at("theta_r"))); + assert(FloatCmp::eq(wc_f(1.0), param_sand.find("theta_s")->second) + or FloatCmp::eq(wc_f(1.0), param_silt.find("theta_s")->second)); + assert(FloatCmp::eq(wc_f(0.0), param_sand.find("theta_r")->second) + or FloatCmp::eq(wc_f(0.0), param_silt.find("theta_r")->second)); } } @@ -143,16 +116,16 @@ void test_parameter_manipulation (FlowParameters& fparam, fparam.bind(cell); // manipulate via returned map (this is the preferred way) - using RP = Dune::Dorie::RichardsParameterization>; + using RP = Dune::Dorie::Parameterization::Richards>; auto par = std::get>(fparam.cache()); auto parameters = par->parameters(); - parameters.at("theta_r") = 0.0; + parameters.find("theta_r")->second = 0.0; assert(par->_theta_r.value == 0.0); // manipulate directly via casting (this is not recommended) - using MvG = Dune::Dorie::MualemVanGenuchtenParameterization>; + using MvG = Dune::Dorie::Parameterization::MualemVanGenuchten>; auto& par_mvg = dynamic_cast(*par); - parameters.at("alpha") = -10.0; + parameters.find("alpha")->second = -10.0; assert(par_mvg._alpha.value = -10.0); } @@ -168,7 +141,7 @@ bool compare_parameters (const FlowParameters& fparam1, // bind parameterization to cell and fetch parameters fparam1.bind(cell); fparam2.bind(cell); - using RP = Dune::Dorie::RichardsParameterization>; + using RP = Dune::Dorie::Parameterization::Richards>; auto par1 = std::get>(fparam1.cache()); auto par2 = std::get>(fparam2.cache()); auto parameters1 = par1->parameters(); @@ -185,7 +158,7 @@ int main (int argc, char** argv) { try{ // Initialize ALL the things! - auto [inifile, log, helper] = Dune::Dorie::Setup::init(argc, argv); + auto [inifile, log, helper] = Dune::Dorie::Setup::init(argc, argv); // create the grid Dune::Dorie::GridCreator::Grid> gc(inifile, helper); @@ -203,7 +176,7 @@ int main (int argc, char** argv) const Dune::Dorie::FlowParameters> const_param(inifile, grid, index_map); - Dune::Dorie::FlowParameters> param(const_param); + Dune::Dorie::FlowParameters> param(const_param); // deep copy // perform the actual tests test_new_parameters(const_param, grid); @@ -228,4 +201,4 @@ int main (int argc, char** argv) std::cerr << "Exception occurred!" << std::endl; throw; } -} +} \ No newline at end of file diff --git a/dune/dorie/test/test-parameterization.mini.in b/dune/dorie/test/test-param-richards.mini.in similarity index 65% rename from dune/dorie/test/test-parameterization.mini.in rename to dune/dorie/test/test-param-richards.mini.in index e8b1f6d2673ae952537f788f6b45fe25849fc45e..073382e6fce3b99edff76714833736e3a5804ed3 100644 --- a/dune/dorie/test/test-parameterization.mini.in +++ b/dune/dorie/test/test-param-richards.mini.in @@ -1,6 +1,6 @@ include ${CMAKE_BINARY_DIR}/doc/default_files/config.ini -__name = param +__name = richards_param _asset_path = "${PROJECT_SOURCE_DIR}/test" [grid] @@ -13,13 +13,14 @@ extensions = 1 1 [grid.mapping] file = "{_asset_path}/maps/cell_ids.h5" -volume = ut/param_test +volume = ut/param_test_2 [richards] -output.fileName = param | unique name -output.outputPath = param | unique name +output.fileName = {__name} +output.outputPath = {__name} +output.logLevel = trace boundary.file = "{_asset_path}/bcs/infiltration_2d.dat" -parameters.file = test-parameterization-unscaled.yml +parameters.file = test-param-richards-unscaled.yml \ No newline at end of file diff --git a/dune/dorie/test/test-param-transport.cc b/dune/dorie/test/test-param-transport.cc new file mode 100644 index 0000000000000000000000000000000000000000..6f63925fa3328f2d3e3e6aaacbbb289ee27c425e --- /dev/null +++ b/dune/dorie/test/test-param-transport.cc @@ -0,0 +1,347 @@ +#ifdef HAVE_CONFIG_H + #include "config.h" +#endif + +#include +#include +#include +#include +#include + +#include +#include +#include +#include +#include + +#include +#include +#include +#include + +#include "test-parameterization.hh" + +/// Return parameter values based on the material index +/** \param index index assigned to the material + * \return Map of parameters and associated values + * \throw runtime_error index is not known + */ +std::multimap get_parameter_map (const int index) +{ + if (index == 0) { + return std::multimap { + {"char_length", 1.5E-11}, + {"mol_diff", 2.E-9}, + {"hydrodynamic_disp_xx", 2.E-8}, + {"hydrodynamic_disp_yy", 2.E-8}, + {"hydrodynamic_disp_xy", 0}, + {"hydrodynamic_disp_yx", 0} + }; + } + else if (index == 1) { + return std::multimap { + {"char_length", 1.5E-11}, + {"mol_diff", 2.E-9}, + {"gamma", 0.8}, + {"alpha", 1.17} + }; + } + else if (index == 2) { + return std::multimap { + {"char_length", 1.5E-11}, + {"mol_diff", 2.E-9}, + {"eff_diff", 2.E-8}, + {"eff_hydromechanic_disp_xx", 2.E-8}, + {"eff_hydromechanic_disp_yy", 2.E-8}, + {"eff_hydromechanic_disp_xy", 0}, + {"eff_hydromechanic_disp_yx", 0} + }; + } + else if (index == 3) { + return std::multimap { + {"char_length", 1.5E-11}, + {"mol_diff", 2.E-9}, + {"phi", 0.34}, + {"lambda_t", 0.0025}, + {"lambda_l", 0.025} + }; + } + else if (index == 4) { + return std::multimap { + {"char_length", 1.5E-11}, + {"mol_diff", 2.E-9}, + {"phi", 0.34}, + {"eff_hydromechanic_disp_xx", 2.E-8}, + {"eff_hydromechanic_disp_yy", 2.E-8}, + {"eff_hydromechanic_disp_xy", 0}, + {"eff_hydromechanic_disp_yx", 0} + }; + } + else { + throw std::runtime_error("Unknown index " + std::to_string(index)); + } +} + +template +bool compare_parameter_maps(const ParamMap& param_map, const ParamGetter& param_getter) +{ + bool pass = true; + for (const auto& [index,param] : param_map) + { + pass &= compare_maps(param->parameters(),param_getter(index)); + } + return pass; +} + +/// Test the parameters and parameterization functions +/** The standard parameters for Sand and Silt are hardcoded here. + * We check the actual values and then the return values of the + * parameterization fuctions. + * \param sparam Dune::Dorie::SoluteParameters object + * \param grid Shared pointer to the grid + */ +template +void test_new_parameters (const SoluteParameters& sparam, + std::shared_ptr grid) +{ + // iterate over grid cells and verify parameterization + auto level_gv = grid->leafGridView(); + for (auto&& cell: elements(level_gv)) + { + constexpr int dim = Grid::dimension; + using Scalar = Dune::FieldVector; + using Vector = Dune::FieldVector; + using Tensor = Dune::FieldMatrix; + + // bind parameterization to cell and fetch parameters + sparam.bind(cell); + auto par = sparam.cache(); + auto name = par->get_name(); + auto parameters = par->parameters(); + + Scalar wc(0.5); + Vector flux(0.); + flux[dim-1] = -1e-8; + flux[0] = 1e-8; + auto fluxn = flux.two_norm(); + // Compute peclet and dispersion with some parameters + Tensor disp = sparam.hydrodynamic_dispersion_f()(flux,wc); + Scalar pe = sparam.peclet_f()(flux,wc); + + namespace FloatCmp = Dune::FloatCmp; + + // compare values to defined parameter sets + if (name == "material_0") { // const hd disp + assert(compare_maps(get_parameter_map(0), parameters)); + + for (int i = 0; i < dim; ++i) + for (int j = 0; j < dim; ++j) + if (i==j){ + assert(FloatCmp::eq(disp[i][j],2.E-8)); + } + else + assert(FloatCmp::eq(disp[i][j],0.)); + auto char_len = parameters.find("char_length")->second; + auto mol_diff = parameters.find("mol_diff")->second; + assert(FloatCmp::eq(pe,char_len*fluxn/(mol_diff*wc))); + } + else if (name == "material_1") // power_law hd disp + { + assert(compare_maps(get_parameter_map(1), parameters)); + + auto alpha = parameters.find("alpha")->second; + auto gamma = parameters.find("gamma")->second; + auto mol_diff = parameters.find("mol_diff")->second; + + auto _disp = mol_diff*gamma*std::pow(pe,alpha); + + for (int i = 0; i < dim; ++i) + for (int j = 0; j < dim; ++j) + if (i==j) + assert(FloatCmp::eq(disp[i][j],_disp)); + else + assert(FloatCmp::eq(disp[i][j],0.)); + + auto char_len = parameters.find("char_length")->second; + assert(FloatCmp::eq(pe,char_len*fluxn/(mol_diff*wc))); + } + else if (name == "material_2") // const diff + const hm disp + { + assert(compare_maps(get_parameter_map(2), parameters)); + for (int i = 0; i < dim; ++i) + for (int j = 0; j < dim; ++j) + if (i==j) + assert(FloatCmp::eq(disp[i][j],2*2.E-8)); + else + assert(FloatCmp::eq(disp[i][j],0.)); + + auto char_len = parameters.find("char_length")->second; + auto mol_diff = parameters.find("mol_diff")->second; + assert(FloatCmp::eq(pe,char_len*fluxn/(mol_diff*wc))); + } + else if (name == "material_3") // MC1 + iso + { + assert(compare_maps(get_parameter_map(3), parameters)); + + // MC1 + auto mol_diff = parameters.find("mol_diff")->second; + auto porosity = parameters.find("phi")->second; + auto diff = mol_diff*std::pow(wc,7./3.)/std::pow(porosity,2./3.); + + // iso + auto lt = parameters.find("lambda_t")->second; + auto ll = parameters.find("lambda_l")->second; + Tensor _disp(0.); + for (int i = 0; i < dim; ++i) { + for (int j = 0; j < dim; ++j) { + _disp[i][j] += (ll-lt)*flux[i]*flux[j]/(fluxn*wc); + if (i==j) + _disp[i][j] += lt*fluxn/wc + diff; + } + } + + for (int i = 0; i < dim; ++i) + for (int j = 0; j < dim; ++j) + assert(FloatCmp::eq(disp[i][j],_disp[i][j])); + + auto char_len = parameters.find("char_length")->second; + + assert(FloatCmp::eq(pe,char_len*fluxn/(mol_diff*wc))); + } + else if (name == "material_4") // MC2 + const hm disp + { + assert(compare_maps(get_parameter_map(4), parameters)); + + // MC2 + auto mol_diff = parameters.find("mol_diff")->second; + auto porosity = parameters.find("phi")->second; + auto diff = mol_diff*wc/std::pow(porosity,2./3.); + + // const hm disp + for (int i = 0; i < dim; ++i) + for (int j = 0; j < dim; ++j) + if (i==j) + assert(FloatCmp::eq(disp[i][j],double(diff)+2.E-8)); + else + assert(FloatCmp::eq(disp[i][j],0.)); + + auto char_len = parameters.find("char_length")->second; + + assert(FloatCmp::eq(pe,char_len*fluxn/(mol_diff*wc))); + } + else + { + DUNE_THROW(Dune::IOError,"Transport parameterization has an error"); + } + } +} + +/// Check if parameters can be manipulated +/** Manipulate parameters using the `parameters()` map or + * a dynamic cast to the derived class. + * \param sparam Dune::Dorie::SoluteParameters object + * \param grid Shared pointer to the grid + */ +template +void test_parameter_manipulation (SoluteParameters& sparam, + std::shared_ptr grid) +{ + // bind to any cell + auto level_gv = grid->leafGridView(); + const auto& cell = *elements(level_gv).begin(); + sparam.bind(cell); + + // manipulate via returned map (this is the preferred way) + auto par = sparam.cache(); + auto parameters = par->parameters(); + + // modify all keys + for (auto it : parameters) + it.second = 0.; + + { // modify specific key + auto [begin,end] = parameters.equal_range("mol_diff"); + for (auto p = begin; p != end; ++p) + p->second = 2.5E-9; + assert(par->_mol_diff.value == 2.5E-9); + } + { // modify specific key + auto [begin,end] = parameters.equal_range("char_length"); + for (auto p = begin; p != end; ++p) + p->second = 1E-11; + assert(par->_char_length.value == 1E-11); + } + +} + +template +bool compare_parameters (const SoluteParameters& sparam1, + const SoluteParameters& sparam2, + const std::shared_ptr& grid) +{ + // iterate over grid cells and verify parameterization + auto level_gv = grid->leafGridView(); + for (auto&& cell: elements(level_gv)) + { + // bind parameterization to cell and fetch parameters + sparam1.bind(cell); + sparam2.bind(cell); + auto par1 = sparam1.cache(); + auto par2 = sparam2.cache(); + auto parameters1 = par1->parameters(); + auto parameters2 = par2->parameters(); + + // compare values between the two parameterizations + if (not compare_maps(parameters1, parameters2)) + return false; + } + return true; +} + +int main (int argc, char** argv) +{ + try{ + // Initialize ALL the things! + auto [inifile, log, helper] = Dune::Dorie::Setup::init(argc, argv); + + // create the grid + Dune::Dorie::GridCreator::Grid> gc(inifile, helper); + auto grid = gc.grid(); + const auto index_map = gc.element_index_map(); + + inifile = Dune::Dorie::Setup::prep_ini_for_transport(inifile); + + // create the Richards logger because SoluteParameters expects it + const auto log_level = inifile.get("output.logLevel"); + log = Dune::Dorie::create_logger(Dune::Dorie::log_transport, + helper, + spdlog::level::from_str(log_level)); + + const Dune::Dorie::SoluteParameters> const_param(inifile, + grid, + index_map); + Dune::Dorie::SoluteParameters> param(const_param); // deep copy + // perform the actual tests + test_new_parameters(const_param, grid); + test_new_parameters(param, grid); + + assert(compare_parameters(const_param, param, grid)); + test_parameter_manipulation(param, grid); + assert(not compare_parameters(const_param, param, grid)); + + return 0; + } + + catch (Dune::Exception &e){ + std::cerr << "Dune reported error: " << e << std::endl; + return 1; + } + catch (std::exception &e) { + std::cerr << "Exception thrown: " << e.what() << std::endl; + return 1; + } + catch(...) { + std::cerr << "Exception occurred!" << std::endl; + throw; + } +} \ No newline at end of file diff --git a/dune/dorie/test/test-param-transport.mini.in b/dune/dorie/test/test-param-transport.mini.in new file mode 100644 index 0000000000000000000000000000000000000000..6c129b625e3719b1ea2629d2dfa1c0cfb4a4e8da --- /dev/null +++ b/dune/dorie/test/test-param-transport.mini.in @@ -0,0 +1,26 @@ +include ${CMAKE_BINARY_DIR}/doc/default_files/config.ini + +__name = transport_param +_asset_path = "${PROJECT_SOURCE_DIR}/test" + +[grid] +dimensions = 2 +initialLevel = 1 +gridType = rectangular +cells = 50 50 +FEorder = 1 +extensions = 1 1 + +[grid.mapping] +file = "{_asset_path}/maps/cell_ids.h5" +volume = ut/param_test_5 + +[transport] + +output.fileName = {__name} +output.outputPath = {__name} +output.logLevel = trace + +boundary.file = "{_asset_path}/bcs/infiltration_2d.dat" + +parameters.file = test-param-transport.yml \ No newline at end of file diff --git a/dune/dorie/test/test-param-transport.yml b/dune/dorie/test/test-param-transport.yml new file mode 100644 index 0000000000000000000000000000000000000000..b3069332576812b0078d25a82e401cf6a228af08 --- /dev/null +++ b/dune/dorie/test/test-param-transport.yml @@ -0,0 +1,53 @@ +volumes: + material_0: + index: 0 + transport: + type: Dhd_const + parameters: + char_length: 1.5E-11 + mol_diff: 2.E-9 + hydrodynamic_disp: [2.E-8, 0., + 0. , 2.E-8] + + material_1: + index: 1 + transport: + type: Dhd_pl + parameters: + char_length: 1.5E-11 + mol_diff: 2.E-9 + gamma: 0.8 + alpha: 1.17 + + material_2: + index: 2 + transport: + type: Deff_const + Dhm_const + parameters: + char_length: 1.5E-11 + mol_diff: 2.E-9 + eff_diff: 2.E-8 + eff_hydromechanic_disp: [2.E-8, 0., + 0. , 2.E-8] + + material_3: + index: 3 + transport: + type: Deff_MQ1 + Dhm_iso + parameters: + char_length: 1.5E-11 + mol_diff: 2.E-9 + phi: 0.34 + lambda_t: 0.0025 + lambda_l: 0.025 + + material_4: + index: 4 + transport: + type: Deff_MQ2 + Dhm_const + parameters: + char_length: 1.5E-11 + mol_diff: 2.E-9 + phi: 0.34 + eff_hydromechanic_disp: [2.E-8, 0., + 0. , 2.E-8] \ No newline at end of file diff --git a/dune/dorie/test/test-parameterization-scaled.yml b/dune/dorie/test/test-parameterization-scaled.yml deleted file mode 100644 index 687d81bac557c9bda1226b36a105ebac9cfe7a70..0000000000000000000000000000000000000000 --- a/dune/dorie/test/test-parameterization-scaled.yml +++ /dev/null @@ -1,32 +0,0 @@ -volumes: - sand: - index: 0 - type: MvG - parameters: - alpha: -2.3 - n: 4.17 - k0: 2.2E-5 - theta_r: 0.03 - theta_s: 0.31 - tau: -1.1 - - silt: - index: 1 - type: MvG - parameters: - alpha: -0.7 - n: 1.3 - k0: 1.0E-5 - theta_r: 0.01 - theta_s: 0.41 - tau: 0.0 - -scaling: - type: Miller - data: - scale_miller: - file: scaling.h5 - dataset: twos - interpolation: nearest - extensions: [1, 1] - offset: [0, 0] diff --git a/dune/dorie/test/test-parameterization-unscaled.yml b/dune/dorie/test/test-parameterization-unscaled.yml deleted file mode 100644 index 905364f13a5eefb3e18da56a16cb6d344ddf768a..0000000000000000000000000000000000000000 --- a/dune/dorie/test/test-parameterization-unscaled.yml +++ /dev/null @@ -1,25 +0,0 @@ -volumes: - sand: - index: 0 - type: MvG - parameters: - alpha: -2.3 - n: 4.17 - k0: 2.2E-5 - theta_r: 0.03 - theta_s: 0.31 - tau: -1.1 - - silt: - index: 1 - type: MvG - parameters: - alpha: -0.7 - n: 1.3 - k0: 1.0E-5 - theta_r: 0.01 - theta_s: 0.41 - tau: 0.0 - -scaling: - type: none diff --git a/dune/dorie/test/test-parameterization.hh b/dune/dorie/test/test-parameterization.hh index 96c385eb9255452a39f2dcc58238929143fc46a5..c5304705f93f8fd11d949cf9db73974c6c2b89c8 100644 --- a/dune/dorie/test/test-parameterization.hh +++ b/dune/dorie/test/test-parameterization.hh @@ -2,16 +2,53 @@ #define DUNE_DORIE_TEST_PARAMETERIZATION_HH #include +#include +#include #include +/// Compare the key-values pairs for two maps +/** This does not check if the set of keys compares equal. + * map_b may have more keys than map_a. + * \return True of the values match + */ +template +bool compare_maps (const MapA& map_a, const MapB& map_b) +{ + auto it_a = map_a.begin(); + auto it_b = map_b.begin(); + while (it_a!=map_a.end() && it_b!=map_b.end()) + { + auto&& [key_a,value_a] = *it_a; + auto&& [key_b,value_b] = *it_b; + + if (key_a struct Traits { static constexpr int dim = dimension; using RF = double; + using Scalar = Dune::FieldVector; using Vector = Dune::FieldVector; + using Tensor = Dune::FieldMatrix; using Grid = Dune::YaspGrid; using DF = typename Grid::ctype; @@ -19,4 +56,4 @@ struct Traits using GridView = typename Grid::LeafGridView; }; -#endif // DUNE_DORIE_TEST_PARAMETERIZATION_HH +#endif // DUNE_DORIE_TEST_PARAMETERIZATION_HH \ No newline at end of file diff --git a/dune/dorie/test/test-scaling-adapters.cc b/dune/dorie/test/test-scaling-adapters.cc index 2323bbf4469ce4acd5ad581741e0f5402f7254b1..0c6ea802c9e13379587c3af0dd88fa7b775d6692 100644 --- a/dune/dorie/test/test-scaling-adapters.cc +++ b/dune/dorie/test/test-scaling-adapters.cc @@ -23,8 +23,8 @@ /// Compare two scaling factor structs template std::enable_if_t, bool> operator== -(const Dune::Dorie::ScalingFactors& a, - const Dune::Dorie::ScalingFactors& b) +(const Dune::Dorie::Parameterization::ScalingFactors& a, + const Dune::Dorie::Parameterization::ScalingFactors& b) { if (a.scale_cond == b.scale_cond and a.scale_head == b. scale_head @@ -56,7 +56,7 @@ void test_scaling_setting (const YAML::Node& cfg, { // create the scaling const auto type = cfg["type"].as(); - using SAF = Dune::Dorie::ScalingAdapterFactory>; + using SAF = Dune::Dorie::Parameterization::ScalingAdapterFactory>; const auto log = Dune::Dorie::get_logger(Dune::Dorie::log_base); auto scaling = SAF::create(type, cfg["data"], grid_view, log); @@ -77,7 +77,7 @@ void evaluate_scaling (const std::string& type, const std::vector& positions) { // create reference scaling - using Scaling = Dune::Dorie::ScalingFactors; + using Scaling = Dune::Dorie::Parameterization::ScalingFactors; Scaling reference; if (type == "none") { reference = Scaling{1.0, 1.0, 0.0}; diff --git a/python/dorie/dorie/__init__.py b/python/dorie/dorie/__init__.py index de40ea7ca058e07a399acf61529d418c07eeee61..e69de29bb2d1d6434b8b29ae775ad8c2e48c5391 100644 --- a/python/dorie/dorie/__init__.py +++ b/python/dorie/dorie/__init__.py @@ -1 +0,0 @@ -__import__('pkg_resources').declare_namespace(__name__) diff --git a/python/dorie/dorie/cli/cmds.py.in b/python/dorie/dorie/cli/cmds.py.in index 061e2d56e31577d7710b4c4895b50e8caec7d45d..6de815d2962e552a8ed0e5e2e0cc5be754dc76ff 100644 --- a/python/dorie/dorie/cli/cmds.py.in +++ b/python/dorie/dorie/cli/cmds.py.in @@ -72,7 +72,7 @@ def run_transport(args): sys.exit(1) def create(args): - for f in ("config.ini", "parfield.ini", "2d_infiltr.bcdat", "3d_infiltr.bcdat", "2d_solute.bcdat", "3d_solute.bcdat", "param.yml"): + for f in ("config.ini", "parfield.ini", "2d_infiltr.bcdat", "3d_infiltr.bcdat", "2d_solute.bcdat", "3d_solute.bcdat", "richards_param.yml", "transport_param.yml"): newfile = os.path.join(os.getcwd(),f) if os.path.exists(newfile) and not args["force"]: override = (input("File {} exists. Override? [Y/n] ".format(f)) or "Y") in ["Y","y"] diff --git a/python/dorie/dorie/parfield/converter/base.py b/python/dorie/dorie/parfield/converter/base.py index 21c9a4566812154ed8be3ba6aa40f139d55175e2..281dcdccf4826a3136aa0a515d7a115763b420e1 100644 --- a/python/dorie/dorie/parfield/converter/base.py +++ b/python/dorie/dorie/parfield/converter/base.py @@ -1,9 +1,10 @@ import h5py import os import sys -import configparser import numpy as np +from dune.testtools.parser import parse_ini_file + from dorie.utilities.text_to_bool import text_to_bool class BaseConverter(object): @@ -20,9 +21,7 @@ class BaseConverter(object): _dataset = None # The field to be written into the target file def __init__(self, param, input_file=None): - self._cp = configparser.ConfigParser(interpolation=None) - with open(param) as file: - self._cp.read_file(file) + self._cp = parse_ini_file(param) outfile = self._read_parameter('general', 'outputFile') dataset = self._read_parameter('general', 'dataset') diff --git a/python/dorie/dorie/test/__init__.py b/python/dorie/dorie/test/__init__.py new file mode 100644 index 0000000000000000000000000000000000000000..e69de29bb2d1d6434b8b29ae775ad8c2e48c5391 diff --git a/python/dorie/dorie/test/test_vtkreader.py b/python/dorie/dorie/test/test_vtkreader.py index c4706aeb2a61bac580a425a871cc0a2166004ffc..b94fadf71e8e1ff5c397b956d27217041c7c44b1 100644 --- a/python/dorie/dorie/test/test_vtkreader.py +++ b/python/dorie/dorie/test/test_vtkreader.py @@ -1,7 +1,7 @@ import pytest import numpy as np -from dorie.utilities.vtktools.vtkreader import VTKReader +from ..utilities.vtktools.vtkreader import VTKReader, VTKDataArray VTK_CELL = "vtkreader-cell.vtu" VTK_VERTEX = "vtkreader-vertex.vtu" @@ -12,65 +12,78 @@ ARRAYS = ["head", "flux", "K_0", "Theta", "theta_w"] @pytest.fixture def vtk_cell (): """Create a VTKReader for the cell VTK test file""" - return VTKReader(VTK_CELL, DEFAULT_ARRAY) + return VTKReader(VTK_CELL) @pytest.fixture def vtk_vertex (): """Create a VTKReader for the vertex VTK test file""" - return VTKReader(VTK_VERTEX, DEFAULT_ARRAY) + return VTKReader(VTK_VERTEX) # Tests ------------------------------------------------------------------ def test_init (): """Test if initialization works""" - VTKReader(VTK_CELL, DEFAULT_ARRAY) - VTKReader(VTK_VERTEX, DEFAULT_ARRAY) + reader = VTKReader(VTK_CELL) + isinstance(reader[DEFAULT_ARRAY], VTKDataArray) + reader = VTKReader(VTK_VERTEX) + isinstance(reader.get(DEFAULT_ARRAY), VTKDataArray) -def test_array_set (vtk_cell, vtk_vertex): +def test_dict_properties (vtk_cell, vtk_vertex): """Test if all included arrays are found""" + assert len(vtk_cell) == len(ARRAYS) + for array in ARRAYS: - vtk_cell.set_data_array(array) - vtk_vertex.set_data_array(array) + assert array in vtk_cell + assert array in vtk_vertex.keys() + + for key in vtk_cell: + assert key in ARRAYS + + for key, array in vtk_cell.items(): + assert key in ARRAYS - with pytest.raises(RuntimeError): - vtk_cell.set_data_array("other") + assert not "other" in vtk_cell + with pytest.raises(KeyError): + array = vtk_vertex["other"] -def test_properties (vtk_cell): - """Test the VTKReader properties""" - assert np.allclose(vtk_cell.bounds, [[0, 1], - [0, 1], - [0, 0]]) - assert vtk_cell.number_of_components == 1 +def test_array_properties (vtk_cell): + """Test the VTKArray properties""" + array = vtk_cell[DEFAULT_ARRAY] + assert array.number_of_components == 1 + assert np.allclose(array.bounds, [[0, 1], + [0, 1], + [0, 0]]) - vtk_cell.set_data_array("flux") - assert vtk_cell.number_of_components == 3 + assert vtk_cell["flux"].number_of_components == 3 def test_evaluate (vtk_cell, vtk_vertex): """Test the evaluation function""" - center = np.mean(vtk_cell.bounds, axis=1) + array_cell = vtk_cell[DEFAULT_ARRAY] + array_vertex = vtk_vertex[DEFAULT_ARRAY] + center = np.mean(array_cell.bounds, axis=1) # should return scalar with pytest.raises(TypeError): - len(vtk_cell.evaluate(center)) + len(array_cell.evaluate(center)) # check return of vector - vtk_vertex.set_data_array("flux") - assert len(vtk_vertex.evaluate(center)) == 3 + array_vertex = vtk_vertex["flux"] + assert len(array_vertex.evaluate(center)) == 3 # check values are equal for array in ["head", "K_0", "flux"]: - vtk_cell.set_data_array(array) - vtk_vertex.set_data_array(array) - assert np.allclose(vtk_cell.evaluate(center), - vtk_vertex.evaluate(center)) + array_cell = vtk_cell[array] + array_vertex = vtk_vertex[array] + assert np.allclose(array_cell.evaluate(center), + array_vertex.evaluate(center)) # check evaluation errors # position has wrong shape with pytest.raises(ValueError): - vtk_vertex.evaluate([0.0, 0.0]) + array_vertex.evaluate([0.0, 0.0]) # position out of bounds - pos = vtk_vertex.bounds[:, 1] + pos = array_vertex.bounds[:, 1] pos = pos + [1.0, 0.0, 1.0] with pytest.raises(ValueError): - vtk_vertex.evaluate(pos) + array_vertex.evaluate(pos) diff --git a/python/dorie/dorie/testtools/dorie_pfg/__init__.py b/python/dorie/dorie/testtools/dorie_pfg/__init__.py new file mode 100644 index 0000000000000000000000000000000000000000..e69de29bb2d1d6434b8b29ae775ad8c2e48c5391 diff --git a/python/dorie/dorie/testtools/dorie_run/ode.py b/python/dorie/dorie/testtools/dorie_run/ode.py index fd774ed3a8c3c516c1316fd2925d94de9ffbe39c..b4d7efda567a2a68f054c312028c8efa4d3f877a 100644 --- a/python/dorie/dorie/testtools/dorie_run/ode.py +++ b/python/dorie/dorie/testtools/dorie_run/ode.py @@ -38,13 +38,14 @@ def evaluate(iniinfo,runtime): if mapping_file in ["None", "none", ""]: medium_index = int(mapping_dataset) param = get_parameter_data(medium_index, volumes=param_data["volumes"]) - if param["type"] != "MvG": + param_type = param["richards"]["type"] + if param_type != "MvG": raise NotImplementedError("This testing function only supports " "Mualem van Genuchten parameterizations. " "Parameterization was: {}".format( - param["type"])) + param_type)) - parfun = lambda y: param["parameters"] + parfun = lambda y: param["richards"]["parameters"] # 1D heterogeneity else: @@ -60,12 +61,13 @@ def evaluate(iniinfo,runtime): cell = math.floor(y * mapping.shape[0] / extensions[-1]) index = mapping[cell, ...].flatten() param = get_parameter_data(index, volumes=param_data["volumes"]) - if param["type"] != "MvG": + param_type = param["richards"]["type"] + if param_type != "MvG": raise NotImplementedError("This testing function only supports " "Mualem van Genuchten parameterizations. " "Parameterization was: {}".format( - param["type"])) - return param["parameters"] + param_type)) + return param["richards"]["parameters"] print("Operating on file: {}".format(vtk_file)) diff --git a/python/dorie/dorie/utilities/vtktools/vtkreader.py b/python/dorie/dorie/utilities/vtktools/vtkreader.py index 8ebb6e6e564a6fea5e21406888b7788773cb0248..f3d03992dea106ed2fc83460b3edc2c6e48a8001 100644 --- a/python/dorie/dorie/utilities/vtktools/vtkreader.py +++ b/python/dorie/dorie/utilities/vtktools/vtkreader.py @@ -1,3 +1,4 @@ +from collections.abc import Mapping from vtk import vtkXMLGenericDataObjectReader, \ vtkCellTreeLocator, \ vtkDataSet, \ @@ -5,114 +6,57 @@ from vtk import vtkXMLGenericDataObjectReader, \ reference import numpy as np -class VTKReader: - """A generic VTK reader for any type of VTK XML files. - - It is intended to evaluate cell or vertex (point) data arrays on any point - of the grid. - - Attributes: - _reader: Generic VTK XML data file reader. - _dataset: The complete data set loaded by the reader. - _locator: Fast locator of grid cells based on global coordinates. - _cell: Generic grid cell implementation. - _cell_id: ID of the cell on the grid. - _eval: Evaluation method based on the selected data array type. - _data_array: The actual data to evaluate. - """ - - def __init__(self, file_name, array_name, data_type=None): - """Open the target VTK file, build a cell locator and open the data array. +class VTKDataArray: + """Wrapper for evaluating a VTK data array""" - See :func:`~set_data_array` for details on the arguments regarding the - data array. + def __init__(self, dataset, locator, data_array, data_type): + """Initialize this wrapper from a dataset, a cell locator, + and the underlying VTK data array. Args: - file_name (str): Path to the VTK file. - array_name (str): Name of the data array to evaluate. - data_type (str): Type of the data array. - """ - self._eval = None - self._data_array = None - - # set up the reader - self._reader = vtkXMLGenericDataObjectReader() - self._reader.SetFileName(file_name) - self._reader.Update() - self._dataset = vtkDataSet.SafeDownCast(self._reader.GetOutput()) - - # build the locator - self._locator = vtkCellTreeLocator() - self._locator.SetDataSet(self._dataset) - self._locator.BuildLocator() - - # cache for cell and ID - self._cell = vtkGenericCell() - self._cell_id = -1 # invalid cell - - # set the evaluation data array - self.set_data_array(array_name, data_type) - - def set_data_array (self, array_name, data_type=None): - """Switch the data array (dataset) to evaluate. - - Args: - array_name (str): The name of the data array to evaluate. - data_type (str, optional): The type of the data array (either - ``cell`` or ``vertex``). If this is ``None``, vertex arrays - take precedence. + dataset: VTK dataset object the data array is stored in. + locator: Locator for finding cells based on global coordinates. + data_array: VTK data array containing the data to evaluate. + data_type (str): Type of data stored in the array (``point`` or + ``cell``). Raises: - RuntimeError: Data array was not found. + ValueError: Unsupported data type. """ - self._data_array = None - - def retrieve_array (data, name): - """Return the array if it exists in the data.""" - if data.HasArray(name) == 1: - return data.GetArray(name) - else: - return None - - if data_type == "cell" or data_type is None: - # retrieve cell data array... - self._data_array = retrieve_array(self._dataset.GetCellData(), - array_name) - self._eval = self.__eval_cell__ - if self._data_array is None \ - and (data_type == "vertex" or data_type is None): - # ...or vertex data array - self._data_array = retrieve_array(self._dataset.GetPointData(), - array_name) + self._locator = locator + + self._data_array = data_array + if data_type == "point": self._eval = self.__eval_vertex__ + elif data_type == "cell": + self._eval = self.__eval_cell__ + else: + raise ValueError("Unsupported VTK array data type: {}".format( + data_type + )) - if self._data_array is None: - if data_type is None: - raise RuntimeError("Data array '{}' not found in VTK file" - .format(array_name)) - else: - raise RuntimeError("Data array '{}' of type '{}' not found in " - "VTK file. Supported types are: 'cell', " - "'vertex'".format(array_name, data_type)) + # store physical dataset bounds + self._bounds = np.reshape(dataset.GetBounds(), (3, 2)) - self._num_comp = self._data_array.GetNumberOfComponents() + # cache for cell and ID + self._cell = vtkGenericCell() + self._cell_id = -1 # invalid cell @property def bounds (self): - """Return the outer bounds of the grid. + """Return the outer bounds of the grid this data array lives on. Returns: np.array: The rectangular bounds in any spatial direction. Shape: ``(3, 2)``. """ - bounds = self._dataset.GetBounds() - return np.reshape(bounds, (3, 2)) + return self._bounds @property def number_of_components (self): """Number of components of the evaluated quantity.""" - return self._num_comp + return self._data_array.GetNumberOfComponents() def evaluate (self, position): """Evaluate the selected data array at a certain position. @@ -174,10 +118,106 @@ class VTKReader: def __eval_vertex__ (self, weights): """Evaluation function for a vertex-type data array.""" - data = np.zeros(self._num_comp) + data = np.zeros(self.number_of_components) for point in range(self._cell.GetNumberOfPoints()): point_id = self._cell.GetPointId(point) tup = np.asarray(self._data_array.GetTuple(point_id)) data += weights[point] * tup return data + +class VTKReader (Mapping): + """A generic VTK reader for any type of VTK XML files. + + It maps the data array names to data array evaluators, see + :py:class:`VTKDataArray`. Use it like any Python dict. + + Warning: + VTK files separately store vertex (point) and cell data and associate + the data arrays by their names. If the loaded file contains one array + of each type with the same name, the *vertex (point)* array takes + precedence. **The cell array cannot be accessed!** The key sequence + will contain the name twice, but the associated value will be the + vertex data set in both cases. + + Args: + file_name (str): Name of the VTK file to open. Any file in XML + format is supported. + + Attributes: + _reader: Generic VTK XML data file reader. + _dataset: The complete data set loaded by the reader. + _locator: Fast locator of grid cells based on global coordinates. + _point_data: Pointer to the VTK point data arrays. + _cell_data: Pointer to the VTK cell data arrays. + """ + + def __init__(self, file_name): + """Open the target VTK file, and build a cell locator.""" + + # set up the reader + self._reader = vtkXMLGenericDataObjectReader() + self._reader.SetFileName(file_name) + self._reader.Update() + self._dataset = vtkDataSet.SafeDownCast(self._reader.GetOutput()) + + # build the locator + self._locator = vtkCellTreeLocator() + self._locator.SetDataSet(self._dataset) + self._locator.BuildLocator() + + # ID of the current array + self._point_data = self._dataset.GetPointData() + self._cell_data = self._dataset.GetCellData() + + def __len__(self): + return self._point_data.GetNumberOfArrays() \ + + self._cell_data.GetNumberOfArrays() + + def __getitem__(self, key): + if self._point_data.HasArray(key) == 1: + return VTKDataArray(self._dataset, + self._locator, + self._point_data.GetArray(key), + "point") + elif self._cell_data.HasArray(key) == 1: + return VTKDataArray(self._dataset, + self._locator, + self._cell_data.GetArray(key), + "cell") + else: + raise KeyError("Data array not found: {}".format(key)) + + def __iter__(self): + for i in range(len(self)): + yield self._get_array_name(i) + + def _get_array_name (self, array_id): + """Return the name of the array at the given ID. + + The VTK data set stores separate IDs for cell and point data arrays. + This method subtracts the number of point data arrays from the given + array ID if it is larger than the number of point data arrays, to + properly access cell data arrays. + + Parameters: + array_id (int): Internal ID of the data array. + + Returns: + str: The name of the array associated with the given ID. + """ + if array_id >= len(self): + raise ValueError("Array id '{}' out of bounds of dataset (size " + "is: {})".format(array_id, len(self))) + + num_point_arrays = self._point_data.GetNumberOfArrays() + + if array_id < num_point_arrays: + data = self._point_data + + # start counting from zero again for cell arrays + else: + array_id = array_id - num_point_arrays + data = self._cell_data + + return data.GetArrayName(array_id) diff --git a/python/dorie/setup.py b/python/dorie/setup.py index 59d3a8223ea8c1d93110185a3f9956dfe0d64136..393eb5bfd84ea9c57dc68297908fc0898710752e 100644 --- a/python/dorie/setup.py +++ b/python/dorie/setup.py @@ -4,7 +4,6 @@ from setuptools import setup, find_packages setup(name='dorie', version='0.1', - namespace_packages=['dorie'], description='Python package providing the DORiE parameter scraper', author='Lukas Riedel', author_email='dorieteam@iup.uni-heidelberg.de', diff --git a/python/dorie/wrapper/pf_from_file.py.in b/python/dorie/wrapper/pf_from_file.py.in index 351ef5104de3b1b3313f8fb13bd3ddb91ecc6ead..16c68f4e4d9f1c3fb55c1f502f270d17ab53ccbd 100644 --- a/python/dorie/wrapper/pf_from_file.py.in +++ b/python/dorie/wrapper/pf_from_file.py.in @@ -8,7 +8,8 @@ import argparse import warnings import subprocess import multiprocessing -import configparser + +from dune.testtools.parser import parse_ini_file from dorie.parfield.converter import BaseConverter, \ BinaryConverter, \ @@ -47,14 +48,13 @@ if __name__ == "__main__": parser.add_argument('--debug',help='Display warnings',action='store_true',required=False) args = vars(parser.parse_args()) - cp = configparser.ConfigParser(interpolation=None) - with open(args["param"]) as file: - cp.read_file(file) + # retrieve the inifile tree + inifile = parse_ini_file(args["param"]) params = {} for par in ("converter", "outputFile"): try: - params[par] = cp["general"][par] + params[par] = inifile["general"][par] except KeyError: raise RuntimeError("Missing option general.{} in " "parameter file {}".format(par, args["param"])) @@ -72,7 +72,7 @@ if __name__ == "__main__": print("FFT Field generator failed") sys.exit(1) - rng_outpath = cp.get("general", "tempDir") + rng_outpath = inifile["general"]["tempDir"] input_file = os.path.join(rng_outpath, "field.stoch.h5") # CALL CONVERTER diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt index 19e4ee9995c505a1c0e646c95c0770bc4d2a26c5..d431c8bb821b104fc1ad3e20329991d6c3aaca28 100644 --- a/test/CMakeLists.txt +++ b/test/CMakeLists.txt @@ -45,6 +45,7 @@ add_custom_target(test_run_parallel add_dependencies(test_run_parallel prepare_testing) # dorie exec tests +dorie_add_metaini_test(TARGET dorie METAINI pfg.mini.in) dorie_add_metaini_test(TARGET dorie METAINI run.mini.in) dorie_add_metaini_test(TARGET dorie METAINI plot.mini.in) @@ -64,7 +65,6 @@ dorie_add_metaini_test( METAINI mass_conservation.mini.in SCRIPT dune_execute.py ) -target_link_libraries(${mc_target} spdlog) add_custom_target(test_mass_conservation COMMAND ctest --output-on-failure --tests-regex ^.+mass-conservation.+$ diff --git a/test/bcs/infiltration_3d.dat b/test/bcs/infiltration_3d.dat index 92227163cf31b2d7ffe72e8eb74d79798caa8562..19fbed050eefb4d61df7a611a0d7b871a21ab8a8 100644 --- a/test/bcs/infiltration_3d.dat +++ b/test/bcs/infiltration_3d.dat @@ -11,4 +11,4 @@ spatial_resolution_front_we -1 spatial_resolution_back_sn -1 spatial_resolution_back_we -1 number_BC_change_times 1 -0 neumann -1e-6 dirichlet 0.0 \ No newline at end of file +0 neumann -1e-6 dirichlet 0 \ No newline at end of file diff --git a/test/conforming_flux_jumps.mini.in b/test/conforming_flux_jumps.mini.in index c6e4df6bc1897d463986659e59f50ae2605fe2f7..b60d66435b608768e2dd79004f258c9458e16673 100644 --- a/test/conforming_flux_jumps.mini.in +++ b/test/conforming_flux_jumps.mini.in @@ -34,7 +34,7 @@ initial.equation = -h boundary.file = "{_asset_path}/bcs/infiltration_2d.dat", "{_asset_path}/bcs/infiltration_3d.dat" | expand dim -parameters.file = "{_asset_path}/param/param.yml" +parameters.file = "{_asset_path}/param/richards_param.yml" time.end = 2E2 time.maxTimestep = 2E2 diff --git a/test/const_solute.mini.in b/test/const_solute.mini.in index 31a511d9b13cd0aeebc1939e9bafcdc5c5d77483..c3d54c62b83d627ad91203d054e1fc0bffc68f52 100644 --- a/test/const_solute.mini.in +++ b/test/const_solute.mini.in @@ -27,7 +27,7 @@ output.policy = none boundary.file = "{_asset_path}/bcs/infiltration_2d.dat" -parameters.file = "{_asset_path}/param/param.yml" +parameters.file = "{_asset_path}/param/richards_param.yml" time.end = 1E6 time.maxTimestep = 1E6 @@ -40,9 +40,12 @@ output.fileName = {__name} output.outputPath = {__name} output.vertexData = false output.policy = endOfRichardsStep +output.writeDispersionTensor = true boundary.file = "{_asset_path}/bcs/solute_2d_const.dat" +parameters.file = "{_asset_path}/param/transport_param.yml" + parameters.effHydrDips = 0, 2E-9 | expand diff time.end = 1E6 diff --git a/test/maps/create_param_maps.py b/test/maps/create_param_maps.py index 9dc551b8846a1ea00cee84a3314c9689a63c9ac9..9bbcbe8107cba8ae946ac837c0248f8988344de1 100644 --- a/test/maps/create_param_maps.py +++ b/test/maps/create_param_maps.py @@ -35,7 +35,8 @@ def create_and_write_datasets(args): ode_layered_20 = np.zeros((20, 1), dtype=np.int_) ode_layered_20[10:, ...] = 1 - param_test = np.random.randint(2, size=(50, 50)) + param_test_2 = np.random.randint(2, size=(50, 50)) + param_test_5 = np.random.randint(2, size=(50, 50)) # have a checkerboard! parallel_reference_2d = np.kron([[1, 0] * 2, [0, 1] * 2] * 2, np.ones((10, 10), dtype=np.int_)) @@ -58,7 +59,8 @@ def create_and_write_datasets(args): f.create_dataset("ode_layered_160", data=ode_layered_160) f.create_dataset("ode_layered_40", data=ode_layered_40) f.create_dataset("ode_layered_20", data=ode_layered_20) - f.create_dataset("ut/param_test", data=param_test) + f.create_dataset("ut/param_test_2", data=param_test_2) + f.create_dataset("ut/param_test_5", data=param_test_5) f.create_dataset("parallel_reference_2d", data=parallel_reference_2d) f.create_dataset("parallel_reference_3d", data=parallel_reference_3d) f.close() diff --git a/test/mass_conservation.mini.in b/test/mass_conservation.mini.in index eb889629886be5dc42297f0630bf71f9c3b8caad..1b1979b477b4efce806102725ebb0790270850bc 100644 --- a/test/mass_conservation.mini.in +++ b/test/mass_conservation.mini.in @@ -33,7 +33,7 @@ output.outputPath = mass_conservation | unique boundary.file = "{_asset_path}/bcs/mass_conservation_2d.dat" -parameters.file = "{_asset_path}/param/param.yml" +parameters.file = "{_asset_path}/param/richards_param.yml" time.end = 1E4 time.maxTimestep = 1E4 diff --git a/test/muphi.mini.in b/test/muphi.mini.in index e666c7e666c1d662ea5d516bc786d4c21cff41b2..ef3eed5cca7cfd66e42ecb56778692a1682f4cae 100644 --- a/test/muphi.mini.in +++ b/test/muphi.mini.in @@ -26,7 +26,7 @@ output.vertexData = false boundary.file = "{_asset_path}/bcs/infiltration_2d.dat" -parameters.file = "{_asset_path}/param/param.yml" +parameters.file = "{_asset_path}/param/richards_param.yml" time.end = 1E4 time.maxTimestep = 1E4 diff --git a/test/ode_homogeneous_sand.mini.in b/test/ode_homogeneous_sand.mini.in index f34bb65fe32ffcad1bd747f0ce7efbc2152609d0..01256d83de424a95c23c03b1ee8c3ad07d304425 100644 --- a/test/ode_homogeneous_sand.mini.in +++ b/test/ode_homogeneous_sand.mini.in @@ -23,7 +23,7 @@ initial.equation = -y boundary.file = "{_asset_path}/bcs/infiltration_2d.dat" -parameters.file = "{_asset_path}/param/param.yml" +parameters.file = "{_asset_path}/param/richards_param.yml" time.end = 1E7 time.maxTimestep = 1E7 diff --git a/test/ode_homogeneous_silt.mini.in b/test/ode_homogeneous_silt.mini.in index 034fdcad64139a4714fb2585f47c98e43b673e53..801d66b145c9fffbe1454d239240acae83c1b507 100644 --- a/test/ode_homogeneous_silt.mini.in +++ b/test/ode_homogeneous_silt.mini.in @@ -27,7 +27,7 @@ initial.equation = -y boundary.file = "{_asset_path}/bcs/infiltration_2d.dat" -parameters.file = "{_asset_path}/param/param.yml" +parameters.file = "{_asset_path}/param/richards_param.yml" time.end = 1E7 time.maxTimestep = 1E7 diff --git a/test/ode_layered.mini.in b/test/ode_layered.mini.in index 1f4c5543bd0be894dc68d496709e55f0cab7eeda..6979cf1fbf90d18c9fbad6ac2031eb664803c774 100644 --- a/test/ode_layered.mini.in +++ b/test/ode_layered.mini.in @@ -27,7 +27,7 @@ initial.equation = -y boundary.file = "{_asset_path}/bcs/infiltration_2d.dat" -parameters.file = "{_asset_path}/param/param.yml" +parameters.file = "{_asset_path}/param/richards_param.yml" time.end = 1E7 time.maxTimestep = 1E7 diff --git a/test/parallel_reference.mini.in b/test/parallel_reference.mini.in index 5ff10977b4c9fe639a6a35aaebd00edb13204a40..50a7cb84a5f2ac6224547e9e8b030d714f7bd744 100644 --- a/test/parallel_reference.mini.in +++ b/test/parallel_reference.mini.in @@ -22,7 +22,7 @@ output.verbose = 0 boundary.file = "{_asset_path}/bcs/infiltration_2d.dat" -parameters.file = "{_asset_path}/param/param.yml" +parameters.file = "{_asset_path}/param/richards_param.yml" time.end = 1E5 time.maxTimestep = 1E5 diff --git a/test/parallel_reference_compare.mini.in b/test/parallel_reference_compare.mini.in index 90dc47b654d6c85ddc8cdf1110c942c901da04df..4aa65f5a3295201c199df1d1af0ad350bc104fa2 100644 --- a/test/parallel_reference_compare.mini.in +++ b/test/parallel_reference_compare.mini.in @@ -31,7 +31,7 @@ output.outputPath = parallel_reference_compare | unique boundary.file = "{_asset_path}/bcs/infiltration_2d.dat" -parameters.file = "{_asset_path}/param/param.yml" +parameters.file = "{_asset_path}/param/richards_param.yml" time.end = 1E5 time.maxTimestep = 1E5 diff --git a/test/param/param.yml b/test/param/param.yml deleted file mode 100644 index b95a2c843e3dbd341874661297d8655098fc54e4..0000000000000000000000000000000000000000 --- a/test/param/param.yml +++ /dev/null @@ -1,25 +0,0 @@ -volumes: - sand: - index: 0 - type: MvG - parameters: - alpha: -2.3 - n: 4.17 - k0: 2.2E-5 - theta_r: 0.03 - theta_s: 0.31 - tau: -1.1 - - silt: - index: 1 - type: MvG - parameters: - alpha: -0.7 - n: 1.3 - k0: 1.0E-5 - theta_r: 0.01 - theta_s: 0.41 - tau: 0.0 - -scaling: - type: None diff --git a/test/param/richards_param.yml b/test/param/richards_param.yml new file mode 100644 index 0000000000000000000000000000000000000000..580a49bcf00824b7ed9f1eeffe26a1a48161cffe --- /dev/null +++ b/test/param/richards_param.yml @@ -0,0 +1,27 @@ +volumes: + sand: + index: 0 + richards: + type: MvG + parameters: + alpha: -2.3 + n: 4.17 + k0: 2.2E-5 + theta_r: 0.03 + theta_s: 0.31 + tau: -1.1 + + silt: + index: 1 + richards: + type: MvG + parameters: + alpha: -0.7 + n: 1.3 + k0: 1.0E-5 + theta_r: 0.01 + theta_s: 0.41 + tau: 0.0 + +scaling: + type: None diff --git a/test/param/transport_param.yml b/test/param/transport_param.yml new file mode 100644 index 0000000000000000000000000000000000000000..08821c65334eca0be177bfe9f34aaa4f371da108 --- /dev/null +++ b/test/param/transport_param.yml @@ -0,0 +1,22 @@ +volumes: + sand: + index: 0 + transport: + type: Deff_MQ1 + Dhm_iso + parameters: + char_length: 1.5E-11 + mol_diff: 2.E-9 + phi: 0.34 + lambda_t: 0.0025 + lambda_l: 0.025 + + silt: + index: 1 + transport: + type: Deff_MQ1 + Dhm_iso + parameters: + char_length: 1.5E-11 + mol_diff: 2.E-9 + phi: 0.34 + lambda_t: 0.0025 + lambda_l: 0.025 \ No newline at end of file diff --git a/test/pfg.mini.in b/test/pfg.mini.in new file mode 100644 index 0000000000000000000000000000000000000000..2489983307dff4f20a8e44c11629e4f8c4459fd8 --- /dev/null +++ b/test/pfg.mini.in @@ -0,0 +1,24 @@ +include ${CMAKE_BINARY_DIR}/doc/default_files/parfield.ini + +_converter = none, binary, exponential | expand +_covariance = exponential, gaussian | expand +_dim = 2, 3 | expand dim + +__name = exec_pfg-{_converter}-{_covariance}-{_dim} +_test_command = pfg + +[general] +outputFile = pfg.h5 +dataset = {_converter}-{_covariance}-{_dim} +tempDir = ./pfg-tmp/ +converter = {_converter} + +[grid] +dimensions = {_dim} +extensions = 1 1, 1 1 1 | expand dim +cells = 50 50, 10 10 10 | expand dim + +[stochastic] +seed = 2 +covariance = {_covariance} +corrLength = .1 .05, .2 .1 .05 | expand dim diff --git a/test/run.mini.in b/test/run.mini.in index a1b9051cda8ee911b75af964acf1dfd00cc71728..b5ff28f3128afe5e96213d5727f662548f2374f5 100644 --- a/test/run.mini.in +++ b/test/run.mini.in @@ -17,7 +17,7 @@ output.outputPath = run | unique name boundary.file = "{_asset_path}/bcs/infiltration_2d.dat" -parameters.file = "{_asset_path}/param/param.yml" +parameters.file = "{_asset_path}/param/richards_param.yml" time.end = 1E4 time.maxTimestep = 1E4