From 23974f94bab1ef1be8ee71e9eea7d280aec43df9 Mon Sep 17 00:00:00 2001 From: Timo Koch Date: Tue, 14 Apr 2020 16:24:15 +0200 Subject: [PATCH 1/3] [fvproblem] Move applyInitialSolution to a separate header --- dumux/assembly/CMakeLists.txt | 1 + dumux/assembly/initialsolution.hh | 111 ++++++++++++++++++++++++++++++ dumux/common/fvproblem.hh | 59 +--------------- 3 files changed, 115 insertions(+), 56 deletions(-) create mode 100644 dumux/assembly/initialsolution.hh diff --git a/dumux/assembly/CMakeLists.txt b/dumux/assembly/CMakeLists.txt index 65a7539064..2543a6631f 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 0000000000..dc8bf7ebf3 --- /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 . * + *****************************************************************************/ +/*! + * \file + * \ingroup Assembly + * \brief Function to create initial solution vectors + */ +#ifndef DUMUX_ASSEMBLY_INITIAL_SOLUTION_HH +#define DUMUX_ASSEMBLY_INITIAL_SOLUTION_HH + +#include +#include + +#include + +namespace Dumux { + +/*! + * \file + * \ingroup Assembly + * \brief Set a solution vector to the initial solution provided by the problem + */ +template +void assembleInitialSolution(SolutionVector& sol, const Problem& problem) +{ + const auto& gg = problem.gridGeometry(); + using GridGeometry = std::decay_t; + + // 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 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(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 +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 47e6f4313e..49dea16866 100644 --- a/dumux/common/fvproblem.hh +++ b/dumux/common/fvproblem.hh @@ -34,6 +34,8 @@ #include #include +#include + namespace Dumux { /*! @@ -510,8 +512,7 @@ public: */ void applyInitialSolution(SolutionVector& sol) const { - // set the initial values by forwarding to a specialized method - applyInitialSolutionImpl_(sol, std::integral_constant()); + assembleInitialSolution(sol, asImp_()); } /*! @@ -596,60 +597,6 @@ protected: { return *static_cast(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 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(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 gridGeometry_; -- GitLab From 17182ac16e25d07169eabeaabae66dfe73a2f39f Mon Sep 17 00:00:00 2001 From: Timo Koch Date: Tue, 14 Apr 2020 16:25:08 +0200 Subject: [PATCH 2/3] [fvproblem] make applyInitialSolution a template --- dumux/common/fvproblem.hh | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/dumux/common/fvproblem.hh b/dumux/common/fvproblem.hh index 49dea16866..6573ecb818 100644 --- a/dumux/common/fvproblem.hh +++ b/dumux/common/fvproblem.hh @@ -72,8 +72,6 @@ class FVProblem using ElementVolumeVariables = typename GridVariables::GridVolumeVariables::LocalView; using VolumeVariables = typename ElementVolumeVariables::VolumeVariables; - using SolutionVector = GetPropType; - static constexpr bool isBox = GridGeometry::discMethod == DiscretizationMethod::box; static constexpr bool isStaggered = GridGeometry::discMethod == DiscretizationMethod::staggered; @@ -510,6 +508,7 @@ public: * \brief Applies the initial solution for all degrees of freedom of the grid. * \param sol the initial solution vector */ + template void applyInitialSolution(SolutionVector& sol) const { assembleInitialSolution(sol, asImp_()); -- GitLab From e396a68ba2819aaf76b97651cab2bd8b1fc65bfd Mon Sep 17 00:00:00 2001 From: Timo Koch Date: Tue, 14 Apr 2020 17:00:55 +0200 Subject: [PATCH 3/3] [fvproblem] Make some functions templates to get rid of property dependencies --- dumux/common/fvproblem.hh | 11 +++++------ 1 file changed, 5 insertions(+), 6 deletions(-) diff --git a/dumux/common/fvproblem.hh b/dumux/common/fvproblem.hh index 6573ecb818..ea3f813b07 100644 --- a/dumux/common/fvproblem.hh +++ b/dumux/common/fvproblem.hh @@ -67,11 +67,6 @@ class FVProblem using PointSourceMap = std::map< std::pair, std::vector >; - using GridVariables = GetPropType; - using ElementFluxVariablesCache = typename GetPropType::LocalView; - using ElementVolumeVariables = typename GridVariables::GridVolumeVariables::LocalView; - using VolumeVariables = typename ElementVolumeVariables::VolumeVariables; - 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 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 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 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 + template 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 NumEqVector scvPointSources(const Element &element, const FVElementGeometry& fvGeometry, const ElementVolumeVariables& elemVolVars, -- GitLab