diff --git a/dumux/assembly/CMakeLists.txt b/dumux/assembly/CMakeLists.txt
index 65a7539064d99a5446047b668a8db18068f110f2..2543a6631f3a61445e1bbb4aeeb8461b93115bd7 100644
--- a/dumux/assembly/CMakeLists.txt
+++ b/dumux/assembly/CMakeLists.txt
@@ -8,6 +8,7 @@ entitycolor.hh
fvassembler.hh
fvlocalassemblerbase.hh
fvlocalresidual.hh
+initialsolution.hh
jacobianpattern.hh
numericepsilon.hh
partialreassembler.hh
diff --git a/dumux/assembly/initialsolution.hh b/dumux/assembly/initialsolution.hh
new file mode 100644
index 0000000000000000000000000000000000000000..dc8bf7ebf3a04615577efcac8d7c40de4c726cd6
--- /dev/null
+++ b/dumux/assembly/initialsolution.hh
@@ -0,0 +1,111 @@
+// -*- 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 <http://www.gnu.org/licenses/>. *
+ *****************************************************************************/
+/*!
+ * \file
+ * \ingroup Assembly
+ * \brief Function to create initial solution vectors
+ */
+#ifndef DUMUX_ASSEMBLY_INITIAL_SOLUTION_HH
+#define DUMUX_ASSEMBLY_INITIAL_SOLUTION_HH
+
+#include <vector>
+#include <type_traits>
+
+#include <dumux/discretization/method.hh>
+
+namespace Dumux {
+
+/*!
+ * \file
+ * \ingroup Assembly
+ * \brief Set a solution vector to the initial solution provided by the problem
+ */
+template<class SolutionVector, class Problem>
+void assembleInitialSolution(SolutionVector& sol, const Problem& problem)
+{
+ const auto& gg = problem.gridGeometry();
+ using GridGeometry = std::decay_t<decltype(gg)>;
+
+ // box method
+ if constexpr (GridGeometry::discMethod == DiscretizationMethod::box)
+ {
+ constexpr int dim = GridGeometry::GridView::dimension;
+ const auto numDofs = gg.vertexMapper().size();
+ const auto numVert = gg.gridView().size(dim);
+ sol.resize(numDofs);
+
+ // if there are more dofs than vertices (enriched nodal dofs), we have to
+ // call initial for all dofs at the nodes, coming from all neighboring elements.
+ if (numDofs != numVert)
+ {
+ std::vector<bool> dofVisited(numDofs, false);
+ for (const auto& element : elements(gg.gridView()))
+ {
+ for (int i = 0; i < element.subEntities(dim); ++i)
+ {
+ const auto dofIdxGlobal = gg.vertexMapper().subIndex(element, i, dim);
+ // forward to implementation if value at dof is not set yet
+ if (!dofVisited[dofIdxGlobal])
+ {
+ sol[dofIdxGlobal] = problem.initial(element.template subEntity<dim>(i));
+ dofVisited[dofIdxGlobal] = true;
+ }
+ }
+ }
+ }
+
+ // otherwise we directly loop over the vertices
+ else
+ {
+ for (const auto& vertex : vertices(gg.gridView()))
+ sol[gg.vertexMapper().index(vertex)] = problem.initial(vertex);
+ }
+ }
+
+ // staggered methods
+ else if constexpr (GridGeometry::discMethod == DiscretizationMethod::staggered)
+ {
+ problem.applyInitialSolution(sol);
+ }
+
+ // default: cell-centered methods
+ else
+ {
+ sol.resize(gg.numDofs());
+ for (const auto& element : elements(gg.gridView()))
+ sol[gg.elementMapper().index(element)] = problem.initial(element);
+ }
+}
+
+/*!
+ * \file
+ * \ingroup Assembly
+ * \brief Create a solution vector filled with the initial solution provided by the problem
+ */
+template<class SolutionVector, class Problem>
+SolutionVector makeInitialSolution(const Problem& problem)
+{
+ SolutionVector sol;
+ assembleInitialSolution(sol, problem);
+ return sol;
+}
+
+} // end namespace Dumux
+
+#endif
diff --git a/dumux/common/fvproblem.hh b/dumux/common/fvproblem.hh
index 47e6f4313e90b308672d9147b2dce9bf18751da5..ea3f813b07a4344fbe54ccc2679ce58254859020 100644
--- a/dumux/common/fvproblem.hh
+++ b/dumux/common/fvproblem.hh
@@ -34,6 +34,8 @@
#include <dumux/common/parameters.hh>
#include <dumux/discretization/method.hh>
+#include <dumux/assembly/initialsolution.hh>
+
namespace Dumux {
/*!
@@ -65,13 +67,6 @@ class FVProblem
using PointSourceMap = std::map< std::pair<std::size_t, std::size_t>,
std::vector<PointSource> >;
- using GridVariables = GetPropType<TypeTag, Properties::GridVariables>;
- using ElementFluxVariablesCache = typename GetPropType<TypeTag, Properties::GridFluxVariablesCache>::LocalView;
- using ElementVolumeVariables = typename GridVariables::GridVolumeVariables::LocalView;
- using VolumeVariables = typename ElementVolumeVariables::VolumeVariables;
-
- using SolutionVector = GetPropType<TypeTag, Properties::SolutionVector>;
-
static constexpr bool isBox = GridGeometry::discMethod == DiscretizationMethod::box;
static constexpr bool isStaggered = GridGeometry::discMethod == DiscretizationMethod::staggered;
@@ -280,6 +275,7 @@ public:
* Negative values mean influx.
* E.g. for the mass balance that would be the mass flux in \f$ [ kg / (m^2 \cdot s)] \f$.
*/
+ template<class ElementVolumeVariables, class ElementFluxVariablesCache>
NumEqVector neumann(const Element& element,
const FVElementGeometry& fvGeometry,
const ElementVolumeVariables& elemVolVars,
@@ -324,6 +320,7 @@ public:
* that the conserved quantity is created, negative ones mean that it vanishes.
* E.g. for the mass balance that would be a mass rate in \f$ [ kg / (m^3 \cdot s)] \f$.
*/
+ template<class ElementVolumeVariables>
NumEqVector source(const Element &element,
const FVElementGeometry& fvGeometry,
const ElementVolumeVariables& elemVolVars,
@@ -387,6 +384,7 @@ public:
* Positive values mean that the conserved quantity is created, negative ones mean that it vanishes.
* E.g. for the mass balance that would be a mass rate in \f$ [ kg / s ] \f$.
*/
+ template<class ElementVolumeVariables>
void pointSource(PointSource& source,
const Element &element,
const FVElementGeometry& fvGeometry,
@@ -419,7 +417,7 @@ public:
* \brief Add source term derivative to the Jacobian
* \note Only needed in case of analytic differentiation and solution dependent sources
*/
- template<class MatrixBlock>
+ template<class MatrixBlock, class VolumeVariables>
void addSourceDerivatives(MatrixBlock& block,
const Element& element,
const FVElementGeometry& fvGeometry,
@@ -432,6 +430,7 @@ public:
* Caution: Only overload this method in the implementation if you know
* what you are doing.
*/
+ template<class ElementVolumeVariables>
NumEqVector scvPointSources(const Element &element,
const FVElementGeometry& fvGeometry,
const ElementVolumeVariables& elemVolVars,
@@ -508,10 +507,10 @@ public:
* \brief Applies the initial solution for all degrees of freedom of the grid.
* \param sol the initial solution vector
*/
+ template<class SolutionVector>
void applyInitialSolution(SolutionVector& sol) const
{
- // set the initial values by forwarding to a specialized method
- applyInitialSolutionImpl_(sol, std::integral_constant<bool, isBox>());
+ assembleInitialSolution(sol, asImp_());
}
/*!
@@ -596,60 +595,6 @@ protected:
{ return *static_cast<const Implementation *>(this); }
private:
- /*!
- * \brief Applies the initial solution for the box method
- */
- void applyInitialSolutionImpl_(SolutionVector& sol, /*isBox=*/std::true_type) const
- {
- const auto numDofs = gridGeometry_->vertexMapper().size();
- const auto numVert = gridGeometry_->gridView().size(dim);
- sol.resize(numDofs);
-
- // if there are more dofs than vertices (enriched nodal dofs), we have to
- // call initial for all dofs at the nodes, coming from all neighboring elements.
- if (numDofs != numVert)
- {
- std::vector<bool> dofVisited(numDofs, false);
- for (const auto& element : elements(gridGeometry_->gridView()))
- {
- for (int i = 0; i < element.subEntities(dim); ++i)
- {
- const auto dofIdxGlobal = gridGeometry_->vertexMapper().subIndex(element, i, dim);
-
- // forward to implementation if value at dof is not set yet
- if (!dofVisited[dofIdxGlobal])
- {
- sol[dofIdxGlobal] = asImp_().initial(element.template subEntity<dim>(i));
- dofVisited[dofIdxGlobal] = true;
- }
- }
- }
- }
-
- // otherwise we directly loop over the vertices
- else
- {
- for (const auto& vertex : vertices(gridGeometry_->gridView()))
- {
- const auto dofIdxGlobal = gridGeometry_->vertexMapper().index(vertex);
- sol[dofIdxGlobal] = asImp_().initial(vertex);
- }
- }
- }
-
- /*!
- * \brief Applies the initial solution for cell-centered methods
- */
- void applyInitialSolutionImpl_(SolutionVector& sol, /*isBox=*/std::false_type) const
- {
- sol.resize(gridGeometry_->numDofs());
- for (const auto& element : elements(gridGeometry_->gridView()))
- {
- const auto dofIdxGlobal = gridGeometry_->elementMapper().index(element);
- sol[dofIdxGlobal] = asImp_().initial(element);
- }
- }
-
//! The finite volume grid geometry
std::shared_ptr<const GridGeometry> gridGeometry_;