diff --git a/.gitlab-ci/default.yml b/.gitlab-ci/default.yml
index ce9026b5e2bcf759dd393f438625afc0de08dd76..bf70c20fd83af11c96295dc068224cfd19f86b26 100644
--- a/.gitlab-ci/default.yml
+++ b/.gitlab-ci/default.yml
@@ -146,6 +146,7 @@ test cpp:
variables:
OMPI_ALLOW_RUN_AS_ROOT: 1
OMPI_ALLOW_RUN_AS_ROOT_CONFIRM: 1
+ DUMUX_NUM_THREADS: 4
script:
- |
pushd build-cmake
@@ -168,6 +169,7 @@ test python:
variables:
OMPI_ALLOW_RUN_AS_ROOT: 1
OMPI_ALLOW_RUN_AS_ROOT_CONFIRM: 1
+ DUMUX_NUM_THREADS: 4
script:
# restore Python virtual env from cache (job:configure artifacts) (Dune 2.9)
- |
diff --git a/CHANGELOG.md b/CHANGELOG.md
index 2d780782266d70ad754924cf54da32d60ebe1951..838ba967d2d243a34794861b9bcc56bd05c27d24 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -13,6 +13,20 @@ Differences Between DuMux 3.5 and DuMux 3.4
### Improvements and Enhancements
+- __ParallelFor and multithreaded assembly__: For single domain applications using Box of CCTpfa discretization,
+ the possibility for multithreaded (thread-parallel) assembly has been added. It is using the newly added
+ `Dumux::parallelFor` interface. This features requires a backend to be available.
+ This can be either external backends (TBB, Kokkos) or a working C++ STL parallel algorithms setup.
+ The backend will be selected automatically, if found. You can also specify the backend by setting the
+ compiler definition `DUMUX_MULTITHREADING_BACKEND=TBB,Cpp,Kokkos,Serial`, where Serial forces to always run in serial.
+ For the assembly, you can explicitly turn multithreading off setting the runtime parameter `Assembly.Multithreading = false`.
+ If a backend is available and the discretization allows it, the default is multithreaded assembly.
+ When using TBB or Kokkos is it required (recommended) to use `Dumux::initialize` in the main file.
+
+- __initialize__: A new control function `Dumux::initialize` has been added which initializes shared and distributed
+ memory parallelism helpers. It is recommended (and may be required for multithreaded applications) to use `Dumux::initialize`
+ instead of manually initializing the `Dune::MPIHelper`.
+
- __Discretization tags__: We introduced tags in the namespace `DiscretizationMethods` (with s) for each discretization method.
These tags replace the `enum class DiscretizationMethod`. Tags have several advantages over the enum. Each tag is a named type
(see `dumux/common/tag.hh`) so they can for example be used in tag dispatch. Moreover specializing with tags is extensible.
diff --git a/CMakeLists.txt b/CMakeLists.txt
index 45da9af9b5417da7cdbae7ae93b371348baee442..560508e4aadb5974eb1adae2a7091207af67a9f5 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -17,6 +17,7 @@ find_package(dune-common)
list(APPEND CMAKE_MODULE_PATH ${dune-common_MODULE_PATH}
"${PROJECT_SOURCE_DIR}/cmake/modules")
+
#include the dune macros
include(DuneMacros)
diff --git a/cmake/modules/AddKokkosFlags.cmake b/cmake/modules/AddKokkosFlags.cmake
new file mode 100644
index 0000000000000000000000000000000000000000..38714cc17db956516f9072c22b11a55451e4b676
--- /dev/null
+++ b/cmake/modules/AddKokkosFlags.cmake
@@ -0,0 +1,12 @@
+include_guard(GLOBAL)
+
+# set variable for config.h
+set(HAVE_KOKKOS ${Kokkos_FOUND})
+
+# perform DUNE-specific setup tasks
+if (Kokkos_FOUND)
+ dune_register_package_flags(
+ COMPILE_DEFINITIONS ENABLE_KOKKOS=1
+ LIBRARIES Kokkos::kokkos
+ )
+endif()
diff --git a/cmake/modules/AddOpenMPFlags.cmake b/cmake/modules/AddOpenMPFlags.cmake
new file mode 100644
index 0000000000000000000000000000000000000000..c96c92b4c375ab0561af0358cda4280cda9e248c
--- /dev/null
+++ b/cmake/modules/AddOpenMPFlags.cmake
@@ -0,0 +1,12 @@
+include_guard(GLOBAL)
+
+# set variable for config.h
+set(HAVE_OPENMP ${OpenMP_FOUND})
+
+# perform DUNE-specific setup tasks
+if (OpenMP_FOUND)
+ dune_register_package_flags(
+ COMPILE_DEFINITIONS ENABLE_OPENMP=1
+ LIBRARIES OpenMP::OpenMP_CXX
+ )
+endif()
diff --git a/cmake/modules/DumuxMacros.cmake b/cmake/modules/DumuxMacros.cmake
index 80738630ae6f5adfcc87cc8774d63458afe16c5e..e298b91752c0fa2c41c19c6805136dbe5c927ea3 100644
--- a/cmake/modules/DumuxMacros.cmake
+++ b/cmake/modules/DumuxMacros.cmake
@@ -13,3 +13,74 @@ find_package(NLOPT QUIET)
find_package(PTScotch QUIET)
include(AddPTScotchFlags)
find_package(PVPython QUIET)
+
+find_package(Kokkos QUIET)
+include(AddKokkosFlags)
+
+# possibly link against TBB
+# even if an older version is found
+# otherwise we get linker errors
+# beacuse of inconsistencies with
+# dune-common's TBB setup
+find_package(TBB)
+include(AddTBBFlags)
+
+# in a second step make sure the
+# minimum TBB version required is found
+set(DUMUX_MIN_TBB_VERSION 2021)
+if(TBB_FOUND)
+ if(TBB_VERSION_MAJOR VERSION_LESS DUMUX_MIN_TBB_VERSION)
+ find_package(TBB ${DUMUX_MIN_TBB_VERSION})
+ # disable TBB manually if required version not found
+ if(NOT TBB_FOUND)
+ message(STATUS "Disabling TBB since version requirement not satisfied (>= ${DUMUX_MIN_TBB_VERSION}).")
+ set(ENABLE_TBB FALSE)
+ set(HAVE_TBB FALSE)
+ endif()
+ endif()
+endif()
+
+find_package(OpenMP QUIET)
+include(AddOpenMPFlags)
+
+# test if we can use parallel algorithms
+check_cxx_symbol_exists(
+ "std::execution::par_unseq"
+ "execution"
+ DUMUX_HAVE_CXX_EXECUTION_POLICY
+)
+
+# setup multithreading backend
+if(NOT DUMUX_MULTITHREADING_BACKEND)
+ if(TBB_FOUND)
+ set(DUMUX_MULTITHREADING_BACKEND "TBB" CACHE STRING "The multithreading backend")
+ message(STATUS "Dumux multithreading backed: TBB")
+ elseif(OpenMP_FOUND)
+ set(DUMUX_MULTITHREADING_BACKEND "OpenMP" CACHE STRING "The multithreading backend")
+ message(STATUS "Dumux multithreading backed: OpenMP")
+ elseif(Kokkos_FOUND)
+ set(DUMUX_MULTITHREADING_BACKEND "Kokkos" CACHE STRING "The multithreading backend")
+ message(STATUS "Dumux multithreading backed: Kokkos")
+ elseif(DUMUX_HAVE_CXX_EXECUTION_POLICY)
+ set(DUMUX_MULTITHREADING_BACKEND "Cpp" CACHE STRING "The multithreading backend")
+ message(STATUS "Dumux multithreading backed: Cpp")
+ else()
+ set(DUMUX_MULTITHREADING_BACKEND "Serial" CACHE STRING "The multithreading backend")
+ message(STATUS "Dumux multithreading backed: Serial")
+ endif()
+
+# abort if a multithreading backend has been manually selected
+# but it is not available
+else()
+ if(DUMUX_MULTITHREADING_BACKEND STREQUAL "TBB" AND NOT TBB_FOUND)
+ message(FATAL_ERROR "Selected TBB as Dumux multithreading backed but TBB has not been found")
+ elseif(DUMUX_MULTITHREADING_BACKEND STREQUAL "OpenMP" AND NOT OpenMP_FOUND)
+ message(FATAL_ERROR "Selected OpenMP as Dumux multithreading backed but OpenMP has not been found")
+ elseif(DUMUX_MULTITHREADING_BACKEND STREQUAL "Kokkos" AND NOT Kokkos_FOUND)
+ message(FATAL_ERROR "Selected Kokkos as Dumux multithreading backed but Kokkos has not been found")
+ elseif(DUMUX_MULTITHREADING_BACKEND STREQUAL "Cpp" AND NOT DUMUX_HAVE_CXX_EXECUTION_POLICY)
+ message(FATAL_ERROR "Selected Cpp as Dumux multithreading backed but your compiler does not implement parallel STL")
+ else()
+ message(STATUS "Dumux multithreading backed: ${DUMUX_MULTITHREADING_BACKEND}")
+ endif()
+endif()
diff --git a/config.h.cmake b/config.h.cmake
index f7f2bc4b053db7fe4ce5ba41320f5226b48b11db..8df45593a1fb170c5f362bb250ec7e919b9c6a72 100644
--- a/config.h.cmake
+++ b/config.h.cmake
@@ -70,6 +70,18 @@
/* Define to 1 if quadmath was found */
#cmakedefine HAVE_QUAD 1
+/* Set if Kokkos was found */
+#cmakedefine HAVE_KOKKOS ENABLE_KOKKOS
+
+/* Set if OpenMP was found */
+#cmakedefine HAVE_OPENMP ENABLE_OPENMP
+
+/* Set the DUMUX_MULTITHREADING_BACKEND */
+#define DUMUX_MULTITHREADING_BACKEND ${DUMUX_MULTITHREADING_BACKEND}
+
+/* Set HAVE_CPP_PARALLEL_ALGORITHMS if available */
+#cmakedefine HAVE_CPP_PARALLEL_ALGORITHMS DUMUX_HAVE_CXX_EXECUTION_POLICY
+
/* end dumux
Everything below here will be overwritten
*/
diff --git a/dumux/assembly/coloring.hh b/dumux/assembly/coloring.hh
new file mode 100644
index 0000000000000000000000000000000000000000..409d102b28dc447dbe67061110077fdab34e2173
--- /dev/null
+++ b/dumux/assembly/coloring.hh
@@ -0,0 +1,252 @@
+// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
+// vi: set et ts=4 sw=4 sts=4:
+/*****************************************************************************
+ * See the file COPYING for full copying permissions. *
+ * *
+ * This program is free software: you can redistribute it and/or modify *
+ * it under the terms of the GNU General Public License as published by *
+ * the Free Software Foundation, either version 3 of the License, or *
+ * (at your option) any later version. *
+ * *
+ * This program is distributed in the hope that it will be useful, *
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of *
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
+ * GNU General Public License for more details. *
+ * *
+ * You should have received a copy of the GNU General Public License *
+ * along with this program. If not, see . *
+ *****************************************************************************/
+/*!
+ * \file
+ * \ingroup Assembly
+ * \brief Coloring schemes for shared-memory-parallel assembly
+ */
+#ifndef DUMUX_ASSEMBLY_COLORING_HH
+#define DUMUX_ASSEMBLY_COLORING_HH
+
+#include
+#include
+#include
+#include
+
+#include
+#include
+
+#include
+#include
+
+#ifndef DOXYGEN // hide from doxygen
+namespace Dumux::Detail {
+
+//! Compute a map from dof indices to element indices (helper data for coloring algorithm)
+template
+std::vector>
+computeDofToElementMap(const GridGeometry& gg)
+{
+ std::vector> dofToElements;
+
+ if constexpr (GridGeometry::discMethod == DiscretizationMethods::cctpfa)
+ {
+ dofToElements.resize(gg.gridView().size(0));
+ const auto& eMapper = gg.elementMapper();
+ for (const auto& element : elements(gg.gridView()))
+ {
+ const auto eIdx = eMapper.index(element);
+ for (const auto& intersection : intersections(gg.gridView(), element))
+ if (intersection.neighbor())
+ dofToElements[eMapper.index(intersection.outside())].push_back(eIdx);
+ }
+ }
+
+ else if constexpr (GridGeometry::discMethod == DiscretizationMethods::box)
+ {
+ static constexpr int dim = GridGeometry::GridView::dimension;
+ dofToElements.resize(gg.gridView().size(dim));
+ const auto& vMapper = gg.vertexMapper();
+ for (const auto& element : elements(gg.gridView()))
+ {
+ const auto eIdx = gg.elementMapper().index(element);
+ for (int i = 0; i < element.subEntities(dim); i++)
+ dofToElements[vMapper.subIndex(element, i, dim)].push_back(eIdx);
+ }
+ }
+
+ else
+ DUNE_THROW(Dune::NotImplemented,
+ "Missing coloring scheme implementation for this discretization method"
+ );
+
+ return dofToElements;
+}
+
+/*!
+ * \brief Compute the colors of neighboring nodes in the dependency graph
+ *
+ * Neighboring nodes are those elements that manipulate the
+ * same data structures (e.g. system matrix, volvars, flux cache) in the same places
+ *
+ * \param gridGeometry the grid geometry
+ * \param element the element we want to color
+ * \param colors a vector of current colors for each element (not assigned: -1)
+ * \param dofToElement a map from dof indices to element indices
+ * \param neighborColors a vector to add the colors of neighbor nodes to
+ */
+template
+void addNeighborColors(const GridGeometry& gg,
+ const typename GridGeometry::LocalView::Element& element,
+ const std::vector& colors,
+ const DofToElementMap& dofToElement,
+ std::vector& neighborColors)
+{
+ if constexpr (GridGeometry::discMethod == DiscretizationMethods::cctpfa)
+ {
+ // we modify neighbor elements during the assembly
+ // check who else modifies these neighbor elements
+ const auto& eMapper = gg.elementMapper();
+ for (const auto& intersection : intersections(gg.gridView(), element))
+ {
+ if (intersection.neighbor())
+ {
+ // direct face neighbors
+ const auto nIdx = eMapper.index(intersection.outside());
+ neighborColors.push_back(colors[nIdx]);
+
+ // neighbor-neighbors
+ for (const auto nnIdx : dofToElement[eMapper.index(intersection.outside())])
+ neighborColors.push_back(colors[nnIdx]);
+ }
+ }
+ }
+
+ else if constexpr (GridGeometry::discMethod == DiscretizationMethods::box)
+ {
+ // we modify the vertex dofs of our element during the assembly
+ // check who else modifies these vertex dofs
+ const auto& vMapper = gg.vertexMapper();
+ static constexpr int dim = GridGeometry::GridView::dimension;
+ // direct vertex neighbors
+ for (int i = 0; i < element.subEntities(dim); i++)
+ for (auto eIdx : dofToElement[vMapper.subIndex(element, i, dim)])
+ neighborColors.push_back(colors[eIdx]);
+ }
+
+ else
+ DUNE_THROW(Dune::NotImplemented,
+ "Missing coloring scheme implementation for this discretization method"
+ );
+}
+
+/*!
+ * \brief Find the smallest color (integer >= 0) _not_ present in the given list of colors
+ * \param colors list of colors which are already taken
+ * \param notAssigned container to store which colors are not yet taken (is resized as required)
+ */
+int smallestAvailableColor(const std::vector& colors,
+ std::vector& colorUsed)
+{
+ const int numColors = colors.size();
+ colorUsed.assign(numColors, false);
+
+ // The worst case for e.g. numColors=3 is colors={0, 1, 2}
+ // in which case we return 3 as smallest available color
+ // That means, we only track candidates in the (half-open) interval [0, numColors)
+ // Mark candidate colors which are present in colors
+ for (int i = 0; i < numColors; i++)
+ if (colors[i] >= 0 && colors[i] < numColors)
+ colorUsed[colors[i]] = true;
+
+ // return smallest color not in colors
+ for (int i = 0; i < numColors; i++)
+ if (!colorUsed[i])
+ return i;
+
+ return numColors;
+}
+
+} // end namespace Dumux::Detail
+#endif // DOXYGEN
+
+namespace Dumux {
+
+/*!
+ * \brief Compute iterable lists of element seeds partitioned by color
+ *
+ * Splits up the elements of a grid view into partitions such that
+ * all elements in one partition do not modify global data structures
+ * at the same place during assembly. This is used to allow for
+ * lock-free thread-parallel (shared memory) assembly routines.
+ *
+ * Implements a simply greedy graph coloring algorithm:
+ * For each node (element), assign the smallest available color
+ * not used by any of the neighboring nodes (element with conflicting memory access)
+ * The greedy algorithm doesn't necessarily return the smallest
+ * possible number of colors (that's a hard problem) but is fast
+ *
+ * Returns a struct with access to the colors of each element (member colors)
+ * and vector of element seed sets of the same color (member sets)
+ *
+ * \param gridGeometry the grid geometry
+ * \param verbosity the verbosity level
+ */
+template
+auto computeColoring(const GridGeometry& gg, int verbosity = 1)
+{
+ Dune::Timer timer;
+
+ using ElementSeed = typename GridGeometry::GridView::Grid::template Codim<0>::EntitySeed;
+ struct Coloring
+ {
+ using Sets = std::deque>;
+ using Colors = std::vector;
+
+ Coloring(std::size_t size) : sets{}, colors(size, -1) {}
+
+ Sets sets;
+ Colors colors;
+ };
+
+ Coloring coloring(gg.gridView().size(0));
+
+ // pre-reserve some memory for helper arrays to avoid reallocation
+ std::vector neighborColors; neighborColors.reserve(30);
+ std::vector colorUsed; colorUsed.reserve(30);
+
+ // dof to element map to speed up neighbor search
+ const auto dofToElement = Detail::computeDofToElementMap(gg);
+
+ for (const auto& element : elements(gg.gridView()))
+ {
+ // compute neighbor colors based on discretization-dependent stencil
+ neighborColors.clear();
+ Detail::addNeighborColors(gg, element, coloring.colors, dofToElement, neighborColors);
+
+ // find smallest color (positive integer) not in neighborColors
+ const auto color = Detail::smallestAvailableColor(neighborColors, colorUsed);
+
+ // assign color to element
+ coloring.colors[gg.elementMapper().index(element)] = color;
+
+ // add element to the set of elements with the same color
+ if (color < coloring.sets.size())
+ coloring.sets[color].push_back(element.seed());
+ else
+ coloring.sets.push_back(std::vector{ element.seed() });
+ }
+
+ if (verbosity > 0)
+ std::cout << Fmt::format("Colored {} elements with {} colors in {} seconds.\n",
+ gg.gridView().size(0), coloring.sets.size(), timer.elapsed());
+
+ return coloring;
+}
+
+//! Traits specifying if a given discretization tag supports coloring
+template
+struct SupportsColoring : public std::false_type {};
+
+template<> struct SupportsColoring : public std::true_type {};
+template<> struct SupportsColoring : public std::true_type {};
+
+} // end namespace Dumux
+
+#endif
diff --git a/dumux/assembly/fvassembler.hh b/dumux/assembly/fvassembler.hh
index e1a6614341bb27bd5395915b5825a435e12d8857..eb5fa3b0f0b4bd45091f230244d7a83bb7aa06be 100644
--- a/dumux/assembly/fvassembler.hh
+++ b/dumux/assembly/fvassembler.hh
@@ -24,7 +24,10 @@
#ifndef DUMUX_FV_ASSEMBLER_HH
#define DUMUX_FV_ASSEMBLER_HH
+#include
+#include
#include
+#include
#include
@@ -33,8 +36,12 @@
#include
#include
-#include "jacobianpattern.hh"
-#include "diffmethod.hh"
+#include
+#include
+#include
+
+#include
+
#include "boxlocalassembler.hh"
#include "cclocalassembler.hh"
#include "fclocalassembler.hh"
@@ -95,6 +102,7 @@ class FVAssembler
using GridView = typename GridGeo::GridView;
using LocalResidual = GetPropType;
using Element = typename GridView::template Codim<0>::Entity;
+ using ElementSeed = typename GridView::Grid::template Codim<0>::EntitySeed;
using TimeLoop = TimeLoopBase>;
using SolutionVector = GetPropType;
@@ -127,6 +135,8 @@ public:
, isStationaryProblem_(true)
{
static_assert(isImplicit, "Explicit assembler for stationary problem doesn't make sense!");
+ enableMultithreading_ = SupportsColoring::value
+ && getParam("Assembly.Multithreading", true);
}
/*!
@@ -145,7 +155,10 @@ public:
, timeLoop_(timeLoop)
, prevSol_(&prevSol)
, isStationaryProblem_(!timeLoop)
- {}
+ {
+ enableMultithreading_ = SupportsColoring::value
+ && getParam("Assembly.Multithreading", true);
+ }
/*!
* \brief Assembles the global Jacobian of the residual
@@ -258,8 +271,7 @@ public:
else if (jacobian_->buildMode() != JacobianMatrix::BuildMode::random)
DUNE_THROW(Dune::NotImplemented, "Only BCRS matrices with random build mode are supported at the moment");
- setJacobianPattern();
- setResidualSize();
+ resize_();
}
/*!
@@ -272,13 +284,21 @@ public:
jacobian_->setBuildMode(JacobianMatrix::random);
residual_ = std::make_shared();
- setJacobianPattern();
- setResidualSize();
+ resize_();
+ }
+
+ /*!
+ * \brief Resizes jacobian and residual and recomputes colors
+ */
+ void updateAfterGridAdaption()
+ {
+ resize_();
}
/*!
* \brief Resizes the jacobian and sets the jacobian' sparsity pattern.
*/
+ [[deprecated("Use updateAfterGridAdaption. Will be removed after release 3.5.")]]
void setJacobianPattern()
{
// resize the jacobian and the residual
@@ -290,9 +310,13 @@ public:
// export pattern to jacobian
occupationPattern.exportIdx(*jacobian_);
+
+ // maybe recompute colors
+ computeColors_();
}
//! Resizes the residual
+ [[deprecated("Use updateAfterGridAdaption. Will be removed after release 3.5.")]]
void setResidualSize()
{ residual_->resize(numDofs()); }
@@ -375,13 +399,48 @@ public:
}
private:
+ /*!
+ * \brief Resizes the jacobian and sets the jacobian' sparsity pattern.
+ */
+ void setJacobianPattern_()
+ {
+ // resize the jacobian and the residual
+ const auto numDofs = this->numDofs();
+ jacobian_->setSize(numDofs, numDofs);
+
+ // create occupation pattern of the jacobian
+ const auto occupationPattern = getJacobianPattern(gridGeometry());
+
+ // export pattern to jacobian
+ occupationPattern.exportIdx(*jacobian_);
+ }
+
+ //! Resizes the residual
+ void setResidualSize_()
+ { residual_->resize(numDofs()); }
+
+ //! Computes the colors
+ void computeColors_()
+ {
+ if (enableMultithreading_)
+ elementSets_ = computeColoring(gridGeometry()).sets;
+ }
+
+ //! Update with resizing the number of elements (e.g. grid adaption)
+ void resize_()
+ {
+ setJacobianPattern_();
+ setResidualSize_();
+ computeColors_();
+ }
+
// reset the residual vector to 0.0
void resetResidual_()
{
if(!residual_)
{
residual_ = std::make_shared();
- setResidualSize();
+ setResidualSize_();
}
(*residual_) = 0.0;
@@ -395,7 +454,7 @@ private:
{
jacobian_ = std::make_shared();
jacobian_->setBuildMode(JacobianMatrix::random);
- setJacobianPattern();
+ setJacobianPattern_();
}
if (partialReassembler)
@@ -425,9 +484,24 @@ private:
// try assembling using the local assembly function
try
{
- // let the local assembler add the element contributions
- for (const auto& element : elements(gridView()))
- assembleElement(element);
+ if (enableMultithreading_)
+ {
+ // make this element loop run in parallel
+ // for this we have to color the elements so that we don't get
+ // race conditions when writing into the global matrix
+ // each color can be assembled using multiple threads
+ for (const auto& elements : elementSets_)
+ {
+ Dumux::parallelFor(elements.size(), [&](const std::size_t i)
+ {
+ const auto element = gridView().grid().entity(elements[i]);
+ assembleElement(element);
+ });
+ }
+ }
+ else
+ for (const auto& element : elements(gridView()))
+ assembleElement(element);
// if we get here, everything worked well on this process
succeeded = true;
@@ -495,6 +569,10 @@ private:
//! shared pointers to the jacobian matrix and residual
std::shared_ptr jacobian_;
std::shared_ptr residual_;
+
+ //! element sets for parallel assembly
+ bool enableMultithreading_ = false;
+ std::deque> elementSets_;
};
} // namespace Dumux
diff --git a/dumux/common/initialize.hh b/dumux/common/initialize.hh
index c069c26148d40d0a4bc3d19a962b09f754bfa5d0..597b58e038828cdcf89ee91a5adc354f4e1b989a 100644
--- a/dumux/common/initialize.hh
+++ b/dumux/common/initialize.hh
@@ -26,6 +26,64 @@
#include
+#if HAVE_TBB
+#include
+#include
+#include
+#include
+
+#ifndef DOXYGEN
+namespace Dumux::Detail {
+
+class TBBGlobalControl
+{
+public:
+ static oneapi::tbb::global_control& instance(int& argc, char* argv[])
+ {
+ int maxNumThreads = oneapi::tbb::info::default_concurrency();
+ if (const char* dumuxNumThreads = std::getenv("DUMUX_NUM_THREADS"))
+ maxNumThreads = std::max(1, std::stoi(std::string{ dumuxNumThreads }));
+
+ static oneapi::tbb::global_control global_limit(
+ oneapi::tbb::global_control::max_allowed_parallelism, maxNumThreads
+ );
+
+ return global_limit;
+ }
+};
+
+} // namespace Dumux::Detail
+#endif // DOXYGEN
+
+#endif // HAVE_TBB
+
+
+#if HAVE_OPENMP
+#include
+#endif // HAVE_OPENMP
+
+
+#if HAVE_KOKKOS
+#include
+
+#ifndef DOXYGEN
+namespace Dumux::Detail {
+
+class KokkosScopeGuard
+{
+public:
+ static Kokkos::ScopeGuard& instance(int& argc, char* argv[])
+ {
+ static Kokkos::ScopeGuard guard(argc, argv);
+ return guard;
+ }
+};
+
+} // namespace Dumux::Detail
+#endif // DOXYGEN
+
+#endif // HAVE_KOKKOS
+
namespace Dumux {
void initialize(int& argc, char* argv[])
@@ -33,6 +91,23 @@ void initialize(int& argc, char* argv[])
// initialize MPI if available
// otherwise this will create a sequential (fake) helper
Dune::MPIHelper::instance(argc, argv);
+
+#if HAVE_TBB
+ // initialize TBB and keep global control alive
+ Detail::TBBGlobalControl::instance(argc, argv);
+#endif
+
+#if HAVE_OPENMP
+ if (const char* dumuxNumThreads = std::getenv("DUMUX_NUM_THREADS"))
+ omp_set_num_threads(
+ std::max(1, std::stoi(std::string{ dumuxNumThreads }))
+ );
+#endif
+
+#if HAVE_KOKKOS
+ // initialize Kokkos (command line / environmental variable interface)
+ Detail::KokkosScopeGuard::instance(argc, argv);
+#endif
}
} // end namespace Dumux
diff --git a/dumux/parallel/multithreading.hh b/dumux/parallel/multithreading.hh
new file mode 100644
index 0000000000000000000000000000000000000000..1efc82f8beda386001fb33f8cf1662255d335083
--- /dev/null
+++ b/dumux/parallel/multithreading.hh
@@ -0,0 +1,49 @@
+// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
+// vi: set et ts=4 sw=4 sts=4:
+/*****************************************************************************
+ * See the file COPYING for full copying permissions. *
+ * *
+ * This program is free software: you can redistribute it and/or modify *
+ * it under the terms of the GNU General Public License as published by *
+ * the Free Software Foundation, either version 3 of the License, or *
+ * (at your option) any later version. *
+ * *
+ * This program is distributed in the hope that it will be useful, *
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of *
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
+ * GNU General Public License for more details. *
+ * *
+ * You should have received a copy of the GNU General Public License *
+ * along with this program. If not, see . *
+ *****************************************************************************/
+
+/*!
+ * \file
+ * \ingroup Parallel
+ * \brief Multithreading in Dumux
+ */
+#ifndef DUMUX_PARALLEL_MULTITHREADING_HH
+#define DUMUX_PARALLEL_MULTITHREADING_HH
+
+#ifndef DUMUX_MULTITHREADING_BACKEND
+#define DUMUX_MULTITHREADING_BACKEND Serial
+#endif
+
+namespace Dumux::Detail::Multithreading {
+
+namespace ExecutionBackends {
+
+struct Serial {};
+struct Cpp {};
+struct TBB {};
+struct Kokkos {};
+struct OpenMP {};
+
+} // end namespace ExecutionBackends
+
+// set the execution backend type
+using ExecutionBackend = ExecutionBackends::DUMUX_MULTITHREADING_BACKEND;
+
+} // end namespace Dumux::Detail::Multithreading
+
+#endif
diff --git a/dumux/parallel/parallel_for.hh b/dumux/parallel/parallel_for.hh
new file mode 100644
index 0000000000000000000000000000000000000000..b109dbf71bab80a58fcddda7478da8990dee30f8
--- /dev/null
+++ b/dumux/parallel/parallel_for.hh
@@ -0,0 +1,182 @@
+// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
+// vi: set et ts=4 sw=4 sts=4:
+/*****************************************************************************
+ * See the file COPYING for full copying permissions. *
+ * *
+ * This program is free software: you can redistribute it and/or modify *
+ * it under the terms of the GNU General Public License as published by *
+ * the Free Software Foundation, either version 3 of the License, or *
+ * (at your option) any later version. *
+ * *
+ * This program is distributed in the hope that it will be useful, *
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of *
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
+ * GNU General Public License for more details. *
+ * *
+ * You should have received a copy of the GNU General Public License *
+ * along with this program. If not, see . *
+ *****************************************************************************/
+
+/*!
+ * \file
+ * \ingroup Parallel
+ * \brief Parallel for loop (multithreading)
+ */
+
+#ifndef DUMUX_PARALLEL_PARALLEL_FOR_HH
+#define DUMUX_PARALLEL_PARALLEL_FOR_HH
+
+#include
+
+#if HAVE_CPP_PARALLEL_ALGORITHMS
+#include
+#include
+#include
+#endif
+
+#if HAVE_TBB
+#include
+#endif
+
+#if HAVE_KOKKOS
+#include
+#endif
+
+// contents of the detail namespace might change
+// any time without prior notice (do not use directly)
+#ifndef DOXYGEN // hide from doxygen
+namespace Dumux::Detail {
+
+// This should be specialized for different ExecutionBackends
+template
+class ParallelFor;
+
+
+// Serial backend implementation
+template
+class ParallelFor
+{
+public:
+ ParallelFor(const std::size_t count, const FunctorType& functor)
+ : functor_(functor), count_(count) {}
+
+ void execute() const
+ {
+ for (std::size_t i = 0; i < count_; ++i)
+ functor_(i);
+ }
+
+private:
+ FunctorType functor_;
+ std::size_t count_;
+};
+
+#if HAVE_CPP_PARALLEL_ALGORITHMS
+// C++ parallel algorithms backend implementation
+template
+class ParallelFor
+{
+public:
+ ParallelFor(const std::size_t count, const FunctorType& functor)
+ : functor_(functor), range_(count) {}
+
+ void execute() const
+ {
+ std::for_each(std::execution::par_unseq, range_.begin(), range_.end(), functor_);
+ }
+
+private:
+ FunctorType functor_;
+ Dune::IntegralRange range_;
+};
+#endif
+
+
+#if HAVE_TBB
+// TBB backend implementation
+template
+class ParallelFor
+{
+public:
+ ParallelFor(const std::size_t count, const FunctorType& functor)
+ : functor_(functor), count_(count) {}
+
+ void execute() const
+ {
+ tbb::parallel_for(std::size_t{0}, count_, [&](const std::size_t i){ functor_(i); });
+ }
+
+private:
+ FunctorType functor_;
+ std::size_t count_;
+};
+#endif // HAVE_TBB
+
+#if HAVE_KOKKOS
+// Kokkos backend implementation
+template
+class ParallelFor
+{
+public:
+ ParallelFor(const std::size_t count, const FunctorType& functor)
+ : functor_(functor), count_(count) {}
+
+ void execute() const
+ {
+ Kokkos::parallel_for(count_, [&](const std::size_t i){ functor_(i); });
+ }
+
+private:
+ FunctorType functor_;
+ std::size_t count_;
+};
+#endif // HAVE_KOKKOS
+
+
+#if HAVE_OPENMP
+// OpenMP backend implementation
+template
+class ParallelFor
+{
+public:
+ ParallelFor(const std::size_t count, const FunctorType& functor)
+ : functor_(functor), count_(count) {}
+
+ void execute() const
+ {
+ #pragma omp parallel for
+ for (std::size_t i = 0; i < count_; ++i)
+ functor_(i);
+ }
+
+private:
+ FunctorType functor_;
+ std::size_t count_;
+};
+#endif // HAVE_OPENMP
+
+
+} // end namespace Detail
+#endif // DOXYGEN
+
+
+namespace Dumux {
+
+/*!
+ * \file
+ * \ingroup Parallel
+ * \brief A parallel for loop (multithreading)
+ * \param count the number of work tasks to perform
+ * \param functor functor executed for each task (get task number as argument)
+ */
+template
+inline void parallelFor(const std::size_t count, const FunctorType& functor)
+{
+ using ExecutionBackend = Detail::Multithreading::ExecutionBackend;
+ Detail::ParallelFor action(count, functor);
+ action.execute();
+}
+
+} // end namespace Dumux
+
+#endif
diff --git a/examples/1protationsymmetry/doc/main.md b/examples/1protationsymmetry/doc/main.md
index 2d475ca49e09862011e5a0cd66b4f4b8511bf22e..2d323df541f4c67d84253bce49bc4f43925f08ec 100644
--- a/examples/1protationsymmetry/doc/main.md
+++ b/examples/1protationsymmetry/doc/main.md
@@ -30,8 +30,8 @@ and compute the convergence rates.
#include
#include
-#include
+#include
#include // for GetPropType
#include // for getParam
#include // for integrateL2Error
@@ -56,8 +56,8 @@ int main(int argc, char** argv) try
{
using namespace Dumux;
- // We initialize MPI. Finalization is done automatically on exit.
- Dune::MPIHelper::instance(argc, argv);
+ // maybe initialize MPI and/or multithreading backend
+ Dumux::initialize(argc, argv);
// We parse the command line arguments.
Parameters::init(argc, argv);
diff --git a/examples/1protationsymmetry/main.cc b/examples/1protationsymmetry/main.cc
index f75766827fcccc4536f59791b1f63692a490cb72..bc8801cdd758a454e796ba211f6645a0de5433b2 100644
--- a/examples/1protationsymmetry/main.cc
+++ b/examples/1protationsymmetry/main.cc
@@ -27,8 +27,8 @@
#include
#include
-#include
+#include
#include // for GetPropType
#include // for getParam
#include // for integrateL2Error
@@ -51,8 +51,8 @@ int main(int argc, char** argv) try
{
using namespace Dumux;
- // We initialize MPI. Finalization is done automatically on exit.
- Dune::MPIHelper::instance(argc, argv);
+ // maybe initialize MPI and/or multithreading backend
+ Dumux::initialize(argc, argv);
// We parse the command line arguments.
Parameters::init(argc, argv);
diff --git a/examples/1ptracer/doc/main.md b/examples/1ptracer/doc/main.md
index 3d9741c64f1020bfde92f41c00e350df0be3c3ef..986760dad8e553b0861a0e2ce83d30ace0c1494f 100644
--- a/examples/1ptracer/doc/main.md
+++ b/examples/1ptracer/doc/main.md
@@ -26,11 +26,11 @@ The code documentation is structured as follows:
### Included header files
Click to show includes
-These are DUNE helper classes related to parallel computations, time measurements and file I/O
+These time measurements and parallel backend initialization
```cpp
-#include
#include
+#include
```
The following headers include functionality related to property definition or retrieval, as well as
@@ -87,8 +87,8 @@ int main(int argc, char** argv) try
{
using namespace Dumux;
- // The Dune MPIHelper must be instantiated for each program using Dune
- Dune::MPIHelper::instance(argc, argv);
+ // maybe initialize MPI and/or multithreading backend
+ Dumux::initialize(argc, argv);
// parse command line arguments and input file
Parameters::init(argc, argv);
diff --git a/examples/1ptracer/main.cc b/examples/1ptracer/main.cc
index 4faeafcdbf54aaf8a012c11a71d50aba9669a515..94a15e82d62b3e136c5cbc20dc84a2831e07e82d 100644
--- a/examples/1ptracer/main.cc
+++ b/examples/1ptracer/main.cc
@@ -28,9 +28,9 @@
#include
// [[/exclude]]
-// These are DUNE helper classes related to parallel computations, time measurements and file I/O
-#include
+// These time measurements and parallel backend initialization
#include
+#include
// The following headers include functionality related to property definition or retrieval, as well as
// the retrieval of input parameters specified in the input file or via the command line.
@@ -69,8 +69,8 @@ int main(int argc, char** argv) try
{
using namespace Dumux;
- // The Dune MPIHelper must be instantiated for each program using Dune
- Dune::MPIHelper::instance(argc, argv);
+ // maybe initialize MPI and/or multithreading backend
+ Dumux::initialize(argc, argv);
// parse command line arguments and input file
Parameters::init(argc, argv);
diff --git a/examples/2pinfiltration/doc/main.md b/examples/2pinfiltration/doc/main.md
index 88b4655dff7451a33c4084126620d4babff3a083..b9e9d0a4d2ebd1391a7758992521e5d33ead1178 100644
--- a/examples/2pinfiltration/doc/main.md
+++ b/examples/2pinfiltration/doc/main.md
@@ -42,6 +42,7 @@ In Dumux, a property system is used to specify the model. For this, different pr
```cpp
#include
#include
+#include
```
We include the linear solver to be used to solve the linear system and the nonlinear Newton's method
@@ -100,8 +101,9 @@ int main(int argc, char** argv) try
// we define the type tag for this problem
using TypeTag = Properties::TTag::PointSourceExample;
- //We initialize MPI, finalize is done automatically on exit
- const auto& mpiHelper = Dune::MPIHelper::instance(argc, argv);
+ // maybe initialize MPI and/or multithreading backend
+ Dumux::initialize(argc, argv);
+ const auto& mpiHelper = Dune::MPIHelper::instance();
//We parse command line arguments and input file
Parameters::init(argc, argv);
@@ -296,12 +298,10 @@ We start the time loop. In each time step before we start calculating a new solu
{
// We overwrite the old solution with the new (resized & interpolated) one
xOld = x;
- // We tell the assembler to resize the matrix and set pattern
- assembler->setJacobianPattern();
- // We tell the assembler to resize the residual
- assembler->setResidualSize();
- // We initialize the secondary variables to the new (and "new old") solution
+ // We initialize the secondary variables to the new (and "new old") solution
gridVariables->updateAfterGridAdaption(x);
+ // We also need to update the assembler (sizes of residual and Jacobian change)
+ assembler->updateAfterGridAdaption();
// We update the point source map
problem->computePointSourceMap();
}
diff --git a/examples/2pinfiltration/main.cc b/examples/2pinfiltration/main.cc
index 3f590b274546bebd1eba4a735fb43f429cde2cc7..6a182be9bd1182f7f6d5e460d4215d7b56c24334 100644
--- a/examples/2pinfiltration/main.cc
+++ b/examples/2pinfiltration/main.cc
@@ -35,6 +35,7 @@
// In Dumux, a property system is used to specify the model. For this, different properties are defined containing type definitions, values and methods. All properties are declared in the file properties.hh. Additionally, we include the parameter class, which manages the definition of input parameters by a default value, the inputfile or the command line.
#include
#include
+#include
//We include the linear solver to be used to solve the linear system and the nonlinear Newton's method
#include
@@ -72,8 +73,9 @@ int main(int argc, char** argv) try
// we define the type tag for this problem
using TypeTag = Properties::TTag::PointSourceExample;
- //We initialize MPI, finalize is done automatically on exit
- const auto& mpiHelper = Dune::MPIHelper::instance(argc, argv);
+ // maybe initialize MPI and/or multithreading backend
+ Dumux::initialize(argc, argv);
+ const auto& mpiHelper = Dune::MPIHelper::instance();
//We parse command line arguments and input file
Parameters::init(argc, argv);
@@ -239,12 +241,10 @@ int main(int argc, char** argv) try
{
// We overwrite the old solution with the new (resized & interpolated) one
xOld = x;
- // We tell the assembler to resize the matrix and set pattern
- assembler->setJacobianPattern();
- // We tell the assembler to resize the residual
- assembler->setResidualSize();
- // We initialize the secondary variables to the new (and "new old") solution
+ // We initialize the secondary variables to the new (and "new old") solution
gridVariables->updateAfterGridAdaption(x);
+ // We also need to update the assembler (sizes of residual and Jacobian change)
+ assembler->updateAfterGridAdaption();
// We update the point source map
problem->computePointSourceMap();
}
diff --git a/examples/biomineralization/doc/fluidmaterial.md b/examples/biomineralization/doc/fluidmaterial.md
index a494709af8925e863207619f4a83d314daeaafa6..1324bb0f64017b4b98fd9d752f41fc57aea307e1 100644
--- a/examples/biomineralization/doc/fluidmaterial.md
+++ b/examples/biomineralization/doc/fluidmaterial.md
@@ -62,7 +62,7 @@ The real work (creating the tables) is done by some external program by Span and
```cpp
#include
-#include
+#include
namespace Dumux::ICP {
#include "co2valueslaboratory.inc"
}// end namespace Dumux::ICP
@@ -154,7 +154,6 @@ which are needed to model biomineralization.
#include
#include
#include
-#include
#include
#include
#include
diff --git a/examples/biomineralization/doc/mainfile.md b/examples/biomineralization/doc/mainfile.md
index 165c57104573c0b925fe6afe9479e5af5e92a109..017295283224685419dd03c440ab388e4a6df190 100644
--- a/examples/biomineralization/doc/mainfile.md
+++ b/examples/biomineralization/doc/mainfile.md
@@ -42,6 +42,7 @@ the retrieval of input parameters specified in the input file or via the command
```cpp
#include
#include
+#include
```
The follwoing files contain the nonlinear Newtown method, the linear solver and the assembler
@@ -87,8 +88,9 @@ int main(int argc, char** argv) try
{
using namespace Dumux;
- // initialize MPI, finalize is done automatically on exit
- const auto& mpiHelper = Dune::MPIHelper::instance(argc, argv);
+ // maybe initialize MPI and/or multithreading backend
+ Dumux::initialize(argc, argv);
+ const auto& mpiHelper = Dune::MPIHelper::instance();
// parse command line arguments and input file
Parameters::init(argc, argv);
diff --git a/examples/biomineralization/main.cc b/examples/biomineralization/main.cc
index 45fa5a95e6cde4bf857163572cdb15e36ab7e2b1..8bccfdb60b5d2fbf5463edb64625324b351b72d5 100644
--- a/examples/biomineralization/main.cc
+++ b/examples/biomineralization/main.cc
@@ -37,6 +37,7 @@
// the retrieval of input parameters specified in the input file or via the command line.
#include
#include
+#include
// The follwoing files contain the nonlinear Newtown method, the linear solver and the assembler
#include
@@ -68,8 +69,9 @@ int main(int argc, char** argv) try
{
using namespace Dumux;
- // initialize MPI, finalize is done automatically on exit
- const auto& mpiHelper = Dune::MPIHelper::instance(argc, argv);
+ // maybe initialize MPI and/or multithreading backend
+ Dumux::initialize(argc, argv);
+ const auto& mpiHelper = Dune::MPIHelper::instance();
// parse command line arguments and input file
Parameters::init(argc, argv);
diff --git a/examples/biomineralization/material/co2tableslaboratory.hh b/examples/biomineralization/material/co2tableslaboratory.hh
index f0ad811d3556273bb8fdd3dbff46ac8748ca7848..972dc17ffb4286a038760b862898f3d1cd0cde3e 100644
--- a/examples/biomineralization/material/co2tableslaboratory.hh
+++ b/examples/biomineralization/material/co2tableslaboratory.hh
@@ -30,7 +30,7 @@
// [[codeblock]]
#include
-#include
+#include
namespace Dumux::ICP {
#include "co2valueslaboratory.inc"
}// end namespace Dumux::ICP
diff --git a/examples/biomineralization/material/fluidsystems/biominsimplechemistry.hh b/examples/biomineralization/material/fluidsystems/biominsimplechemistry.hh
index 006f1c7e1ad9d242823bf1930146fbe988dda04d..ee407715cb00975a772acd7467384a0920203f51 100644
--- a/examples/biomineralization/material/fluidsystems/biominsimplechemistry.hh
+++ b/examples/biomineralization/material/fluidsystems/biominsimplechemistry.hh
@@ -42,7 +42,6 @@
#include
#include
#include
-#include
#include
#include
#include
diff --git a/examples/freeflowchannel/README.md b/examples/freeflowchannel/README.md
index 00c7fce8980b48aad9e61888912fa6f57794aa0a..e361dd3a4c3c21720c97a093f8ea5c33e2c770e0 100644
--- a/examples/freeflowchannel/README.md
+++ b/examples/freeflowchannel/README.md
@@ -375,6 +375,7 @@ the retrieval of input parameters specified in the input file or via the command
```cpp
#include
#include
+#include
```
The following files contain the non-linear Newton solver, the available linear solver backends and the assembler for the linear
@@ -430,8 +431,9 @@ int main(int argc, char** argv) try
{
using namespace Dumux;
- // The Dune MPIHelper must be instantiated for each program using Dune
- const auto& mpiHelper = Dune::MPIHelper::instance(argc, argv);
+ // maybe initialize MPI and/or multithreading backend
+ Dumux::initialize(argc, argv);
+ const auto& mpiHelper = Dune::MPIHelper::instance();
// parse command line arguments and input file
Parameters::init(argc, argv);
diff --git a/examples/freeflowchannel/main.cc b/examples/freeflowchannel/main.cc
index abe4b4e163b70ecf368d6ef6acfd126ac9aa8f7a..eda6cc688f8857ca07a87186daa15195962960a1 100644
--- a/examples/freeflowchannel/main.cc
+++ b/examples/freeflowchannel/main.cc
@@ -37,6 +37,7 @@
// the retrieval of input parameters specified in the input file or via the command line.
#include
#include
+#include
// The following files contain the non-linear Newton solver, the available linear solver backends and the assembler for the linear
// systems arising from the staggered-grid discretization.
@@ -75,8 +76,9 @@ int main(int argc, char** argv) try
{
using namespace Dumux;
- // The Dune MPIHelper must be instantiated for each program using Dune
- const auto& mpiHelper = Dune::MPIHelper::instance(argc, argv);
+ // maybe initialize MPI and/or multithreading backend
+ Dumux::initialize(argc, argv);
+ const auto& mpiHelper = Dune::MPIHelper::instance();
// parse command line arguments and input file
Parameters::init(argc, argv);
diff --git a/examples/liddrivencavity/doc/problem.md b/examples/liddrivencavity/doc/problem.md
index d8f4725da4de7f2713de27aeec68f777d9faba56..a76c6fcad560d22d40b0a471fd2d8e4611f5770e 100644
--- a/examples/liddrivencavity/doc/problem.md
+++ b/examples/liddrivencavity/doc/problem.md
@@ -158,7 +158,6 @@ conditions for the Navier-Stokes single-phase flow simulation.
```cpp
#include
#include
-#include
```
Include the `NavierStokesProblem` class, the base
@@ -191,8 +190,10 @@ class LidDrivenCavityExampleProblem : public NavierStokesProblem
using SubControlVolume = typename GridGeometry::SubControlVolume;
using SubControlVolumeFace = typename GridGeometry::SubControlVolumeFace;
using Indices = typename GetPropType::Indices;
- using PrimaryVariables = typename ParentType::PrimaryVariables;
- using NumEqVector = typename ParentType::NumEqVector;
+ using InitialValues = typename ParentType::InitialValues;
+ using Sources = typename ParentType::Sources;
+ using DirichletValues = typename ParentType::DirichletValues;
+ using BoundaryFluxes = typename ParentType::BoundaryFluxes;
using Scalar = GetPropType;
static constexpr auto dimWorld = GridGeometry::GridView::dimensionworld;
@@ -235,9 +236,9 @@ The following function specifies the __values on Dirichlet boundaries__.
We need to define values for the primary variables (velocity).
```cpp
- PrimaryVariables dirichletAtPos(const GlobalPosition &globalPos) const
+ DirichletValues dirichletAtPos(const GlobalPosition &globalPos) const
{
- PrimaryVariables values(0.0);
+ DirichletValues values(0.0);
if constexpr (ParentType::isMomentumProblem())
{
@@ -254,13 +255,13 @@ We define a (zero) mass flux here.
```cpp
template
- NumEqVector neumann(const Element& element,
- const FVElementGeometry& fvGeometry,
- const ElementVolumeVariables& elemVolVars,
- const ElementFluxVariablesCache& elemFluxVarsCache,
- const SubControlVolumeFace& scvf) const
+ BoundaryFluxes neumann(const Element& element,
+ const FVElementGeometry& fvGeometry,
+ const ElementVolumeVariables& elemVolVars,
+ const ElementFluxVariablesCache& elemFluxVarsCache,
+ const SubControlVolumeFace& scvf) const
{
- NumEqVector values(0.0);
+ BoundaryFluxes values(0.0);
if constexpr (!ParentType::isMomentumProblem())
{
@@ -286,9 +287,9 @@ constraint for pressure__ in a single cell.
{ return !ParentType::isMomentumProblem(); }
// Set a fixed pressure a the lower-left cell.
- std::bitset hasInternalDirichletConstraint(const Element& element, const SubControlVolume& scv) const
+ std::bitset hasInternalDirichletConstraint(const Element& element, const SubControlVolume& scv) const
{
- std::bitset values;
+ std::bitset values;
if constexpr (!ParentType::isMomentumProblem())
{
@@ -301,8 +302,8 @@ constraint for pressure__ in a single cell.
}
// Specify the pressure value in the internal Dirichlet cell.
- PrimaryVariables internalDirichlet(const Element& element, const SubControlVolume& scv) const
- { return PrimaryVariables(1.1e5); }
+ DirichletValues internalDirichlet(const Element& element, const SubControlVolume& scv) const
+ { return DirichletValues(1.1e5); }
```
Setting a __reference pressure__ can help to improve the Newton convergence rate by making the numerical derivatives more exact.
@@ -318,9 +319,9 @@ This is related to floating point arithmetic as pressure values are usually much
The following function defines the initial conditions.
```cpp
- PrimaryVariables initialAtPos(const GlobalPosition &globalPos) const
+ InitialValues initialAtPos(const GlobalPosition &globalPos) const
{
- PrimaryVariables values(0.0);
+ InitialValues values(0.0);
if constexpr (!ParentType::isMomentumProblem())
values[Indices::pressureIdx] = 1.0e+5;
@@ -365,6 +366,7 @@ the retrieval of input parameters specified in the input file or via the command
```cpp
#include
#include
+#include
```
The following files contain the multi-domain Newton solver, the available linear solver backends and the assembler for the linear
@@ -449,8 +451,9 @@ int main(int argc, char** argv)
{
using namespace Dumux;
- // The Dune MPIHelper must be instantiated for each program using Dune, it is finalized automatically on exit
- const auto& mpiHelper = Dune::MPIHelper::instance(argc, argv);
+ // maybe initialize MPI and/or multithreading backend
+ Dumux::initialize(argc, argv);
+ const auto& mpiHelper = Dune::MPIHelper::instance();
// parse command line arguments and input file
Parameters::init(argc, argv);
diff --git a/examples/liddrivencavity/main.cc b/examples/liddrivencavity/main.cc
index b0e1d8ba2859a48615c4d3756a965e955e7dd1d7..9e99893bf10f884bd1d2a9d5cf4665d9294412fb 100644
--- a/examples/liddrivencavity/main.cc
+++ b/examples/liddrivencavity/main.cc
@@ -37,6 +37,7 @@
// the retrieval of input parameters specified in the input file or via the command line.
#include
#include
+#include
// The following files contain the multi-domain Newton solver, the available linear solver backends and the assembler for the linear
// systems arising from the staggered-grid discretization.
@@ -106,8 +107,9 @@ int main(int argc, char** argv)
{
using namespace Dumux;
- // The Dune MPIHelper must be instantiated for each program using Dune, it is finalized automatically on exit
- const auto& mpiHelper = Dune::MPIHelper::instance(argc, argv);
+ // maybe initialize MPI and/or multithreading backend
+ Dumux::initialize(argc, argv);
+ const auto& mpiHelper = Dune::MPIHelper::instance();
// parse command line arguments and input file
Parameters::init(argc, argv);
diff --git a/examples/porenetwork_upscaling/doc/main.md b/examples/porenetwork_upscaling/doc/main.md
index 4433880183bc7c5797fec360cdddf27f36ad13a7..7c85039daa388ce1d3efb57681d2ec576d0ce455 100644
--- a/examples/porenetwork_upscaling/doc/main.md
+++ b/examples/porenetwork_upscaling/doc/main.md
@@ -34,6 +34,7 @@ Pore-Network-Model to evaluate the upscaled Darcy permeability of a given networ
#include // for GetPropType
#include // for getParam
+#include
#include // for ILU0BiCGSTABBackend
#include // for LinearPDESolver
@@ -60,6 +61,9 @@ int main(int argc, char** argv) try
{
using namespace Dumux;
+ // maybe initialize MPI and/or multithreading backend
+ Dumux::initialize(argc, argv);
+
// We parse the command line arguments.
Parameters::init(argc, argv);
diff --git a/examples/porenetwork_upscaling/doc/problem.md b/examples/porenetwork_upscaling/doc/problem.md
index 1c4289a2de558d498ce44906a1ab1b092f662bf2..7a5514983565fb2146b2ecdd2a5a9fb743550727 100644
--- a/examples/porenetwork_upscaling/doc/problem.md
+++ b/examples/porenetwork_upscaling/doc/problem.md
@@ -183,7 +183,6 @@ class UpscalingProblem : public PorousMediumFlowProblem
```