From cd4765e88104dbf6063a9f5d7886af1f82a15656 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Dennis=20Gl=C3=A4ser?= Date: Thu, 4 Feb 2021 21:30:49 +0100 Subject: [PATCH 01/32] introduce time level class --- dumux/timestepping/CMakeLists.txt | 3 ++ dumux/timestepping/timelevel.hh | 78 +++++++++++++++++++++++++++++++ 2 files changed, 81 insertions(+) create mode 100644 dumux/timestepping/CMakeLists.txt create mode 100644 dumux/timestepping/timelevel.hh diff --git a/dumux/timestepping/CMakeLists.txt b/dumux/timestepping/CMakeLists.txt new file mode 100644 index 0000000000..3e7bc52b16 --- /dev/null +++ b/dumux/timestepping/CMakeLists.txt @@ -0,0 +1,3 @@ +install(FILES +timelevel.hh +DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/dumux/timestepping) diff --git a/dumux/timestepping/timelevel.hh b/dumux/timestepping/timelevel.hh new file mode 100644 index 0000000000..cf9e1fd375 --- /dev/null +++ b/dumux/timestepping/timelevel.hh @@ -0,0 +1,78 @@ +// -*- 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 + * \brief Class that represents a time level during time integration. + */ +#ifndef DUMUX_TIME_LEVEL_HH +#define DUMUX_TIME_LEVEL_HH + +namespace Dumux::Experimental { + +/*! + * \brief Class that represents a time level during time integration. + */ +template +class TimeLevel +{ +public: + + /*! + * \brief Construct a time level with a time. + * \note This can be used in contexts outside of time integration, + * where no information on a previous time or time step size is needed. + */ + explicit TimeLevel(Scalar curTime) + : curTime_(curTime) + , prevTime_(curTime) + , timeStepFraction_(1.0) + {} + + /*! + * \brief Construct a time level with information on an ongoing time step. + * \param curTime The current time level + * \param prevTime The previous time level + * \param dtFraction The fraction of a time step this level corresponds to. + * \note Within a time integration step, several time levels might occur + * when multi-stage methods are used. The argument dtFraction allows + * for determining the time that will be reached at the end of the + * time integration step. + */ + explicit TimeLevel(Scalar curTime, Scalar prevTime, Scalar dtFraction) + : curTime_(curTime) + , prevTime_(prevTime) + , timeStepFraction_(dtFraction) + {} + + //! Return the current time + Scalar current() const { return curTime_; } + //! Return the time at the beginning of time integration + Scalar previous() const { return prevTime_; } + //! Return the fraction of the time step this level corresponds to + Scalar timeStepFraction() const { return timeStepFraction_; } + +private: + Scalar curTime_; + Scalar prevTime_; + Scalar timeStepFraction_; +}; + +} // end namespace Dumux::Experimental + +#endif -- GitLab From 1fe2c061a9158ce06217aa6fc3738cc83a538991 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Dennis=20Gl=C3=A4ser?= Date: Thu, 4 Feb 2021 21:33:18 +0100 Subject: [PATCH 02/32] [common] add base clase for (grid)variables --- dumux/common/CMakeLists.txt | 1 + dumux/common/variables.hh | 137 ++++++++++++++++++++++++++++++++++++ 2 files changed, 138 insertions(+) create mode 100644 dumux/common/variables.hh diff --git a/dumux/common/CMakeLists.txt b/dumux/common/CMakeLists.txt index 135b7f0edc..9d03d03b99 100644 --- a/dumux/common/CMakeLists.txt +++ b/dumux/common/CMakeLists.txt @@ -50,4 +50,5 @@ timeloop.hh timemanager.hh valgrind.hh variablelengthspline_.hh +variables.hh DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/dumux/common) diff --git a/dumux/common/variables.hh b/dumux/common/variables.hh new file mode 100644 index 0000000000..d37443183b --- /dev/null +++ b/dumux/common/variables.hh @@ -0,0 +1,137 @@ +// -*- 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 Common + * \copydoc Dumux::Variables + */ +#ifndef DUMUX_VARIABLES_HH +#define DUMUX_VARIABLES_HH + +#include + +#include +#include + +namespace Dumux::Experimental { + +/*! + * \ingroup Discretization + * \brief Class that represents the variables of a model. + * We assume that models are formulated on the basis of primary and + * possibly secondary variables, where the latter may non-linearly + * depend on the former. Variables in Dumux represent the state of + * a numerical solution of a model, consisting of all primary/secondary + * variables and, if the a transient problem is modeled, of time information. + * + * This class defines the interface that is expected of variable classes, + * and it provides the implementation for models that do not require storing + * any additional information besides the primary variables and (optionally) + * time. + * \tparam X The type used for solution vectors, i.e. all primary variables. + */ +template +class Variables +{ + template + struct ScalarExtractorHelper; + + template + struct ScalarExtractorHelper + { using Type = T; }; + + template + struct ScalarExtractorHelper + { + private: + using ValueType = std::decay_t()[0])>; + static constexpr bool indexable = IsIndexable::value; + public: + using Type = typename ScalarExtractorHelper::Type; + }; + +public: + //! export the type of solution vector + using SolutionVector = X; + + //! export the underlying scalar type + using Scalar = typename ScalarExtractorHelper::value>::Type; + + //! export the time representation + using TimeLevel = Dumux::Experimental::TimeLevel; + + //! Default constructor + explicit Variables() : x_(), t_(0.0) {} + + //! Construction from a solution + explicit Variables(const SolutionVector& x, + const TimeLevel& t = TimeLevel{0.0}) + : x_(x), t_(t) + {} + + //! Construction from a solution + explicit Variables(SolutionVector&& x, + const TimeLevel& t = TimeLevel{0.0}) + : x_(std::move(x)), t_(t) + {} + + //! Construction from initializer lambda + template), int> = 0> + explicit Variables(const Initializer& initializeSolution, + const TimeLevel& timeLevel = TimeLevel{0.0}) + : t_(timeLevel) + { + initializeSolution(x_); + } + + //! Return the time level + const TimeLevel& timeLevel() const + { return t_; } + + //! Return reference to the solution + const SolutionVector& dofs() const { return x_; } + + //! Non-const access still required for privar switch (TODO: Remove dependency) + SolutionVector& dofs() { return x_; } + + //! Update the state to a new solution + void update(const SolutionVector& x) + { x_ = x; } + + //! Update the time level only + void updateTime(const TimeLevel& t) + { t_ = t; } + + //! Update the state to a new solution & time level + void update(const SolutionVector& x, + const TimeLevel& t) + { + x_ = x; + t_ = t; + } + +private: + SolutionVector x_; + TimeLevel t_; +}; + +} // end namespace Dumux::Experimental + +#endif -- GitLab From ce6b66113cee3c19a9223ec830af9ceb0bb14972 Mon Sep 17 00:00:00 2001 From: "Dennis.Glaeser" Date: Fri, 16 Oct 2020 15:04:07 +0200 Subject: [PATCH 03/32] [disc] introduce base grid variables --- dumux/discretization/CMakeLists.txt | 1 + dumux/discretization/gridvariables.hh | 70 +++++++++++++++++++++++++++ 2 files changed, 71 insertions(+) create mode 100644 dumux/discretization/gridvariables.hh diff --git a/dumux/discretization/CMakeLists.txt b/dumux/discretization/CMakeLists.txt index 6fbdf151ab..3b145a92fb 100644 --- a/dumux/discretization/CMakeLists.txt +++ b/dumux/discretization/CMakeLists.txt @@ -19,6 +19,7 @@ fluxstencil.hh functionspacebasis.hh fvgridvariables.hh fvproperties.hh +gridvariables.hh localview.hh method.hh rotationpolicy.hh diff --git a/dumux/discretization/gridvariables.hh b/dumux/discretization/gridvariables.hh new file mode 100644 index 0000000000..b8d040c465 --- /dev/null +++ b/dumux/discretization/gridvariables.hh @@ -0,0 +1,70 @@ +// -*- 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 Discretization + * \brief Base class for grid variables + */ +#ifndef DUMUX_GRID_VARIABLES_HH +#define DUMUX_GRID_VARIABLES_HH + +#include + +#include + +namespace Dumux::Experimental { + +/*! + * \ingroup Discretization + * \brief Base class for grid variables. + * \tparam GG The grid geometry type + * \tparam X The type used for solution vectors + */ +template +class GridVariables +: public Variables +{ + using ParentType = Variables; + +public: + //! export the grid geometry type + using GridGeometry = GG; + + /*! + * \brief Constructor from a grid geometry. The remaining arguments must + * be valid arguments for the construction of the Variables class. + */ + template + explicit GridVariables(std::shared_ptr gridGeometry, + Args&&... args) + : ParentType(std::forward(args)...) + , gridGeometry_(gridGeometry) + {} + + //! Return a reference to the grid geometry + const GridGeometry& gridGeometry() const + { return *gridGeometry_; } + +private: + std::shared_ptr gridGeometry_; +}; + +} // end namespace Dumux::Experimental + +#endif -- GitLab From b763bf1cf904e6d595102a511d24570b75e7c49f Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Dennis=20Gl=C3=A4ser?= Date: Thu, 4 Feb 2021 21:38:14 +0100 Subject: [PATCH 04/32] [fvgridvars] add implementation for new grid vars concept --- dumux/discretization/fvgridvariables.hh | 206 +++++++++++++++++++++++- 1 file changed, 204 insertions(+), 2 deletions(-) diff --git a/dumux/discretization/fvgridvariables.hh b/dumux/discretization/fvgridvariables.hh index ce693bd283..92656b7d45 100644 --- a/dumux/discretization/fvgridvariables.hh +++ b/dumux/discretization/fvgridvariables.hh @@ -19,14 +19,14 @@ /*! * \file * \ingroup Discretization - * \brief The grid variable class for finite volume schemes storing variables on scv and scvf (volume and flux variables) + * \brief The grid variable class for finite volume schemes, + * storing variables on scv and scvf (volume and flux variables) */ #ifndef DUMUX_FV_GRID_VARIABLES_HH #define DUMUX_FV_GRID_VARIABLES_HH #include #include -#include namespace Dumux { @@ -169,4 +169,206 @@ private: } // end namespace Dumux +////////////////////////////////////////////////////////////// +// Experimental implementation of new grid variables layout // +////////////////////////////////////////////////////////////// + +#include +#include "gridvariables.hh" + +namespace Dumux::Experimental { + + /*! + * \ingroup Discretization + * \brief Finite volume-specific local view on grid variables. + * \tparam GV The grid variables class + */ + template + class FVGridVariablesLocalView + { + using GridGeometry = typename GV::GridGeometry; + using FVElementGeometry = typename GridGeometry::LocalView; + + using GridView = typename GridGeometry::GridView; + using Element = typename GridView::template Codim<0>::Entity; + + using ElementVolumeVariables = typename GV::GridVolumeVariables::LocalView; + using ElementFluxVariablesCache = typename GV::GridFluxVariablesCache::LocalView; + + public: + //! export corresponding grid-wide class + using GridVariables = GV; + + //! Constructor + FVGridVariablesLocalView(const GridVariables& gridVariables) + : gridVariables_(&gridVariables) + , elemVolVars_(gridVariables.gridVolVars()) + , elemFluxVarsCache_(gridVariables.gridFluxVarsCache()) + {} + + /*! + * \brief Bind this local view to a grid element. + * \param element The grid element + * \param fvGeometry Local view on the grid geometry + */ + void bind(const Element& element, + const FVElementGeometry& fvGeometry) + { + const auto& x = gridVariables().dofs(); + elemVolVars_.bind(element, fvGeometry, x); + elemFluxVarsCache_.bind(element, fvGeometry, elemVolVars_); + } + + /*! + * \brief Bind only the volume variables local view to a grid element. + * \param element The grid element + * \param fvGeometry Local view on the grid geometry + */ + void bindElemVolVars(const Element& element, + const FVElementGeometry& fvGeometry) + { + elemVolVars_.bind(element, fvGeometry, gridVariables().dofs()); + + // unbind flux variables cache + elemFluxVarsCache_ = localView(gridVariables().gridFluxVarsCache()); + } + + //! return reference to the elem vol vars + const ElementVolumeVariables& elemVolVars() const { return elemVolVars_; } + ElementVolumeVariables& elemVolVars() { return elemVolVars_; } + + //! return reference to the flux variables cache + const ElementFluxVariablesCache& elemFluxVarsCache() const { return elemFluxVarsCache_; } + ElementFluxVariablesCache& elemFluxVarsCache() { return elemFluxVarsCache_; } + + //! Return reference to the grid variables + const GridVariables& gridVariables() const + { return *gridVariables_; } + + private: + const GridVariables* gridVariables_; + ElementVolumeVariables elemVolVars_; + ElementFluxVariablesCache elemFluxVarsCache_; + }; + + /*! + * \ingroup Discretization + * \brief The grid variable class for finite volume schemes, storing + * variables on scv and scvf (volume and flux variables). + * \tparam GG the type of the grid geometry + * \tparam GVV the type of the grid volume variables + * \tparam GFVC the type of the grid flux variables cache + * \tparam X the type used for solution vectors + * \todo TODO: GG is an obsolete (or redundant) template argument? + */ + template + class FVGridVariables : public GridVariables + { + using ParentType = GridVariables; + using ThisType = FVGridVariables; + + public: + using typename ParentType::SolutionVector; + + //! export type of the finite volume grid geometry + using GridGeometry = GG; + + //! export type of the grid volume variables + using GridVolumeVariables = GVV; + + //! export type of the volume variables + using VolumeVariables = typename GridVolumeVariables::VolumeVariables; + + //! export primary variable type + using PrimaryVariables = typename VolumeVariables::PrimaryVariables; + + //! export cache type for flux variables + using GridFluxVariablesCache = GFVC; + + //! export the local view on this class + using LocalView = FVGridVariablesLocalView; + + /*! + * \brief Constructor + * \param problem The problem to be solved + * \param gridGeometry The geometry of the computational grid + * \note This constructor initializes the solution using the + * initializer function in the given problem, and thus, + * this only compiles if the problem implements it. + */ + template + FVGridVariables(std::shared_ptr problem, + std::shared_ptr gridGeometry) + : ParentType(gridGeometry, [problem] (auto& x) { problem->applyInitialSolution(x); }) + , gridVolVars_(*problem) + , gridFluxVarsCache_(*problem) + {} + + /*! + * \brief Constructor with custom initialization of the solution. + * \param problem The problem to be solved + * \param gridGeometry The geometry of the computational grid + * \param solOrInitializer This can be either a reference to a solution + * vector, or an initializer lambda. + * See Dumux::Experimental::Variables. + */ + template + FVGridVariables(std::shared_ptr problem, + std::shared_ptr gridGeometry, + SolOrInitializer&& solOrInitializer) + : ParentType(gridGeometry, std::forward(solOrInitializer)) + , gridVolVars_(*problem) + , gridFluxVarsCache_(*problem) + { + gridVolVars_.update(this->gridGeometry(), this->dofs()); + gridFluxVarsCache_.update(this->gridGeometry(), gridVolVars_, this->dofs(), true); + } + + //! Update all variables that may be affected by a change in solution + void update(const SolutionVector& curSol) + { + ParentType::update(curSol); + + // resize and update the volVars with the initial solution + gridVolVars_.update(this->gridGeometry(), curSol); + + // update the flux variables caches + gridFluxVarsCache_.update(this->gridGeometry(), gridVolVars_, curSol); + } + + //! Force the update of all variables + void forceUpdateAll(const SolutionVector& curSol) + { + ParentType::update(curSol); + + // resize and update the volVars with the initial solution + gridVolVars_.update(this->gridGeometry(), curSol); + + // update the flux variables caches + gridFluxVarsCache_.update(this->gridGeometry(), gridVolVars_, curSol, true); + } + + //! return the flux variables cache + const GridFluxVariablesCache& gridFluxVarsCache() const + { return gridFluxVarsCache_; } + + //! return the flux variables cache + GridFluxVariablesCache& gridFluxVarsCache() + { return gridFluxVarsCache_; } + + //! return the current volume variables + const GridVolumeVariables& gridVolVars() const + { return gridVolVars_; } + + //! return the current volume variables + GridVolumeVariables& gridVolVars() + { return gridVolVars_; } + + private: + GridVolumeVariables gridVolVars_; //!< the current volume variables (primary and secondary variables) + GridFluxVariablesCache gridFluxVarsCache_; //!< the flux variables cache + }; + +} // end namespace Dumux::Experimental + #endif -- GitLab From 184509c4a4acc846ebd084cdc37e752893bb314d Mon Sep 17 00:00:00 2001 From: "Dennis.Glaeser" Date: Tue, 20 Oct 2020 09:05:55 +0200 Subject: [PATCH 05/32] [test] add test for grid variables construction --- test/discretization/CMakeLists.txt | 6 + test/discretization/test_fvgridvariables.cc | 162 ++++++++++++++++++++ 2 files changed, 168 insertions(+) create mode 100644 test/discretization/test_fvgridvariables.cc diff --git a/test/discretization/CMakeLists.txt b/test/discretization/CMakeLists.txt index 25c1bdd297..c9e5b31c6e 100644 --- a/test/discretization/CMakeLists.txt +++ b/test/discretization/CMakeLists.txt @@ -3,3 +3,9 @@ add_subdirectory(staggered) add_subdirectory(box) add_subdirectory(projection) add_subdirectory(rotationsymmetry) + +dune_add_test(NAME test_fvgridvariables + LABELS unit + SOURCES test_fvgridvariables.cc + COMMAND ./test_fvgridvariables + CMD_ARGS -Problem.Name gridvarstest) diff --git a/test/discretization/test_fvgridvariables.cc b/test/discretization/test_fvgridvariables.cc new file mode 100644 index 0000000000..39ed7a4eba --- /dev/null +++ b/test/discretization/test_fvgridvariables.cc @@ -0,0 +1,162 @@ +// -*- 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 + * \brief Test for the new experimental finite volume grid variables. + */ +#include +#include + +#include +#include +#include +#include + +// we use the 1p type tag here in order not to be obliged +// to define grid flux vars cache & vol vars cache... +#include +#include +#include + +#include +#include + +#include +#include +#include +#include + +namespace Dumux { + +template +class MockSpatialParams +: public FVSpatialParamsOneP> +{ +public: + using PermeabilityType = Scalar; +}; + +template +class MockProblem : public FVProblem +{ + using Parent = FVProblem; +public: + using Parent::Parent; +}; + +namespace Properties { + +// new type tags +namespace TTag { +struct GridVariablesTest { using InheritsFrom = std::tuple; }; +struct GridVariablesTestBox { using InheritsFrom = std::tuple; }; +} // end namespace TTag + +// property definitions +template +struct Grid +{ using type = Dune::YaspGrid<2>; }; + +template +struct Problem +{ using type = MockProblem; }; + +template +struct SpatialParams +{ +private: + using Scalar = GetPropType; + using GG = GetPropType; +public: + using type = MockSpatialParams; +}; + +template +struct FluidSystem +{ +private: + using Scalar = GetPropType; +public: + using type = FluidSystems::OnePLiquid>; +}; + +// use the new experimental grid variables +template +struct GridVariables +{ +private: + using GG = GetPropType; + using GVV = GetPropType; + using GFC = GetPropType; + using X = GetPropType; + +public: + using type = Dumux::Experimental::FVGridVariables; +}; + +} // end namespace Properties +} // end namespace Dumux + +int main (int argc, char *argv[]) +{ + Dune::MPIHelper::instance(argc, argv); + + using namespace Dumux; + Dumux::Parameters::init(argc, argv); + + using TypeTag = Properties::TTag::GridVariablesTestBox; + using Grid = GetPropType; + using GridFactory = Dune::StructuredGridFactory; + auto grid = GridFactory::createCubeGrid({0.0, 0.0}, {1.0, 1.0}, {2, 2}); + + using GridGeometry = GetPropType; + auto gridGeometry = std::make_shared(grid->leafGridView()); + + using Problem = GetPropType; + auto problem = std::make_shared(gridGeometry); + + using GridVariables = GetPropType; + + // This constructor should fail as the problem does not implement initial() + bool caught = false; + try { auto gridVariables = std::make_shared(problem, gridGeometry); } + catch (...) { caught = true; } + if (!caught) + DUNE_THROW(Dune::Exception, "Expected construction to fail"); + + // Construction with a solution + using SolutionVector = GetPropType; + SolutionVector x; x.resize(gridGeometry->numDofs()); x = 0.0; + auto gridVariables = std::make_shared(problem, gridGeometry, x); + + // Construction from a moved solution (TODO: how to check if move succeeded?) + gridVariables = std::make_shared(problem, gridGeometry, std::move(x)); + + // Construction from initializer lambda + const auto init = [gridGeometry] (auto& x) { x.resize(gridGeometry->numDofs()); x = 2.25; }; + gridVariables = std::make_shared(problem, gridGeometry, init); + + const auto& dofs = gridVariables->dofs(); + if (std::any_of(dofs.begin(), dofs.end(), + [] (auto d) { return Dune::FloatCmp::ne(2.25, d[0]); })) + DUNE_THROW(Dune::Exception, "Unexpected dof value"); + + std::cout << "\nAll tests passed" << std::endl; + return 0; +} -- GitLab From dba2f08e31dc75429e2fd9367cf667cc4fbbacae Mon Sep 17 00:00:00 2001 From: Timo Koch Date: Tue, 20 Oct 2020 10:18:59 +0200 Subject: [PATCH 06/32] [newton] introduce variables backend The newton solver linearizes the system of equations around a given state, currently represented by the primary variables. In order to support the state being represented by grid variables, this backend can be used. --- dumux/nonlinear/newtonvariablesbackend.hh | 86 +++++++++++++++++++++++ 1 file changed, 86 insertions(+) create mode 100644 dumux/nonlinear/newtonvariablesbackend.hh diff --git a/dumux/nonlinear/newtonvariablesbackend.hh b/dumux/nonlinear/newtonvariablesbackend.hh new file mode 100644 index 0000000000..a379eb8ec1 --- /dev/null +++ b/dumux/nonlinear/newtonvariablesbackend.hh @@ -0,0 +1,86 @@ +// -*- 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 Nonlinear + * \brief Backends for the Newton for different variable types + */ +#ifndef DUMUX_NEWTON_VARIABLES_BACKEND_HH +#define DUMUX_NEWTON_VARIABLES_BACKEND_HH + +#include +#include + +namespace Dumux { + +/*! + * \file + * \ingroup Nonlinear + * \brief Class providing operations with variables needed in the Newton solver + */ +template +class NewtonVariablesBackend; + +/*! + * \file + * \ingroup Nonlinear + * \brief Class providing Newton operations for scalar/number types + */ +template +class NewtonVariablesBackend::value, Scalar>> +{ +public: + using Variables = Scalar; //!< the type of the variables object + using DofVector = Scalar; //!< the type of the dofs parametrizing the variables object + + static std::size_t size(const DofVector& d) + { return 1; } + + static DofVector makeDofVector(const DofVector& d) + { return d; } + + static DofVector makeZeroDofVector(std::size_t size) + { return 0.0; } + + static Scalar makeDofVectorForSolver(const DofVector& d) + { return d; } + + static Scalar makeZeroDofVectorForSolver(std::size_t size) + { return 0.0; } + + static Scalar reconstructDofVectorFromSolver(const Scalar& d) + { return d; } + + static Scalar maxRelativeShift(const DofVector& previous, const DofVector& delta) + { + const auto current = previous - delta; + using std::abs; using std::max; + auto error = abs(previous - current); + error /= max(1.0, abs(previous + current)/2.0); + return error; + } + + // //! operations on variables + // static DofVector& getDofVector(Variables& v) + // { return v; } +}; + +} // end namespace Dumux + +#endif -- GitLab From 8788322fd3153c24a94c99811d68097b1cddfbe5 Mon Sep 17 00:00:00 2001 From: "Dennis.Glaeser" Date: Wed, 21 Oct 2020 10:40:12 +0200 Subject: [PATCH 07/32] [newton][varsbackend] add backends for vectors & generic variables This commit also removes the interfaces related to the construction of vectors for the solver, as this is currently done in Detail::assign in newtonsolver.hh. The same holds for maxRelativeShift(), which is also currently implemented in the Detail namespace in newtonsolver.hh- In the future, we could think about outsourcing putting these two functions into the backend. --- dumux/nonlinear/newtonvariablesbackend.hh | 157 +++++++++++++++++++--- 1 file changed, 136 insertions(+), 21 deletions(-) diff --git a/dumux/nonlinear/newtonvariablesbackend.hh b/dumux/nonlinear/newtonvariablesbackend.hh index a379eb8ec1..9cb617140b 100644 --- a/dumux/nonlinear/newtonvariablesbackend.hh +++ b/dumux/nonlinear/newtonvariablesbackend.hh @@ -24,29 +24,41 @@ #ifndef DUMUX_NEWTON_VARIABLES_BACKEND_HH #define DUMUX_NEWTON_VARIABLES_BACKEND_HH +#include +#include #include + +#include #include +#include +#include +#include + +// forward declaration +namespace Dune { + +template +class MultiTypeBlockVector; + +} // end namespace Dune namespace Dumux { /*! - * \file * \ingroup Nonlinear - * \brief Class providing operations with variables needed in the Newton solver + * \brief Class providing operations with primary variable vectors */ -template -class NewtonVariablesBackend; +template +class NewtonDofBackend; /*! - * \file * \ingroup Nonlinear - * \brief Class providing Newton operations for scalar/number types + * \brief Specialization providing Newton operations for scalar/number types */ template -class NewtonVariablesBackend::value, Scalar>> +class NewtonDofBackend::value, Scalar>> { public: - using Variables = Scalar; //!< the type of the variables object using DofVector = Scalar; //!< the type of the dofs parametrizing the variables object static std::size_t size(const DofVector& d) @@ -57,30 +69,133 @@ public: static DofVector makeZeroDofVector(std::size_t size) { return 0.0; } +}; + +/*! + * \ingroup Nonlinear + * \brief Specialization providing Newton operations for block vectors + */ +template +class NewtonDofBackend> +{ + +public: + using DofVector = Dune::BlockVector; //!< the type of the dofs parametrizing the variables object - static Scalar makeDofVectorForSolver(const DofVector& d) + static std::size_t size(const DofVector& d) + { return d.size(); } + + static DofVector makeDofVector(const DofVector& d) { return d; } - static Scalar makeZeroDofVectorForSolver(std::size_t size) - { return 0.0; } + static DofVector makeZeroDofVector(std::size_t size) + { DofVector d; d.resize(size); return d; } +}; + +/*! + * \ingroup Nonlinear + * \brief Specialization providing Newton operations for multitype block vectors + */ +template +class NewtonDofBackend> +{ + using DV = Dune::MultiTypeBlockVector; + static constexpr auto numBlocks = DV::size(); + + using VectorSizeInfo = std::array; + +public: + using DofVector = DV; //!< the type of the dofs parametrizing the variables object + + static VectorSizeInfo size(const DofVector& d) + { + VectorSizeInfo result; + using namespace Dune::Hybrid; + forEach(std::make_index_sequence{}, [&](auto i) { + result[i] = d[Dune::index_constant{}].size(); + }); + return result; + } - static Scalar reconstructDofVectorFromSolver(const Scalar& d) + static DofVector makeDofVector(const DofVector& d) { return d; } - static Scalar maxRelativeShift(const DofVector& previous, const DofVector& delta) + static DofVector makeZeroDofVector(const VectorSizeInfo& size) { - const auto current = previous - delta; - using std::abs; using std::max; - auto error = abs(previous - current); - error /= max(1.0, abs(previous + current)/2.0); - return error; + DofVector result; + using namespace Dune::Hybrid; + forEach(std::make_index_sequence{}, [&](auto i) { + result[Dune::index_constant{}].resize(size[i]); + }); + return result; } +}; + +namespace Impl { + +template +using SolutionVectorType = typename Vars::SolutionVector; + +template +class NewtonVariablesBackend; + +/*! + * \ingroup Nonlinear + * \brief Class providing Newton operations for scalar/number types + * \note We assume the variables being simply a dof vector if we + * do not find the variables class to export `SolutionVector`. + */ +template +class NewtonVariablesBackend +: public NewtonDofBackend +{ + using ParentType = NewtonDofBackend; + +public: + using Variables = Vars; + using typename ParentType::DofVector; - // //! operations on variables - // static DofVector& getDofVector(Variables& v) - // { return v; } + static void update(Variables& v, const DofVector& dofs) + { v = dofs; } + + //! operations on variables + static const DofVector& getDofVector(const Variables& v) + { return v; } }; +/*! + * \file + * \ingroup Nonlinear + * \brief Class providing Newton operations for generic variable + * types, possibly containing also secondary variables. + */ +template +class NewtonVariablesBackend +: public NewtonDofBackend +{ +public: + using DofVector = typename Vars::SolutionVector; + using Variables = Vars; //!< the type of the variables object + + static void update(Variables& v, const DofVector& dofs) + { v.update(dofs); } + + //! operations on variables + static const DofVector& getDofVector(const Variables& v) + { return v.dofs(); } +}; +} // end namespace Impl + +/*! + * \ingroup Nonlinear + * \brief Class providing Newton operations for generic variable classes + * that represent the state of a numerical solution, possibly + * consisting of primary/secondary variables and information on + * the time level. + */ +template +using NewtonVariablesBackend = Impl::NewtonVariablesBackend>; + } // end namespace Dumux #endif -- GitLab From 394d5d8f3326053e9323d8404f1cb437ef06f648 Mon Sep 17 00:00:00 2001 From: "Dennis.Glaeser" Date: Wed, 21 Oct 2020 12:02:10 +0200 Subject: [PATCH 08/32] [newtonbackend] provide access to non-const dofvector --- dumux/nonlinear/newtonvariablesbackend.hh | 14 ++++++++++++-- 1 file changed, 12 insertions(+), 2 deletions(-) diff --git a/dumux/nonlinear/newtonvariablesbackend.hh b/dumux/nonlinear/newtonvariablesbackend.hh index 9cb617140b..48d6faab06 100644 --- a/dumux/nonlinear/newtonvariablesbackend.hh +++ b/dumux/nonlinear/newtonvariablesbackend.hh @@ -155,12 +155,17 @@ public: using Variables = Vars; using typename ParentType::DofVector; + //! update to new solution vector static void update(Variables& v, const DofVector& dofs) { v = dofs; } - //! operations on variables + //! return const reference to dof vector static const DofVector& getDofVector(const Variables& v) { return v; } + + //! return reference to dof vector + static DofVector& getDofVector(Variables& v) + { return v; } }; /*! @@ -177,12 +182,17 @@ public: using DofVector = typename Vars::SolutionVector; using Variables = Vars; //!< the type of the variables object + //! update to new solution vector static void update(Variables& v, const DofVector& dofs) { v.update(dofs); } - //! operations on variables + //! return const reference to dof vector static const DofVector& getDofVector(const Variables& v) { return v.dofs(); } + + //! return reference to dof vector + static DofVector& getDofVector(Variables& v) + { return v.dofs(); } }; } // end namespace Impl -- GitLab From 336fca13874a5e029fc3c0c6a6348069ab9bfe09 Mon Sep 17 00:00:00 2001 From: "Dennis.Glaeser" Date: Wed, 21 Oct 2020 11:02:30 +0200 Subject: [PATCH 09/32] [pdesolver] detect and export assembler variables type --- dumux/common/pdesolver.hh | 42 ++++++++++++++++++++++++++++++++------- 1 file changed, 35 insertions(+), 7 deletions(-) diff --git a/dumux/common/pdesolver.hh b/dumux/common/pdesolver.hh index 99f61bcaa6..d11285c7ab 100644 --- a/dumux/common/pdesolver.hh +++ b/dumux/common/pdesolver.hh @@ -28,6 +28,7 @@ #include #include +#include #include @@ -35,9 +36,25 @@ namespace Dune { template class MultiTypeBlockMatrix; -} +} // end namespace Dune namespace Dumux { +namespace Impl { + + template + using AssemblerVariablesType = typename Assembler::Variables; + + template + inline constexpr bool exportsVariables = Dune::Std::is_detected_v; + + template> struct VariablesChooser; + template struct VariablesChooser { using Type = AssemblerVariablesType; }; + template struct VariablesChooser { using Type = typename A::ResidualType; }; + + template + using AssemblerVariables = typename VariablesChooser::Type; + +} // end namespace Impl /*! * \ingroup Common @@ -53,11 +70,19 @@ namespace Dumux { template class PDESolver { - using SolutionVector = typename Assembler::ResidualType; using Scalar = typename Assembler::Scalar; using TimeLoop = TimeLoopBase; public: + + //! export the type of variables that represent a numerical solution + using Variables = Impl::AssemblerVariables; + + /*! + * \brief Constructor + * \param assembler pointer to the assembler of the linear system + * \param linearSolver pointer to the solver of the resulting linear system + */ PDESolver(std::shared_ptr assembler, std::shared_ptr linearSolver) : assembler_(assembler) @@ -68,24 +93,27 @@ public: /*! * \brief Solve the given PDE system (usually assemble + solve linear system + update) - * \param sol a solution vector possbilty containing an initial solution + * \param vars instance of the `Variables` class representing a numerical + * solution, defining primary and possibly secondary variables + * and information on the time level. */ - virtual void solve(SolutionVector& sol) = 0; + virtual void solve(Variables& vars) = 0; /*! * \brief Solve the given PDE system with time step control * \note This is used for solvers that are allowed to e.g. automatically reduce the * time step if the solve was not successful - * \param sol a solution vector possbilty containing an initial solution + * \param vars instance of the `Variables` class representing a numerical solution * \param timeLoop a reference to the current time loop */ - virtual void solve(SolutionVector& sol, TimeLoop& timeLoop) + virtual void solve(Variables& vars, TimeLoop& timeLoop) { // per default we just forward to the method without time step control - solve(sol); + solve(vars); } protected: + /*! * \brief Access the assembler */ -- GitLab From e7d8af35b1511b14ba34214f4e10dbd6a7cd6039 Mon Sep 17 00:00:00 2001 From: "Dennis.Glaeser" Date: Wed, 21 Oct 2020 11:34:57 +0200 Subject: [PATCH 10/32] [newton] use variables backend --- dumux/nonlinear/newtonsolver.hh | 228 ++++++++++++++++++++------------ 1 file changed, 146 insertions(+), 82 deletions(-) diff --git a/dumux/nonlinear/newtonsolver.hh b/dumux/nonlinear/newtonsolver.hh index 3265ef8687..3777cdc778 100644 --- a/dumux/nonlinear/newtonsolver.hh +++ b/dumux/nonlinear/newtonsolver.hh @@ -54,6 +54,7 @@ #include #include "newtonconvergencewriter.hh" +#include "newtonvariablesbackend.hh" namespace Dumux { namespace Detail { @@ -176,18 +177,21 @@ template { using ParentType = PDESolver; + using Backend = NewtonVariablesBackend; + using SolutionVector = typename Backend::DofVector; + using Scalar = typename Assembler::Scalar; using JacobianMatrix = typename Assembler::JacobianMatrix; - using SolutionVector = typename Assembler::ResidualType; using ConvergenceWriter = ConvergenceWriterInterface; using TimeLoop = TimeLoopBase; using PrimaryVariableSwitch = typename Detail::GetPVSwitch::type; using HasPriVarsSwitch = typename Detail::GetPVSwitch::value_t; // std::true_type or std::false_type static constexpr bool hasPriVarsSwitch() { return HasPriVarsSwitch::value; }; + static constexpr bool assemblerExportsVariables = Impl::exportsVariables; public: - + using typename ParentType::Variables; using Communication = Comm; /*! @@ -288,8 +292,10 @@ public: /*! * \brief Run the Newton method to solve a non-linear system. * Does time step control when the Newton fails to converge + * \param vars The variables object representing the current state of the + * numerical solution (primary and possibly secondary variables). */ - void solve(SolutionVector& uCurrentIter, TimeLoop& timeLoop) override + void solve(Variables& vars, TimeLoop& timeLoop) override { if (this->assembler().isStationaryProblem()) DUNE_THROW(Dune::InvalidStateException, "Using time step control with stationary problem makes no sense!"); @@ -298,7 +304,7 @@ public: for (std::size_t i = 0; i <= maxTimeStepDivisions_; ++i) { // linearize & solve - const bool converged = solve_(uCurrentIter); + const bool converged = solve_(vars); if (converged) return; @@ -306,10 +312,17 @@ public: else if (!converged && i < maxTimeStepDivisions_) { // set solution to previous solution - uCurrentIter = this->assembler().prevSol(); - - // reset the grid variables to the previous solution - this->assembler().resetTimeStep(uCurrentIter); + Backend::update(vars, this->assembler().prevSol()); + + // if this is true, we assume old-style assemblers/grid variables + // TODO: reset time step is more efficient as it simply copies + // prevVolVars into curVolVars. Above (in the new style) + // there is an update. We should think about how to achieve + // this more efficiently. + // TODO: Is there a test with time step reductions? It should + // probably be tested if this works properly. + if (!assemblerExportsVariables) + this->assembler().resetTimeStep(Backend::getDofVector(vars)); if (verbosity_ >= 1) { @@ -334,10 +347,12 @@ public: /*! * \brief Run the Newton method to solve a non-linear system. * The solver is responsible for all the strategic decisions. + * \param vars The variables object representing the current state of the + * numerical solution (primary and possibly secondary variables). */ - void solve(SolutionVector& uCurrentIter) override + void solve(Variables& vars) override { - const bool converged = solve_(uCurrentIter); + const bool converged = solve_(vars); if (!converged) DUNE_THROW(NumericalProblem, Fmt::format("Newton solver didn't converge after {} iterations.\n", numSteps_)); @@ -347,36 +362,36 @@ public: * \brief Called before the Newton method is applied to an * non-linear system of equations. * - * \param u The initial solution + * \param initVars The variables representing the initial solution */ - virtual void newtonBegin(SolutionVector& u) + virtual void newtonBegin(Variables& initVars) { numSteps_ = 0; - initPriVarSwitch_(u, HasPriVarsSwitch{}); + initPriVarSwitch_(initVars, HasPriVarsSwitch{}); + const auto& initSol = Backend::getDofVector(initVars); // write the initial residual if a convergence writer was set if (convergenceWriter_) { - this->assembler().assembleResidual(u); - SolutionVector delta(u); - delta = 0; // dummy vector, there is no delta before solving the linear system - convergenceWriter_->write(u, delta, this->assembler().residual()); + this->assembler().assembleResidual(initVars); + auto delta = Backend::makeZeroDofVector(Backend::size(initSol)); + convergenceWriter_->write(initSol, delta, this->assembler().residual()); } if (enablePartialReassembly_) { partialReassembler_->resetColors(); - resizeDistanceFromLastLinearization_(u, distanceFromLastLinearization_); + resizeDistanceFromLastLinearization_(initSol, distanceFromLastLinearization_); } } /*! * \brief Returns true if another iteration should be done. * - * \param uCurrentIter The solution of the current Newton iteration + * \param varsCurrentIter The variables of the current Newton iteration * \param converged if the Newton method's convergence criterion was met in this step */ - virtual bool newtonProceed(const SolutionVector &uCurrentIter, bool converged) + virtual bool newtonProceed(const Variables &varsCurrentIter, bool converged) { if (numSteps_ < minSteps_) return true; @@ -399,7 +414,7 @@ public: /*! * \brief Indicates the beginning of a Newton iteration. */ - virtual void newtonBeginStep(const SolutionVector& u) + virtual void newtonBeginStep(const Variables& vars) { lastShift_ = shift_; if (numSteps_ == 0) @@ -415,11 +430,11 @@ public: /*! * \brief Assemble the linear system of equations \f$\mathbf{A}x - b = 0\f$. * - * \param uCurrentIter The current iteration's solution vector + * \param vars The current iteration's variables */ - virtual void assembleLinearSystem(const SolutionVector& uCurrentIter) + virtual void assembleLinearSystem(const Variables& vars) { - assembleLinearSystem_(this->assembler(), uCurrentIter); + assembleLinearSystem_(this->assembler(), vars); if (enablePartialReassembly_) partialReassembler_->report(comm_, endIterMsgStream_); @@ -499,15 +514,15 @@ public: * subtract deltaU from uLastIter, i.e. * \f[ u^{k+1} = u^k - \Delta u^k \f] * - * \param uCurrentIter The solution vector after the current iteration + * \param vars The variables after the current iteration * \param uLastIter The solution vector after the last iteration * \param deltaU The delta as calculated from solving the linear * system of equations. This parameter also stores * the updated solution. */ - void newtonUpdate(SolutionVector &uCurrentIter, - const SolutionVector &uLastIter, - const SolutionVector &deltaU) + void newtonUpdate(Variables& vars, + const SolutionVector& uLastIter, + const SolutionVector& deltaU) { if (enableShiftCriterion_ || enablePartialReassembly_) newtonUpdateShift_(uLastIter, deltaU); @@ -548,32 +563,32 @@ public: } if (useLineSearch_) - lineSearchUpdate_(uCurrentIter, uLastIter, deltaU); + lineSearchUpdate_(vars, uLastIter, deltaU); else if (useChop_) - choppedUpdate_(uCurrentIter, uLastIter, deltaU); + choppedUpdate_(vars, uLastIter, deltaU); else { - uCurrentIter = uLastIter; + auto uCurrentIter = Backend::makeDofVector(uLastIter); uCurrentIter -= deltaU; - solutionChanged_(uCurrentIter); + solutionChanged_(vars, uCurrentIter); if (enableResidualCriterion_) - computeResidualReduction_(uCurrentIter); + computeResidualReduction_(vars); } } /*! * \brief Indicates that one Newton iteration was finished. * - * \param uCurrentIter The solution after the current Newton iteration + * \param vars The variables after the current Newton iteration * \param uLastIter The solution at the beginning of the current Newton iteration */ - virtual void newtonEndStep(SolutionVector &uCurrentIter, + virtual void newtonEndStep(Variables &vars, const SolutionVector &uLastIter) { - invokePriVarSwitch_(uCurrentIter, HasPriVarsSwitch{}); + invokePriVarSwitch_(vars, HasPriVarsSwitch{}); ++numSteps_; @@ -785,22 +800,37 @@ public: protected: /*! - * \brief Update solution-depended quantities like grid variables after the solution has changed. + * \brief Update solution-dependent quantities like grid variables after the solution has changed. + * \todo TODO: In case we stop support for old-style grid variables / assemblers at one point, + * this would become obsolete as only the update call to the backend would remain. */ - virtual void solutionChanged_(const SolutionVector &uCurrentIter) + virtual void solutionChanged_(Variables& vars, const SolutionVector& uCurrentIter) { - this->assembler().updateGridVariables(uCurrentIter); + Backend::update(vars, uCurrentIter); + + if constexpr (!assemblerExportsVariables) + this->assembler().updateGridVariables(Backend::getDofVector(vars)); } - void computeResidualReduction_(const SolutionVector &uCurrentIter) + void computeResidualReduction_(const Variables& vars) { + // we assume that the assembler works on solution vectors + // if it doesn't export the variables type if constexpr (Detail::hasNorm()) { - this->assembler().assembleResidual(uCurrentIter); + if constexpr (!assemblerExportsVariables) + this->assembler().assembleResidual(Backend::getDofVector(vars)); + else + this->assembler().assembleResidual(vars); residualNorm_ = this->linearSolver().norm(this->assembler().residual()); } else - residualNorm_ = this->assembler().residualNorm(uCurrentIter); + { + if constexpr (!assemblerExportsVariables) + residualNorm_ = this->assembler().residualNorm(Backend::getDofVector(vars)); + else + residualNorm_ = this->assembler().residualNorm(vars); + } reduction_ = residualNorm_; reduction_ /= initialResidual_; @@ -812,51 +842,80 @@ protected: /*! * \brief Initialize the privar switch, noop if there is no priVarSwitch */ - void initPriVarSwitch_(SolutionVector&, std::false_type) {} + void initPriVarSwitch_(Variables&, std::false_type) {} /*! * \brief Initialize the privar switch */ - void initPriVarSwitch_(SolutionVector& sol, std::true_type) + void initPriVarSwitch_(Variables& vars, std::true_type) { - priVarSwitch_->reset(sol.size()); + priVarSwitch_->reset(Backend::size(Backend::getDofVector(vars))); priVarsSwitchedInLastIteration_ = false; const auto& problem = this->assembler().problem(); const auto& gridGeometry = this->assembler().gridGeometry(); - auto& gridVariables = this->assembler().gridVariables(); - priVarSwitch_->updateDirichletConstraints(problem, gridGeometry, gridVariables, sol); + + // if the assembler does not export the variables type + // we assume old-style assemblers which operate on solution + // vectors and return the grid variables + if constexpr (!assemblerExportsVariables) + priVarSwitch_->updateDirichletConstraints(problem, gridGeometry, + this->assembler().gridVariables(), + vars); + else + priVarSwitch_->updateDirichletConstraints(problem, gridGeometry, vars, Backend::getDofVector(vars)); } /*! * \brief Switch primary variables if necessary, noop if there is no priVarSwitch */ - void invokePriVarSwitch_(SolutionVector&, std::false_type) {} + void invokePriVarSwitch_(Variables&, std::false_type) {} /*! * \brief Switch primary variables if necessary */ - void invokePriVarSwitch_(SolutionVector& uCurrentIter, std::true_type) + void invokePriVarSwitch_(Variables& vars, std::true_type) { // update the variable switch (returns true if the pri vars at at least one dof were switched) // for disabled grid variable caching const auto& gridGeometry = this->assembler().gridGeometry(); const auto& problem = this->assembler().problem(); - auto& gridVariables = this->assembler().gridVariables(); - // invoke the primary variable switch - priVarsSwitchedInLastIteration_ = priVarSwitch_->update(uCurrentIter, gridVariables, - problem, gridGeometry); + // we assume an old-style assembler if it doesn't export the Variable type + if constexpr (!assemblerExportsVariables) + { + auto& gridVariables = this->assembler().gridVariables(); - if (priVarsSwitchedInLastIteration_) + // invoke the primary variable switch + priVarsSwitchedInLastIteration_ = priVarSwitch_->update(vars, gridVariables, problem, gridGeometry); + + if (priVarsSwitchedInLastIteration_) + { + for (const auto& element : elements(gridGeometry.gridView())) + { + // if the volume variables are cached globally, we need to update those where the primary variables have been switched + priVarSwitch_->updateSwitchedVolVars(problem, element, gridGeometry, gridVariables, Backend::getDofVector(vars)); + + // if the flux variables are cached globally, we need to update those where the primary variables have been switched + priVarSwitch_->updateSwitchedFluxVarsCache(problem, element, gridGeometry, gridVariables, Backend::getDofVector(vars)); + } + } + } + else { - for (const auto& element : elements(gridGeometry.gridView())) + // invoke the primary variable switch + priVarsSwitchedInLastIteration_ = priVarSwitch_->update(Backend::getDofVector(vars), vars, problem, gridGeometry); + + if (priVarsSwitchedInLastIteration_) { - // if the volume variables are cached globally, we need to update those where the primary variables have been switched - priVarSwitch_->updateSwitchedVolVars(problem, element, gridGeometry, gridVariables, uCurrentIter); + for (const auto& element : elements(gridGeometry.gridView())) + { + // if the volume variables are cached globally, we need to update those where the primary variables have been switched + priVarSwitch_->updateSwitchedVolVars(problem, element, gridGeometry, vars, Backend::getDofVector(vars)); - // if the flux variables are cached globally, we need to update those where the primary variables have been switched - priVarSwitch_->updateSwitchedFluxVarsCache(problem, element, gridGeometry, gridVariables, uCurrentIter); + // if the flux variables are cached globally, we need to update those where the primary variables have been switched + priVarSwitch_->updateSwitchedFluxVarsCache(problem, element, gridGeometry, vars, Backend::getDofVector(vars)); + } } } } @@ -890,16 +949,16 @@ private: * \brief Run the Newton method to solve a non-linear system. * The solver is responsible for all the strategic decisions. */ - bool solve_(SolutionVector& uCurrentIter) + bool solve_(Variables& vars) { try { // newtonBegin may manipulate the solution - newtonBegin(uCurrentIter); + newtonBegin(vars); // the given solution is the initial guess - SolutionVector uLastIter(uCurrentIter); - SolutionVector deltaU(uCurrentIter); + auto uLastIter = Backend::makeDofVector(Backend::getDofVector(vars)); + auto deltaU = Backend::makeDofVector(Backend::getDofVector(vars)); // setup timers Dune::Timer assembleTimer(false); @@ -909,15 +968,15 @@ private: // execute the method as long as the solver thinks // that we should do another iteration bool converged = false; - while (newtonProceed(uCurrentIter, converged)) + while (newtonProceed(vars, converged)) { // notify the solver that we're about to start // a new iteration - newtonBeginStep(uCurrentIter); + newtonBeginStep(vars); // make the current solution to the old one if (numSteps_ > 0) - uLastIter = uCurrentIter; + uLastIter = Backend::getDofVector(vars); if (verbosity_ >= 1 && enableDynamicOutput_) std::cout << "Assemble: r(x^k) = dS/dt + div F - q; M = grad r" @@ -929,7 +988,7 @@ private: // linearize the problem at the current solution assembleTimer.start(); - assembleLinearSystem(uCurrentIter); + assembleLinearSystem(vars); assembleTimer.stop(); /////////////// @@ -964,17 +1023,17 @@ private: updateTimer.start(); // update the current solution (i.e. uOld) with the delta // (i.e. u). The result is stored in u - newtonUpdate(uCurrentIter, uLastIter, deltaU); + newtonUpdate(vars, uLastIter, deltaU); updateTimer.stop(); // tell the solver that we're done with this iteration - newtonEndStep(uCurrentIter, uLastIter); + newtonEndStep(vars, uLastIter); // if a convergence writer was specified compute residual and write output if (convergenceWriter_) { - this->assembler().assembleResidual(uCurrentIter); - convergenceWriter_->write(uCurrentIter, deltaU, this->assembler().residual()); + this->assembler().assembleResidual(vars); + convergenceWriter_->write(Backend::getDofVector(vars), deltaU, this->assembler().residual()); } // detect if the method has converged @@ -988,6 +1047,8 @@ private: if (!newtonConverged()) { totalWastedIter_ += numSteps_; + // TODO: what should NewtonFail receive as arg? + auto uCurrentIter = Backend::getDofVector(vars); newtonFail(uCurrentIter); return false; } @@ -1014,6 +1075,9 @@ private: std::cout << "Newton: Caught exception: \"" << e.what() << "\"\n"; totalWastedIter_ += numSteps_; + + // TODO: what should NewtonFail receive as arg? + auto uCurrentIter = Backend::getDofVector(vars); newtonFail(uCurrentIter); return false; } @@ -1021,18 +1085,18 @@ private: //! assembleLinearSystem_ for assemblers that support partial reassembly template - auto assembleLinearSystem_(const A& assembler, const SolutionVector& uCurrentIter) + auto assembleLinearSystem_(const A& assembler, const Variables& vars) -> typename std::enable_if_t { - this->assembler().assembleJacobianAndResidual(uCurrentIter, partialReassembler_.get()); + this->assembler().assembleJacobianAndResidual(vars, partialReassembler_.get()); } //! assembleLinearSystem_ for assemblers that don't support partial reassembly template - auto assembleLinearSystem_(const A& assembler, const SolutionVector& uCurrentIter) + auto assembleLinearSystem_(const A& assembler, const Variables& vars) -> typename std::enable_if_t { - this->assembler().assembleJacobianAndResidual(uCurrentIter); + this->assembler().assembleJacobianAndResidual(vars); } /*! @@ -1053,7 +1117,7 @@ private: shift_ = comm_.max(shift_); } - virtual void lineSearchUpdate_(SolutionVector &uCurrentIter, + virtual void lineSearchUpdate_(Variables &vars, const SolutionVector &uLastIter, const SolutionVector &deltaU) { @@ -1061,12 +1125,12 @@ private: while (true) { - uCurrentIter = deltaU; + auto uCurrentIter = Backend::makeDofVector(deltaU); uCurrentIter *= -lambda; uCurrentIter += uLastIter; - solutionChanged_(uCurrentIter); + solutionChanged_(vars, uCurrentIter); - computeResidualReduction_(uCurrentIter); + computeResidualReduction_(vars); if (reduction_ < lastReduction_ || lambda <= 0.125) { endIterMsgStream_ << Fmt::format(", residual reduction {:.4e}->{:.4e}@lambda={:.4f}", lastReduction_, reduction_, lambda); @@ -1079,9 +1143,9 @@ private: } //! \note method must update the gridVariables, too! - virtual void choppedUpdate_(SolutionVector &uCurrentIter, - const SolutionVector &uLastIter, - const SolutionVector &deltaU) + virtual void choppedUpdate_(Variables& vars, + const SolutionVector& uLastIter, + const SolutionVector& deltaU) { DUNE_THROW(Dune::NotImplemented, "Chopped Newton update strategy not implemented."); -- GitLab From 6b8956923c67ddbc8165538fd04f52815ca80068 Mon Sep 17 00:00:00 2001 From: "Dennis.Glaeser" Date: Wed, 21 Oct 2020 13:29:33 +0200 Subject: [PATCH 11/32] [newton] allow access to Backend type in derived classes --- dumux/nonlinear/newtonsolver.hh | 3 +++ 1 file changed, 3 insertions(+) diff --git a/dumux/nonlinear/newtonsolver.hh b/dumux/nonlinear/newtonsolver.hh index 3777cdc778..04077aad08 100644 --- a/dumux/nonlinear/newtonsolver.hh +++ b/dumux/nonlinear/newtonsolver.hh @@ -177,9 +177,12 @@ template { using ParentType = PDESolver; + +protected: using Backend = NewtonVariablesBackend; using SolutionVector = typename Backend::DofVector; +private: using Scalar = typename Assembler::Scalar; using JacobianMatrix = typename Assembler::JacobianMatrix; using ConvergenceWriter = ConvergenceWriterInterface; -- GitLab From a86baa9d80fe8dbbad27f8b1a908430ffb88d33e Mon Sep 17 00:00:00 2001 From: "Dennis.Glaeser" Date: Wed, 21 Oct 2020 13:29:58 +0200 Subject: [PATCH 12/32] [richards][newton] use new backend --- dumux/porousmediumflow/richards/newtonsolver.hh | 16 ++++++++++------ 1 file changed, 10 insertions(+), 6 deletions(-) diff --git a/dumux/porousmediumflow/richards/newtonsolver.hh b/dumux/porousmediumflow/richards/newtonsolver.hh index 007506df95..455d78e6b7 100644 --- a/dumux/porousmediumflow/richards/newtonsolver.hh +++ b/dumux/porousmediumflow/richards/newtonsolver.hh @@ -45,27 +45,31 @@ class RichardsNewtonSolver : public NewtonSolver { using Scalar = typename Assembler::Scalar; using ParentType = NewtonSolver; - using SolutionVector = typename Assembler::ResidualType; using Indices = typename Assembler::GridVariables::VolumeVariables::Indices; enum { pressureIdx = Indices::pressureIdx }; + using typename ParentType::Backend; + using typename ParentType::SolutionVector; + public: using ParentType::ParentType; + using typename ParentType::Variables; private: /*! * \brief Update the current solution of the Newton method * - * \param uCurrentIter The solution after the current Newton iteration \f$ u^{k+1} \f$ + * \param varsCurrentIter The variables after the current Newton iteration \f$ u^{k+1} \f$ * \param uLastIter The solution after the last Newton iteration \f$ u^k \f$ * \param deltaU The vector of differences between the last * iterative solution and the next one \f$ \Delta u^k \f$ */ - void choppedUpdate_(SolutionVector &uCurrentIter, + void choppedUpdate_(Variables &varsCurrentIter, const SolutionVector &uLastIter, const SolutionVector &deltaU) final { + auto uCurrentIter = Backend::getDofVector(varsCurrentIter); uCurrentIter = uLastIter; uCurrentIter -= deltaU; @@ -110,11 +114,11 @@ private: } } - // update the grid variables - this->solutionChanged_(uCurrentIter); + // update the variables + this->solutionChanged_(varsCurrentIter, uCurrentIter); if (this->enableResidualCriterion()) - this->computeResidualReduction_(uCurrentIter); + this->computeResidualReduction_(varsCurrentIter); } }; -- GitLab From 85b6c76ad8843a5f2706cfb838ef9c273207a1e9 Mon Sep 17 00:00:00 2001 From: "Dennis.Glaeser" Date: Wed, 21 Oct 2020 13:30:14 +0200 Subject: [PATCH 13/32] [nonequil][newton] use new backend --- .../nonequilibrium/newtonsolver.hh | 17 ++++++++++++----- 1 file changed, 12 insertions(+), 5 deletions(-) diff --git a/dumux/porousmediumflow/nonequilibrium/newtonsolver.hh b/dumux/porousmediumflow/nonequilibrium/newtonsolver.hh index 87b34a8e4b..e0c5c51e17 100644 --- a/dumux/porousmediumflow/nonequilibrium/newtonsolver.hh +++ b/dumux/porousmediumflow/nonequilibrium/newtonsolver.hh @@ -40,20 +40,27 @@ template class NonEquilibriumNewtonSolver : public NewtonSolver { using ParentType = NewtonSolver; - using SolutionVector = typename Assembler::ResidualType; + + using typename ParentType::Backend; + using typename ParentType::SolutionVector; + static constexpr bool assemblerExportsVariables = Impl::exportsVariables; public: using ParentType::ParentType; + using typename ParentType::Variables; - void newtonEndStep(SolutionVector &uCurrentIter, + void newtonEndStep(Variables &varsCurrentIter, const SolutionVector &uLastIter) final { - ParentType::newtonEndStep(uCurrentIter, uLastIter); + ParentType::newtonEndStep(varsCurrentIter, uLastIter); + const auto& uCurrentIter = Backend::getDofVector(varsCurrentIter); - auto& gridVariables = this->assembler().gridVariables(); // Averages the face velocities of a vertex. Implemented in the model. // The velocities are stored in the model. - gridVariables.calcVelocityAverage(uCurrentIter); + if constexpr(!assemblerExportsVariables) + this->assembler().gridVariables().calcVelocityAverage(uCurrentIter); + else + varsCurrentIter.calcVelocityAverage(uCurrentIter); } }; -- GitLab From 0f202cce2831fa94231527b846f6cf04bb17cb4c Mon Sep 17 00:00:00 2001 From: "Dennis.Glaeser" Date: Wed, 21 Oct 2020 13:30:32 +0200 Subject: [PATCH 14/32] [md][newton] use new backend --- dumux/multidomain/newtonsolver.hh | 62 ++++++++++++++++++------------- 1 file changed, 36 insertions(+), 26 deletions(-) diff --git a/dumux/multidomain/newtonsolver.hh b/dumux/multidomain/newtonsolver.hh index dba8c58a41..898b263e89 100644 --- a/dumux/multidomain/newtonsolver.hh +++ b/dumux/multidomain/newtonsolver.hh @@ -50,7 +50,10 @@ template { using ParentType = NewtonSolver; - using SolutionVector = typename Assembler::ResidualType; + using typename ParentType::Backend; + using typename ParentType::SolutionVector; + + static constexpr bool assemblerExportsVariables = Impl::exportsVariables; template using PrimaryVariableSwitch = typename Detail::GetPVSwitchMultiDomain::type; @@ -63,6 +66,7 @@ class MultiDomainNewtonSolver: public NewtonSolver; public: + using typename ParentType::Variables; /*! * \brief The constructor @@ -89,10 +93,10 @@ public: /*! * \brief Indicates the beginning of a Newton iteration. */ - void newtonBeginStep(const SolutionVector& uCurrentIter) override + void newtonBeginStep(const Variables& varsCurrentIter) override { - ParentType::newtonBeginStep(uCurrentIter); - couplingManager_->updateSolution(uCurrentIter); + ParentType::newtonBeginStep(varsCurrentIter); + couplingManager_->updateSolution(Backend::getDofVector(varsCurrentIter)); } /*! @@ -102,14 +106,14 @@ public: * * \param u The initial solution */ - void newtonBegin(SolutionVector& u) override + void newtonBegin(Variables& vars) override { - ParentType::newtonBegin(u); + ParentType::newtonBegin(vars); using namespace Dune::Hybrid; forEach(std::make_index_sequence{}, [&](auto&& id) { - this->initPriVarSwitch_(u, id, HasPriVarsSwitch::value>{}); + this->initPriVarSwitch_(vars, id, HasPriVarsSwitch::value>{}); }); } @@ -129,19 +133,24 @@ public: /*! * \brief Indicates that one Newton iteration was finished. * - * \param uCurrentIter The solution after the current Newton iteration + * \param varsCurrentIter The variables after the current Newton iteration * \param uLastIter The solution at the beginning of the current Newton iteration */ - void newtonEndStep(SolutionVector &uCurrentIter, const SolutionVector &uLastIter) override + void newtonEndStep(Variables& varsCurrentIter, const SolutionVector& uLastIter) override { using namespace Dune::Hybrid; forEach(std::make_index_sequence{}, [&](auto&& id) { - this->invokePriVarSwitch_(uCurrentIter[id], id, HasPriVarsSwitch::value>{}); + auto& uCurrentIter = Backend::getDofVector(varsCurrentIter)[id]; + if constexpr (!assemblerExportsVariables) + this->invokePriVarSwitch_(this->assembler().gridVariables(id), + uCurrentIter, id, HasPriVarsSwitch::value>{}); + else + this->invokePriVarSwitch_(varsCurrentIter[id], uCurrentIter, id, HasPriVarsSwitch::value>{}); }); - ParentType::newtonEndStep(uCurrentIter, uLastIter); - couplingManager_->updateSolution(uCurrentIter); + ParentType::newtonEndStep(varsCurrentIter, uLastIter); + couplingManager_->updateSolution(Backend::getDofVector(varsCurrentIter)); } private: @@ -150,60 +159,61 @@ private: * \brief Reset the privar switch state, noop if there is no priVarSwitch */ template - void initPriVarSwitch_(SolutionVector&, Dune::index_constant id, std::false_type) {} + void initPriVarSwitch_(Variables&, Dune::index_constant id, std::false_type) {} /*! * \brief Switch primary variables if necessary */ template - void initPriVarSwitch_(SolutionVector& sol, Dune::index_constant id, std::true_type) + void initPriVarSwitch_(Variables& vars, Dune::index_constant id, std::true_type) { using namespace Dune::Hybrid; auto& priVarSwitch = *elementAt(priVarSwitches_, id); + auto& sol = Backend::getDofVector(vars)[id]; - priVarSwitch.reset(sol[id].size()); + priVarSwitch.reset(sol.size()); priVarsSwitchedInLastIteration_[i] = false; const auto& problem = this->assembler().problem(id); const auto& gridGeometry = this->assembler().gridGeometry(id); - auto& gridVariables = this->assembler().gridVariables(id); - priVarSwitch.updateDirichletConstraints(problem, gridGeometry, gridVariables, sol[id]); + if constexpr (!assemblerExportsVariables) + priVarSwitch.updateDirichletConstraints(problem, gridGeometry, this->assembler().gridVariables(id), sol); + else // This expects usage of MultiDomainGridVariables + priVarSwitch.updateDirichletConstraints(problem, gridGeometry, vars[id], sol[id]); } /*! * \brief Switch primary variables if necessary, noop if there is no priVarSwitch */ - template - void invokePriVarSwitch_(SubSol&, Dune::index_constant id, std::false_type) {} + template + void invokePriVarSwitch_(SubVars&, SubSol&, Dune::index_constant id, std::false_type) {} /*! * \brief Switch primary variables if necessary */ - template - void invokePriVarSwitch_(SubSol& uCurrentIter, Dune::index_constant id, std::true_type) + template + void invokePriVarSwitch_(SubVars& subVars, SubSol& uCurrentIter, Dune::index_constant id, std::true_type) { // update the variable switch (returns true if the pri vars at at least one dof were switched) // for disabled grid variable caching const auto& gridGeometry = this->assembler().gridGeometry(id); const auto& problem = this->assembler().problem(id); - auto& gridVariables = this->assembler().gridVariables(id); using namespace Dune::Hybrid; auto& priVarSwitch = *elementAt(priVarSwitches_, id); // invoke the primary variable switch - priVarsSwitchedInLastIteration_[i] = priVarSwitch.update(uCurrentIter, gridVariables, - problem, gridGeometry); + priVarsSwitchedInLastIteration_[i] = priVarSwitch.update(uCurrentIter, subVars, problem, gridGeometry); if (priVarsSwitchedInLastIteration_[i]) { for (const auto& element : elements(gridGeometry.gridView())) { // if the volume variables are cached globally, we need to update those where the primary variables have been switched - priVarSwitch.updateSwitchedVolVars(problem, element, gridGeometry, gridVariables, uCurrentIter); + priVarSwitch.updateSwitchedVolVars(problem, element, gridGeometry, subVars, uCurrentIter); // if the flux variables are cached globally, we need to update those where the primary variables have been switched - priVarSwitch.updateSwitchedFluxVarsCache(problem, element, gridGeometry, gridVariables, uCurrentIter); + priVarSwitch.updateSwitchedFluxVarsCache(problem, element, gridGeometry, subVars, uCurrentIter); } } } -- GitLab From 397e2c39a1018ed923fcc70af76a8e58910d6f7a Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Dennis=20Gl=C3=A4ser?= Date: Fri, 5 Feb 2021 12:40:55 +0100 Subject: [PATCH 15/32] [newton][varsbackend] remove makeDofVector() all implementation rely on the copy constructor, anyway. In case support for some solution vector type that is not copy constructible is required at one point, such logic can be reintroduced. --- dumux/nonlinear/newtonsolver.hh | 8 ++++---- dumux/nonlinear/newtonvariablesbackend.hh | 9 --------- 2 files changed, 4 insertions(+), 13 deletions(-) diff --git a/dumux/nonlinear/newtonsolver.hh b/dumux/nonlinear/newtonsolver.hh index 04077aad08..4db4b094f1 100644 --- a/dumux/nonlinear/newtonsolver.hh +++ b/dumux/nonlinear/newtonsolver.hh @@ -573,7 +573,7 @@ public: else { - auto uCurrentIter = Backend::makeDofVector(uLastIter); + auto uCurrentIter = uLastIter; uCurrentIter -= deltaU; solutionChanged_(vars, uCurrentIter); @@ -960,8 +960,8 @@ private: newtonBegin(vars); // the given solution is the initial guess - auto uLastIter = Backend::makeDofVector(Backend::getDofVector(vars)); - auto deltaU = Backend::makeDofVector(Backend::getDofVector(vars)); + auto uLastIter = Backend::getDofVector(vars); + auto deltaU = Backend::getDofVector(vars); // setup timers Dune::Timer assembleTimer(false); @@ -1128,7 +1128,7 @@ private: while (true) { - auto uCurrentIter = Backend::makeDofVector(deltaU); + auto uCurrentIter = deltaU; uCurrentIter *= -lambda; uCurrentIter += uLastIter; solutionChanged_(vars, uCurrentIter); diff --git a/dumux/nonlinear/newtonvariablesbackend.hh b/dumux/nonlinear/newtonvariablesbackend.hh index 48d6faab06..a154671f8c 100644 --- a/dumux/nonlinear/newtonvariablesbackend.hh +++ b/dumux/nonlinear/newtonvariablesbackend.hh @@ -64,9 +64,6 @@ public: static std::size_t size(const DofVector& d) { return 1; } - static DofVector makeDofVector(const DofVector& d) - { return d; } - static DofVector makeZeroDofVector(std::size_t size) { return 0.0; } }; @@ -85,9 +82,6 @@ public: static std::size_t size(const DofVector& d) { return d.size(); } - static DofVector makeDofVector(const DofVector& d) - { return d; } - static DofVector makeZeroDofVector(std::size_t size) { DofVector d; d.resize(size); return d; } }; @@ -117,9 +111,6 @@ public: return result; } - static DofVector makeDofVector(const DofVector& d) - { return d; } - static DofVector makeZeroDofVector(const VectorSizeInfo& size) { DofVector result; -- GitLab From d07494dc682724cd8e984763e556e477acd81c05 Mon Sep 17 00:00:00 2001 From: Timo Koch Date: Mon, 25 May 2020 17:32:28 +0200 Subject: [PATCH 16/32] [newton] Extract privarswitch implementation to separate class --- dumux/nonlinear/CMakeLists.txt | 1 + dumux/nonlinear/newtonsolver.hh | 125 +++------------- .../nonlinear/primaryvariableswitchadapter.hh | 137 ++++++++++++++++++ 3 files changed, 160 insertions(+), 103 deletions(-) create mode 100644 dumux/nonlinear/primaryvariableswitchadapter.hh diff --git a/dumux/nonlinear/CMakeLists.txt b/dumux/nonlinear/CMakeLists.txt index 2c5c2050a3..0cdeabcdf5 100644 --- a/dumux/nonlinear/CMakeLists.txt +++ b/dumux/nonlinear/CMakeLists.txt @@ -2,5 +2,6 @@ install(FILES findscalarroot.hh newtonconvergencewriter.hh newtonsolver.hh +primaryvariableswitchadapter.hh staggerednewtonconvergencewriter.hh DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/dumux/nonlinear) diff --git a/dumux/nonlinear/newtonsolver.hh b/dumux/nonlinear/newtonsolver.hh index 4db4b094f1..e12661ff3b 100644 --- a/dumux/nonlinear/newtonsolver.hh +++ b/dumux/nonlinear/newtonsolver.hh @@ -55,6 +55,7 @@ #include "newtonconvergencewriter.hh" #include "newtonvariablesbackend.hh" +#include "primaryvariableswitchadapter.hh" namespace Dumux { namespace Detail { @@ -69,13 +70,6 @@ struct supportsPartialReassembly {} }; -//! helper aliases to extract a primary variable switch from the VolumeVariables (if defined, yields int otherwise) -template -using DetectPVSwitch = typename Assembler::GridVariables::VolumeVariables::PrimaryVariableSwitch; - -template -using GetPVSwitch = Dune::Std::detected_or; - // helper struct and function detecting if the linear solver features a norm() function template using NormDetector = decltype(std::declval().norm(std::declval())); @@ -188,10 +182,13 @@ private: using ConvergenceWriter = ConvergenceWriterInterface; using TimeLoop = TimeLoopBase; - using PrimaryVariableSwitch = typename Detail::GetPVSwitch::type; - using HasPriVarsSwitch = typename Detail::GetPVSwitch::value_t; // std::true_type or std::false_type - static constexpr bool hasPriVarsSwitch() { return HasPriVarsSwitch::value; }; + // enable models with primary variable switch + // TODO: Always use ParentType::Variables once we require assemblers to export variables static constexpr bool assemblerExportsVariables = Impl::exportsVariables; + using Vars = std::conditional_t; + using PrimaryVariableSwitchAdapter = Dumux::PrimaryVariableSwitchAdapter; public: using typename ParentType::Variables; @@ -208,6 +205,7 @@ public: , endIterMsgStream_(std::ostringstream::out) , comm_(comm) , paramGroup_(paramGroup) + , priVarSwitchAdapter_(std::make_unique(paramGroup)) { initParams_(paramGroup); @@ -221,12 +219,6 @@ public: // initialize the partial reassembler if (enablePartialReassembly_) partialReassembler_ = std::make_unique(this->assembler()); - - if (hasPriVarsSwitch()) - { - const int priVarSwitchVerbosity = getParamFromGroup(paramGroup, "PrimaryVariableSwitch.Verbosity", 1); - priVarSwitch_ = std::make_unique(priVarSwitchVerbosity); - } } //! the communicator for parallel runs @@ -370,13 +362,20 @@ public: virtual void newtonBegin(Variables& initVars) { numSteps_ = 0; - initPriVarSwitch_(initVars, HasPriVarsSwitch{}); + + if constexpr (assemblerExportsVariables) + priVarSwitchAdapter_->initialize(Backend::getDofVector(initVars), initVars); + else // this assumes assembly with solution (i.e. Variables=SolutionVector) + priVarSwitchAdapter_->initialize(initVars, this->assembler().gridVariables()); + const auto& initSol = Backend::getDofVector(initVars); // write the initial residual if a convergence writer was set if (convergenceWriter_) { this->assembler().assembleResidual(initVars); + + // dummy vector, there is no delta before solving the linear system auto delta = Backend::makeZeroDofVector(Backend::size(initSol)); convergenceWriter_->write(initSol, delta, this->assembler().residual()); } @@ -591,7 +590,10 @@ public: virtual void newtonEndStep(Variables &vars, const SolutionVector &uLastIter) { - invokePriVarSwitch_(vars, HasPriVarsSwitch{}); + if constexpr (assemblerExportsVariables) + priVarSwitchAdapter_->invoke(Backend::getDofVector(vars), vars); + else // this assumes assembly with solution (i.e. Variables=SolutionVector) + priVarSwitchAdapter_->invoke(vars, this->assembler().gridVariables()); ++numSteps_; @@ -633,7 +635,7 @@ public: { // in case the model has a priVar switch and some some primary variables // actually switched their state in the last iteration, enforce another iteration - if (priVarsSwitchedInLastIteration_) + if (priVarSwitchAdapter_->switched()) return false; if (enableShiftCriterion_ && !enableResidualCriterion_) @@ -842,87 +844,6 @@ protected: bool enableResidualCriterion() const { return enableResidualCriterion_; } - /*! - * \brief Initialize the privar switch, noop if there is no priVarSwitch - */ - void initPriVarSwitch_(Variables&, std::false_type) {} - - /*! - * \brief Initialize the privar switch - */ - void initPriVarSwitch_(Variables& vars, std::true_type) - { - priVarSwitch_->reset(Backend::size(Backend::getDofVector(vars))); - priVarsSwitchedInLastIteration_ = false; - - const auto& problem = this->assembler().problem(); - const auto& gridGeometry = this->assembler().gridGeometry(); - - // if the assembler does not export the variables type - // we assume old-style assemblers which operate on solution - // vectors and return the grid variables - if constexpr (!assemblerExportsVariables) - priVarSwitch_->updateDirichletConstraints(problem, gridGeometry, - this->assembler().gridVariables(), - vars); - else - priVarSwitch_->updateDirichletConstraints(problem, gridGeometry, vars, Backend::getDofVector(vars)); - } - - /*! - * \brief Switch primary variables if necessary, noop if there is no priVarSwitch - */ - void invokePriVarSwitch_(Variables&, std::false_type) {} - - /*! - * \brief Switch primary variables if necessary - */ - void invokePriVarSwitch_(Variables& vars, std::true_type) - { - // update the variable switch (returns true if the pri vars at at least one dof were switched) - // for disabled grid variable caching - const auto& gridGeometry = this->assembler().gridGeometry(); - const auto& problem = this->assembler().problem(); - - // we assume an old-style assembler if it doesn't export the Variable type - if constexpr (!assemblerExportsVariables) - { - auto& gridVariables = this->assembler().gridVariables(); - - // invoke the primary variable switch - priVarsSwitchedInLastIteration_ = priVarSwitch_->update(vars, gridVariables, problem, gridGeometry); - - if (priVarsSwitchedInLastIteration_) - { - for (const auto& element : elements(gridGeometry.gridView())) - { - // if the volume variables are cached globally, we need to update those where the primary variables have been switched - priVarSwitch_->updateSwitchedVolVars(problem, element, gridGeometry, gridVariables, Backend::getDofVector(vars)); - - // if the flux variables are cached globally, we need to update those where the primary variables have been switched - priVarSwitch_->updateSwitchedFluxVarsCache(problem, element, gridGeometry, gridVariables, Backend::getDofVector(vars)); - } - } - } - else - { - // invoke the primary variable switch - priVarsSwitchedInLastIteration_ = priVarSwitch_->update(Backend::getDofVector(vars), vars, problem, gridGeometry); - - if (priVarsSwitchedInLastIteration_) - { - for (const auto& element : elements(gridGeometry.gridView())) - { - // if the volume variables are cached globally, we need to update those where the primary variables have been switched - priVarSwitch_->updateSwitchedVolVars(problem, element, gridGeometry, vars, Backend::getDofVector(vars)); - - // if the flux variables are cached globally, we need to update those where the primary variables have been switched - priVarSwitch_->updateSwitchedFluxVarsCache(problem, element, gridGeometry, vars, Backend::getDofVector(vars)); - } - } - } - } - //! optimal number of iterations we want to achieve int targetSteps_; //! minimum number of iterations we do @@ -1383,9 +1304,7 @@ private: std::size_t numLinearSolverBreakdowns_ = 0; //! total number of linear solves that failed //! the class handling the primary variable switch - std::unique_ptr priVarSwitch_; - //! if we switched primary variables in the last iteration - bool priVarsSwitchedInLastIteration_ = false; + std::unique_ptr priVarSwitchAdapter_; //! convergence writer std::shared_ptr convergenceWriter_ = nullptr; diff --git a/dumux/nonlinear/primaryvariableswitchadapter.hh b/dumux/nonlinear/primaryvariableswitchadapter.hh new file mode 100644 index 0000000000..df202314c6 --- /dev/null +++ b/dumux/nonlinear/primaryvariableswitchadapter.hh @@ -0,0 +1,137 @@ +// -*- 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 Nonlinear + * \brief An adapter for the Newton to manage models with primary variable switch + */ +#ifndef DUMUX_NONLINEAR_PRIMARY_VARIABLE_SWITCH_ADAPTER_HH +#define DUMUX_NONLINEAR_PRIMARY_VARIABLE_SWITCH_ADAPTER_HH + +#include +#include + +namespace Dumux { +namespace Impl { + +//! helper aliases to extract a primary variable switch from the VolumeVariables (if defined, yields int otherwise) +template +using DetectPVSwitch = typename Variables::VolumeVariables::PrimaryVariableSwitch; + +template +using GetPVSwitch = Dune::Std::detected_or; + +template +using PrimaryVariableSwitch = typename GetPVSwitch::type; + +template +constexpr bool hasPriVarsSwitch() { return typename GetPVSwitch::value_t(); }; + +} // end namespace Impl + +/*! + * \ingroup Nonlinear + * \brief An adapter for the Newton to manage models with primary variable switch + */ +template ()> +class PrimaryVariableSwitchAdapter +{ + using PrimaryVariableSwitch = typename Impl::PrimaryVariableSwitch; + +public: + PrimaryVariableSwitchAdapter(const std::string& paramGroup = "") + { + const int priVarSwitchVerbosity = getParamFromGroup(paramGroup, "PrimaryVariableSwitch.Verbosity", 1); + priVarSwitch_ = std::make_unique(priVarSwitchVerbosity); + } + + /*! + * \brief Initialize the privar switch + */ + template + void initialize(SolutionVector& sol, Variables& vars) + { + priVarSwitch_->reset(sol.size()); + priVarsSwitchedInLastIteration_ = false; + const auto& problem = vars.curGridVolVars().problem(); + const auto& gridGeometry = problem.gridGeometry(); + priVarSwitch_->updateDirichletConstraints(problem, gridGeometry, vars, sol); + } + + /*! + * \brief Switch primary variables if necessary + */ + template + void invoke(SolutionVector& uCurrentIter, Variables& vars) + { + // update the variable switch (returns true if the pri vars at at least one dof were switched) + // for disabled grid variable caching + const auto& problem = vars.curGridVolVars().problem(); + const auto& gridGeometry = problem.gridGeometry(); + + // invoke the primary variable switch + priVarsSwitchedInLastIteration_ = priVarSwitch_->update(uCurrentIter, vars, problem, gridGeometry); + if (priVarsSwitchedInLastIteration_) + { + for (const auto& element : elements(gridGeometry.gridView())) + { + // if the volume variables are cached globally, we need to update those where the primary variables have been switched + priVarSwitch_->updateSwitchedVolVars(problem, element, gridGeometry, vars, uCurrentIter); + + // if the flux variables are cached globally, we need to update those where the primary variables have been switched + priVarSwitch_->updateSwitchedFluxVarsCache(problem, element, gridGeometry, vars, uCurrentIter); + } + } + } + + /*! + * \brief Whether the primary variables have been switched in the last call to invoke + */ + bool switched() const + { return priVarsSwitchedInLastIteration_; } + +private: + //! the class handling the primary variable switch + std::unique_ptr priVarSwitch_; + //! if we switched primary variables in the last iteration + bool priVarsSwitchedInLastIteration_ = false; +}; + +/*! + * \ingroup Nonlinear + * \brief An empty adapter for the Newton for models without primary variable switch + */ +template +class PrimaryVariableSwitchAdapter +{ +public: + PrimaryVariableSwitchAdapter(const std::string& paramGroup = "") {} + + template + void initialize(SolutionVector&, Variables&) {} + + template + void invoke(SolutionVector&, Variables&) {} + + bool switched() const { return false; } +}; + +} // end namespace Dumux + +#endif -- GitLab From 49961eb0430bfc16a30c28b9b4c586f7e2d7c97c Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Dennis=20Gl=C3=A4ser?= Date: Mon, 8 Feb 2021 12:45:10 +0100 Subject: [PATCH 17/32] [newton] more distinction between new&old assemblers --- dumux/nonlinear/newtonsolver.hh | 52 +++++++++++++++++++-------------- 1 file changed, 30 insertions(+), 22 deletions(-) diff --git a/dumux/nonlinear/newtonsolver.hh b/dumux/nonlinear/newtonsolver.hh index e12661ff3b..902c2f2abb 100644 --- a/dumux/nonlinear/newtonsolver.hh +++ b/dumux/nonlinear/newtonsolver.hh @@ -292,8 +292,11 @@ public: */ void solve(Variables& vars, TimeLoop& timeLoop) override { - if (this->assembler().isStationaryProblem()) - DUNE_THROW(Dune::InvalidStateException, "Using time step control with stationary problem makes no sense!"); + if constexpr (!assemblerExportsVariables) + { + if (this->assembler().isStationaryProblem()) + DUNE_THROW(Dune::InvalidStateException, "Using time step control with stationary problem makes no sense!"); + } // try solving the non-linear system for (std::size_t i = 0; i <= maxTimeStepDivisions_; ++i) @@ -306,28 +309,33 @@ public: else if (!converged && i < maxTimeStepDivisions_) { - // set solution to previous solution - Backend::update(vars, this->assembler().prevSol()); - - // if this is true, we assume old-style assemblers/grid variables - // TODO: reset time step is more efficient as it simply copies - // prevVolVars into curVolVars. Above (in the new style) - // there is an update. We should think about how to achieve - // this more efficiently. - // TODO: Is there a test with time step reductions? It should - // probably be tested if this works properly. - if (!assemblerExportsVariables) - this->assembler().resetTimeStep(Backend::getDofVector(vars)); - - if (verbosity_ >= 1) + if constexpr (assemblerExportsVariables) + DUNE_THROW(Dune::NotImplemented, "Time step reset for new assembly methods"); + else { - const auto dt = timeLoop.timeStepSize(); - std::cout << Fmt::format("Newton solver did not converge with dt = {} seconds. ", dt) - << Fmt::format("Retrying with time step of dt = {} seconds.\n", dt*retryTimeStepReductionFactor_); + // set solution to previous solution + Backend::update(vars, this->assembler().prevSol()); + + // if this is true, we assume old-style assemblers/grid variables + // TODO: reset time step is more efficient as it simply copies + // prevVolVars into curVolVars. Above (in the new style) + // there will probably be an update. We should think about + // how to achieve this efficiently. + // TODO: Is there a test with time step reductions? It should + // probably be tested if this works properly. + if (!assemblerExportsVariables) + this->assembler().resetTimeStep(Backend::getDofVector(vars)); + + if (verbosity_ >= 1) + { + const auto dt = timeLoop.timeStepSize(); + std::cout << Fmt::format("Newton solver did not converge with dt = {} seconds. ", dt) + << Fmt::format("Retrying with time step of dt = {} seconds.\n", dt*retryTimeStepReductionFactor_); + } + + // try again with dt = dt * retryTimeStepReductionFactor_ + timeLoop.setTimeStepSize(timeLoop.timeStepSize() * retryTimeStepReductionFactor_); } - - // try again with dt = dt * retryTimeStepReductionFactor_ - timeLoop.setTimeStepSize(timeLoop.timeStepSize() * retryTimeStepReductionFactor_); } else -- GitLab From 486f4416af70d67cf11515b498e71c3073085e87 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Dennis=20Gl=C3=A4ser?= Date: Fri, 5 Feb 2021 12:52:59 +0100 Subject: [PATCH 18/32] [test][newton] avoid unnecessary mock interfaces --- test/nonlinear/newton/test_newton.cc | 5 ++--- 1 file changed, 2 insertions(+), 3 deletions(-) diff --git a/test/nonlinear/newton/test_newton.cc b/test/nonlinear/newton/test_newton.cc index 3a0d48a0d7..5a08387777 100644 --- a/test/nonlinear/newton/test_newton.cc +++ b/test/nonlinear/newton/test_newton.cc @@ -32,8 +32,9 @@ class MockScalarAssembler { public: using ResidualType = Dune::BlockVector; - using Scalar = double; + using Variables = ResidualType; using JacobianMatrix = double; + using Scalar = double; void setLinearSystem() {} @@ -65,8 +66,6 @@ public: return res_[0]; } - void updateGridVariables(const ResidualType& sol) {} - private: JacobianMatrix jac_; ResidualType res_; -- GitLab From e6730368e37b0c50f98eeaae673ef79348c170f0 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Dennis=20Gl=C3=A4ser?= Date: Fri, 5 Feb 2021 13:35:06 +0100 Subject: [PATCH 19/32] [test][newton] use norm() of linear solver --- test/nonlinear/newton/test_newton.cc | 17 ++++++++++------- 1 file changed, 10 insertions(+), 7 deletions(-) diff --git a/test/nonlinear/newton/test_newton.cc b/test/nonlinear/newton/test_newton.cc index 5a08387777..248314baf7 100644 --- a/test/nonlinear/newton/test_newton.cc +++ b/test/nonlinear/newton/test_newton.cc @@ -2,6 +2,7 @@ #include #include +#include #include #include @@ -60,12 +61,6 @@ public: ResidualType& residual() { return res_; } - double residualNorm(const ResidualType& sol) - { - assembleResidual(sol); - return res_[0]; - } - private: JacobianMatrix jac_; ResidualType res_; @@ -77,11 +72,19 @@ public: void setResidualReduction(double residualReduction) {} template - bool solve (const double& A, Vector& x, const Vector& b) const + bool solve(const double& A, Vector& x, const Vector& b) const { x[0] = b[0]/A; return true; } + + double norm(const Dune::BlockVector& residual) const + { + assert(residual.size() == 1); + + using std::abs; + return abs(residual[0]); + } }; } // end namespace Dumux -- GitLab From b16965c20211cd453984e237bb54fccc899910af Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Dennis=20Gl=C3=A4ser?= Date: Mon, 8 Feb 2021 12:48:23 +0100 Subject: [PATCH 20/32] [test][newton] remove obsolete assembler interface --- test/nonlinear/newton/test_newton.cc | 2 -- 1 file changed, 2 deletions(-) diff --git a/test/nonlinear/newton/test_newton.cc b/test/nonlinear/newton/test_newton.cc index 248314baf7..c1d81048da 100644 --- a/test/nonlinear/newton/test_newton.cc +++ b/test/nonlinear/newton/test_newton.cc @@ -39,8 +39,6 @@ public: void setLinearSystem() {} - bool isStationaryProblem() { return false; } - ResidualType prevSol() { return ResidualType(0.0); } void resetTimeStep(const ResidualType& sol) {} -- GitLab From f16c7769aa56021046e1136294dea678b3696e42 Mon Sep 17 00:00:00 2001 From: "Dennis.Glaeser" Date: Wed, 21 Oct 2020 12:06:32 +0200 Subject: [PATCH 21/32] [test][1p] test assembly with grid variables This implementation should be replaced as soon as the new assembly is available. For now, this tests if NewtonSolver supports the new grid variables-based assembly. --- .../1p/compressible/instationary/assembler.hh | 53 +++++++++++++++ .../instationary/gridvariables.hh | 68 +++++++++++++++++++ .../1p/compressible/instationary/main.cc | 14 ++-- .../compressible/instationary/properties.hh | 23 ++++++- 4 files changed, 148 insertions(+), 10 deletions(-) create mode 100644 test/porousmediumflow/1p/compressible/instationary/assembler.hh create mode 100644 test/porousmediumflow/1p/compressible/instationary/gridvariables.hh diff --git a/test/porousmediumflow/1p/compressible/instationary/assembler.hh b/test/porousmediumflow/1p/compressible/instationary/assembler.hh new file mode 100644 index 0000000000..6eccb4c9c8 --- /dev/null +++ b/test/porousmediumflow/1p/compressible/instationary/assembler.hh @@ -0,0 +1,53 @@ +// -*- 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 OnePTests + * \brief Dummy assembler that fulfills the new experimental assembly interface. + */ +#ifndef DUMUX_COMPRESSIBLE_ONEP_TEST_ASSEMBLER_HH +#define DUMUX_COMPRESSIBLE_ONEP_TEST_ASSEMBLER_HH + +namespace Dumux::OnePCompressibleTest { + +// Custom assembler to test assembly with grid variables, +// fulfilling the foreseen required interface +template +class GridVarsAssembler : public Assembler +{ +public: + using Assembler::Assembler; + using typename Assembler::GridVariables; + using typename Assembler::ResidualType; + + using Variables = GridVariables; + + void assembleJacobianAndResidual(const GridVariables& gridVars) + { Assembler::assembleJacobianAndResidual(gridVars.dofs()); } + + void assembleResidual(const GridVariables& gridVars) + { Assembler::assembleResidual(gridVars.dofs()); } + + auto residualNorm(const GridVariables& gridVars) + { return Assembler::residualNorm(gridVars.dofs()); } +}; + +} // end namespace Dumux::OnePCompressibleTest + +#endif diff --git a/test/porousmediumflow/1p/compressible/instationary/gridvariables.hh b/test/porousmediumflow/1p/compressible/instationary/gridvariables.hh new file mode 100644 index 0000000000..c9234e9c05 --- /dev/null +++ b/test/porousmediumflow/1p/compressible/instationary/gridvariables.hh @@ -0,0 +1,68 @@ +// -*- 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 OnePTests + * \brief Wrapper around the current FVGridVariables to fulfill the layout + * of the new grid variables to test grid variables-based assembly. + */ +#ifndef DUMUX_COMPRESSIBLE_ONEP_TEST_GRID_VARIABLES_HH +#define DUMUX_COMPRESSIBLE_ONEP_TEST_GRID_VARIABLES_HH + +#include + +namespace Dumux::OnePCompressibleTest { + +template +class TestGridVariables +: public BaseGridVariables +, public Dumux::Experimental::GridVariables +{ + using ExperimentalBase = Dumux::Experimental::GridVariables; + +public: + // export some types to avoid ambiguity + using GridGeometry = GG; + using Scalar = typename BaseGridVariables::Scalar; + + template + TestGridVariables(std::shared_ptr problem, + std::shared_ptr gridGeometry, + const SolutionVector& x) + : BaseGridVariables(problem, gridGeometry) + , ExperimentalBase(gridGeometry, x) + { + BaseGridVariables::init(x); + } + + // update to a new solution + void update(const SolutionVector& x) + { + BaseGridVariables::update(x); + ExperimentalBase::update(x); + } + + // overload some functions to avoid ambiguity + decltype(auto) gridGeometry() const + { return ExperimentalBase::gridGeometry(); } +}; + +} // end namespace Dumux::OnePCompressibleTest + +#endif diff --git a/test/porousmediumflow/1p/compressible/instationary/main.cc b/test/porousmediumflow/1p/compressible/instationary/main.cc index d9b60a4b39..f9371770a3 100644 --- a/test/porousmediumflow/1p/compressible/instationary/main.cc +++ b/test/porousmediumflow/1p/compressible/instationary/main.cc @@ -46,6 +46,7 @@ #include #include "properties.hh" +#include "assembler.hh" int main(int argc, char** argv) { @@ -98,8 +99,7 @@ int main(int argc, char** argv) // the grid variables using GridVariables = GetPropType; - auto gridVariables = std::make_shared(problem, gridGeometry); - gridVariables->init(x); + auto gridVariables = std::make_shared(problem, gridGeometry, x); // get some time loop parameters using Scalar = GetPropType; @@ -108,7 +108,7 @@ int main(int argc, char** argv) auto maxDt = getParam("TimeLoop.MaxTimeStepSize"); // intialize the vtk output module - VtkOutputModule vtkWriter(*gridVariables, x, problem->name()); + VtkOutputModule vtkWriter(*gridVariables, gridVariables->dofs(), problem->name()); using VelocityOutput = GetPropType; vtkWriter.addVelocityOutput(std::make_shared(*gridVariables)); using IOFields = GetPropType; @@ -120,7 +120,8 @@ int main(int argc, char** argv) timeLoop->setMaxTimeStepSize(maxDt); // the assembler with time loop for instationary problem - using Assembler = FVAssembler; + using BaseAssembler = FVAssembler; + using Assembler = Dumux::OnePCompressibleTest::GridVarsAssembler; auto assembler = std::make_shared(problem, gridGeometry, gridVariables, timeLoop, xOld); // the linear solver @@ -138,11 +139,10 @@ int main(int argc, char** argv) timeLoop->start(); do { // linearize & solve - nonLinearSolver.solve(x, *timeLoop); + nonLinearSolver.solve(*gridVariables, *timeLoop); // make the new solution the old solution - xOld = x; - gridVariables->advanceTimeStep(); + xOld = gridVariables->dofs(); // advance to the time loop to the next step timeLoop->advanceTimeStep(); diff --git a/test/porousmediumflow/1p/compressible/instationary/properties.hh b/test/porousmediumflow/1p/compressible/instationary/properties.hh index c3e71c6a9b..43e6604f9e 100644 --- a/test/porousmediumflow/1p/compressible/instationary/properties.hh +++ b/test/porousmediumflow/1p/compressible/instationary/properties.hh @@ -27,17 +27,18 @@ #include - #include #include #include #include +#include #include #include "problem.hh" #include "spatialparams.hh" +#include "gridvariables.hh" namespace Dumux::Properties { // create the type tag nodes @@ -76,6 +77,21 @@ public: using type = FluidSystems::OnePLiquid>>; }; +// use the new experimental grid variables as we test the new experimental gridvars-based-assembly here +template +struct GridVariables +{ +private: + using GG = GetPropType; + using GVV = GetPropType; + using GFC = GetPropType; + using Base = Dumux::FVGridVariables; + using X = GetPropType; + +public: + using type = Dumux::OnePCompressibleTest::TestGridVariables; +}; + // Disable caching (for testing purposes) template struct EnableGridVolumeVariablesCache { static constexpr bool value = false; }; @@ -83,6 +99,7 @@ template struct EnableGridFluxVariablesCache { static constexpr bool value = false; }; template struct EnableGridGeometryCache { static constexpr bool value = false; }; -} -#endif \ No newline at end of file +} // end namespace Dumux::Properties + +#endif -- GitLab From 4466e7bc5c2a16311cf74f35a1ffb6a0e9a2dd78 Mon Sep 17 00:00:00 2001 From: "Dennis.Glaeser" Date: Thu, 22 Oct 2020 16:48:14 +0200 Subject: [PATCH 22/32] [newton] make varsbackend newton-agnostic --- dumux/common/CMakeLists.txt | 1 + .../variablesbackend.hh} | 44 ++++++++++--------- dumux/nonlinear/newtonsolver.hh | 5 ++- 3 files changed, 27 insertions(+), 23 deletions(-) rename dumux/{nonlinear/newtonvariablesbackend.hh => common/variablesbackend.hh} (80%) diff --git a/dumux/common/CMakeLists.txt b/dumux/common/CMakeLists.txt index 9d03d03b99..63abd68aa2 100644 --- a/dumux/common/CMakeLists.txt +++ b/dumux/common/CMakeLists.txt @@ -51,4 +51,5 @@ timemanager.hh valgrind.hh variablelengthspline_.hh variables.hh +variablesbackend.hh DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/dumux/common) diff --git a/dumux/nonlinear/newtonvariablesbackend.hh b/dumux/common/variablesbackend.hh similarity index 80% rename from dumux/nonlinear/newtonvariablesbackend.hh rename to dumux/common/variablesbackend.hh index a154671f8c..52e506e4ae 100644 --- a/dumux/nonlinear/newtonvariablesbackend.hh +++ b/dumux/common/variablesbackend.hh @@ -19,10 +19,12 @@ /*! * \file * \ingroup Nonlinear - * \brief Backends for the Newton for different variable types + * \brief Backends for operations on different solution vector types + * or more generic variable classes to be used in places where + * several different types/layouts should be supported. */ -#ifndef DUMUX_NEWTON_VARIABLES_BACKEND_HH -#define DUMUX_NEWTON_VARIABLES_BACKEND_HH +#ifndef DUMUX_VARIABLES_BACKEND_HH +#define DUMUX_VARIABLES_BACKEND_HH #include #include @@ -49,14 +51,14 @@ namespace Dumux { * \brief Class providing operations with primary variable vectors */ template -class NewtonDofBackend; +class DofBackend; /*! * \ingroup Nonlinear - * \brief Specialization providing Newton operations for scalar/number types + * \brief Specialization providing operations for scalar/number types */ template -class NewtonDofBackend::value, Scalar>> +class DofBackend::value, Scalar>> { public: using DofVector = Scalar; //!< the type of the dofs parametrizing the variables object @@ -70,10 +72,10 @@ public: /*! * \ingroup Nonlinear - * \brief Specialization providing Newton operations for block vectors + * \brief Specialization providing operations for block vectors */ template -class NewtonDofBackend> +class DofBackend> { public: @@ -88,10 +90,10 @@ public: /*! * \ingroup Nonlinear - * \brief Specialization providing Newton operations for multitype block vectors + * \brief Specialization providing operations for multitype block vectors */ template -class NewtonDofBackend> +class DofBackend> { using DV = Dune::MultiTypeBlockVector; static constexpr auto numBlocks = DV::size(); @@ -128,19 +130,19 @@ template using SolutionVectorType = typename Vars::SolutionVector; template -class NewtonVariablesBackend; +class VariablesBackend; /*! * \ingroup Nonlinear - * \brief Class providing Newton operations for scalar/number types + * \brief Class providing operations for primary variable vector/scalar types * \note We assume the variables being simply a dof vector if we * do not find the variables class to export `SolutionVector`. */ template -class NewtonVariablesBackend -: public NewtonDofBackend +class VariablesBackend +: public DofBackend { - using ParentType = NewtonDofBackend; + using ParentType = DofBackend; public: using Variables = Vars; @@ -162,12 +164,12 @@ public: /*! * \file * \ingroup Nonlinear - * \brief Class providing Newton operations for generic variable - * types, possibly containing also secondary variables. + * \brief Class providing operations for generic variable classes, + * containing primary and possibly also secondary variables. */ template -class NewtonVariablesBackend -: public NewtonDofBackend +class VariablesBackend +: public DofBackend { public: using DofVector = typename Vars::SolutionVector; @@ -189,13 +191,13 @@ public: /*! * \ingroup Nonlinear - * \brief Class providing Newton operations for generic variable classes + * \brief Class providing operations for generic variable classes * that represent the state of a numerical solution, possibly * consisting of primary/secondary variables and information on * the time level. */ template -using NewtonVariablesBackend = Impl::NewtonVariablesBackend>; +using VariablesBackend = Impl::VariablesBackend>; } // end namespace Dumux diff --git a/dumux/nonlinear/newtonsolver.hh b/dumux/nonlinear/newtonsolver.hh index 902c2f2abb..f4f986a553 100644 --- a/dumux/nonlinear/newtonsolver.hh +++ b/dumux/nonlinear/newtonsolver.hh @@ -48,13 +48,14 @@ #include #include #include +#include + #include #include #include #include #include "newtonconvergencewriter.hh" -#include "newtonvariablesbackend.hh" #include "primaryvariableswitchadapter.hh" namespace Dumux { @@ -173,7 +174,7 @@ class NewtonSolver : public PDESolver using ParentType = PDESolver; protected: - using Backend = NewtonVariablesBackend; + using Backend = VariablesBackend; using SolutionVector = typename Backend::DofVector; private: -- GitLab From 702c760cf1dd19e0157f47a877c898b08802479c Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Dennis=20Gl=C3=A4ser?= Date: Mon, 8 Mar 2021 12:56:20 +0100 Subject: [PATCH 23/32] [disc][fvelemvars] export underlying views --- dumux/discretization/fvgridvariables.hh | 9 ++++++--- 1 file changed, 6 insertions(+), 3 deletions(-) diff --git a/dumux/discretization/fvgridvariables.hh b/dumux/discretization/fvgridvariables.hh index 92656b7d45..5a95f76a96 100644 --- a/dumux/discretization/fvgridvariables.hh +++ b/dumux/discretization/fvgridvariables.hh @@ -192,13 +192,16 @@ namespace Dumux::Experimental { using GridView = typename GridGeometry::GridView; using Element = typename GridView::template Codim<0>::Entity; - using ElementVolumeVariables = typename GV::GridVolumeVariables::LocalView; - using ElementFluxVariablesCache = typename GV::GridFluxVariablesCache::LocalView; - public: //! export corresponding grid-wide class using GridVariables = GV; + //! export underlying volume variables cache + using ElementVolumeVariables = typename GV::GridVolumeVariables::LocalView; + + //! export underlying flux variables cache + using ElementFluxVariablesCache = typename GV::GridFluxVariablesCache::LocalView; + //! Constructor FVGridVariablesLocalView(const GridVariables& gridVariables) : gridVariables_(&gridVariables) -- GitLab From d17c85efc510750bf2bf74313ab7950a5a698a92 Mon Sep 17 00:00:00 2001 From: "Dennis.Glaeser" Date: Thu, 22 Oct 2020 16:59:30 +0200 Subject: [PATCH 24/32] [linear][pdesolver] support variables-based assembly --- dumux/linear/pdesolver.hh | 19 ++++++++++++++----- 1 file changed, 14 insertions(+), 5 deletions(-) diff --git a/dumux/linear/pdesolver.hh b/dumux/linear/pdesolver.hh index 2303b8cbc5..82e7803282 100644 --- a/dumux/linear/pdesolver.hh +++ b/dumux/linear/pdesolver.hh @@ -42,6 +42,8 @@ #include #include #include +#include + #include #include @@ -67,6 +69,8 @@ class LinearPDESolver : public PDESolver using TimeLoop = TimeLoopBase; public: + using typename ParentType::Variables; + /*! * \brief The Constructor */ @@ -86,7 +90,7 @@ public: /*! * \brief Solve a linear PDE system */ - void solve(SolutionVector& uCurrentIter) override + void solve(Variables& vars) override { Dune::Timer assembleTimer(false); Dune::Timer solveTimer(false); @@ -104,9 +108,9 @@ public: // linearize the problem at the current solution assembleTimer.start(); if (reuseMatrix_) - this->assembler().assembleResidual(uCurrentIter); + this->assembler().assembleResidual(vars); else - this->assembler().assembleJacobianAndResidual(uCurrentIter); + this->assembler().assembleJacobianAndResidual(vars); assembleTimer.stop(); /////////////// @@ -128,7 +132,8 @@ public: solveTimer.start(); // set the delta vector to zero before solving the linear system! - SolutionVector deltaU(uCurrentIter); + using Backend = VariablesBackend; + auto deltaU = Backend::getDofVector(vars); deltaU = 0; // solve by calling the appropriate implementation depending on whether the linear solver @@ -150,8 +155,12 @@ public: // update the current solution and secondary variables updateTimer.start(); + + // TODO: This currently does one additional copy in case assembly from solution is used + auto uCurrentIter = Backend::getDofVector(vars); uCurrentIter -= deltaU; - this->assembler().updateGridVariables(uCurrentIter); + Backend::update(vars, uCurrentIter); + updateTimer.stop(); if (verbose_) { -- GitLab From 9c1be091fa92d699402325e14a14608107c5a4a1 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Dennis=20Gl=C3=A4ser?= Date: Mon, 8 Mar 2021 16:00:55 +0100 Subject: [PATCH 25/32] [disc] introduce local context --- dumux/discretization/CMakeLists.txt | 1 + dumux/discretization/localcontext.hh | 124 +++++++++++++++++++++++++++ 2 files changed, 125 insertions(+) create mode 100644 dumux/discretization/localcontext.hh diff --git a/dumux/discretization/CMakeLists.txt b/dumux/discretization/CMakeLists.txt index 3b145a92fb..625fe8b3cc 100644 --- a/dumux/discretization/CMakeLists.txt +++ b/dumux/discretization/CMakeLists.txt @@ -20,6 +20,7 @@ functionspacebasis.hh fvgridvariables.hh fvproperties.hh gridvariables.hh +localcontext.hh localview.hh method.hh rotationpolicy.hh diff --git a/dumux/discretization/localcontext.hh b/dumux/discretization/localcontext.hh new file mode 100644 index 0000000000..285c510748 --- /dev/null +++ b/dumux/discretization/localcontext.hh @@ -0,0 +1,124 @@ +// -*- 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 Discretization + * \brief Class that contains the element-local (or element stencil-local) data + * required to evaluate the terms of discrete equations. + */ +#ifndef DUMUX_LOCAL_CONTEXT_HH +#define DUMUX_LOCAL_CONTEXT_HH + +#include + +namespace Dumux { +namespace Experimental { + +/*! + * \ingroup Discretization + * \brief Implementation of element-stencil-local contexts, which solely store + * the local geometry and primary/secondary variables. This implementation + * defines the minimum interface required for contexts in single-domain + * applications to work with the default assembly mechanism. + * \tparam EV The element-local view on the grid variables + */ +template +class LocalContext +{ + using GridVariables = typename EV::GridVariables; + using GridGeometry = typename GridVariables::GridGeometry; + +public: + using ElementGridGeometry = typename GridGeometry::LocalView; + using ElementVariables = EV; + + //! Constructor + LocalContext(const ElementGridGeometry& eg, + const ElementVariables& ev) + : elemGeom_(&eg) + , elemVars_(&ev) + {} + + //! Return the element-local view on the grid geometry + const ElementGridGeometry elementGridGeometry() const + { return *elemGeom_; } + + //! Return the element-local view on the grid variables + const ElementVariables& elementVariables() const + { return *elemVars_; } + +private: + const ElementGridGeometry* elemGeom_; + const ElementVariables* elemVars_; +}; + +/*! + * \ingroup Discretization + * \brief Implementation of element-stencil-local contexts for multidomain simulations, + * which additionally provide access to coupling data within the local scope. + * \tparam EV The element-local view on the grid variables + * \tparam CD The type containing the local coupling data. + */ +template +class MultiDomainLocalContext +{ + using ParentType = LocalContext +LocalContext +makeLocalContext(const typename EV::GridVariables::GridGeometry::LocalView& gglocalView, + const EV& elemVariables) +{ return {gglocalView, elemVariables}; } + +/*! + * \ingroup Discretization + * \brief Factory function to create a multidomain context from local views. + */ +template +MultiDomainLocalContext +makeMultiDomainLocalContext(const typename EV::GridVariables::GridGeometry::LocalView& gglocalView, + const EV& elemVariables, + const CD& couplingData) +{ return {gglocalView, elemVariables, couplingData}; } + +} // end namespace Experimental +} // end namespace Dumux + +#endif -- GitLab From 3e75815401b0b6b757084cb09db860f56ba0ea0b Mon Sep 17 00:00:00 2001 From: "Dennis.Glaeser" Date: Thu, 22 Oct 2020 17:02:20 +0200 Subject: [PATCH 26/32] [assembly] introduce fv operators --- dumux/assembly/fv/CMakeLists.txt | 3 + dumux/assembly/fv/operators.hh | 125 +++++++++++++++++++++++++++++++ 2 files changed, 128 insertions(+) create mode 100644 dumux/assembly/fv/CMakeLists.txt create mode 100644 dumux/assembly/fv/operators.hh diff --git a/dumux/assembly/fv/CMakeLists.txt b/dumux/assembly/fv/CMakeLists.txt new file mode 100644 index 0000000000..72883ac6ae --- /dev/null +++ b/dumux/assembly/fv/CMakeLists.txt @@ -0,0 +1,3 @@ +install(FILES +operators.hh +DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/dumux/assembly/fv) diff --git a/dumux/assembly/fv/operators.hh b/dumux/assembly/fv/operators.hh new file mode 100644 index 0000000000..e73de3d124 --- /dev/null +++ b/dumux/assembly/fv/operators.hh @@ -0,0 +1,125 @@ +// -*- 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 The base class for the sub-control entity-local evaluation of + * the terms of equations in the context of finite-volume schemes + */ +#ifndef DUMUX_FV_OPERATORS_HH +#define DUMUX_FV_OPERATORS_HH + +#include +#include + +namespace Dumux::Experimental { + +/*! + * \ingroup Assembly + * \brief The base class for the sub-control entity-local evaluation of + * the terms of equations in the context of finite-volume schemes + * \tparam Context the element-stencil-local data required to evaluate the terms + */ +template +class FVOperators +{ + // The grid geometry on which the scheme operates + using FVElementGeometry = typename LocalContext::ElementGridGeometry; + using GridGeometry = typename GridVariables::GridGeometry; + using SubControlVolume = typename GridGeometry::SubControlVolume; + using SubControlVolumeFace = typename GridGeometry::SubControlVolumeFace; + using Extrusion = Extrusion_t; + +public: + //! export the type used to store scalar values for all equations + using NumEqVector = Dumux::NumEqVector; + + // export the types of the flux/source/storage terms + using FluxTerm = NumEqVector; + using SourceTerm = NumEqVector; + using StorageTerm = NumEqVector; + + /*! + * \name Model specific interfaces + * \note The following method are the model specific implementations of the + * operators appearing in the equation and should be overloaded by the + * implementations. + */ + // \{ + + /*! + * \brief Compute the storage term of the equations for the given sub-control volume + * \param problem The problem to be solved (could store additionally required quantities) + * \param context The element-stencil-local required data + * \param scv The sub-control volume + * \note This must be overloaded by the implementation + */ + template + static StorageTerm storage(const Problem& problem, + const LocalContext& context, + const SubControlVolume& scv) + { DUNE_THROW(Dune::NotImplemented, "Storage operator not implemented!"); } + + /*! + * \brief Compute the flux term of the equations for the given sub-control volume face + * \param problem The problem to be solved (could store additionally required quantities) + * \param context The element-stencil-local required data + * \param scvf The sub-control volume face for which the flux term is to be computed + * \note This must be overloaded by the implementation + */ + template + static FluxTerm flux(const Problem& problem, + const LocalContext& context, + const SubControlVolumeFace& scvf) + { DUNE_THROW(Dune::NotImplemented, "This model does not implement a flux method!"); } + + /*! + * \brief Compute the source term of the equations for the given sub-control volume + * \param problem The problem to be solved (could store additionally required quantities) + * \param context The element-stencil-local required data + * \param scv The sub-control volume for which the source term is to be computed + * \note This is a default implementation forwarding to interfaces in the problem + */ + template + static SourceTerm source(const Problem& problem, + const LocalContext& context, + const SubControlVolume& scv) + { + SourceTerm source(0.0); + + // TODO: The problem-interfaces should be adapted to receive context? + const auto& fvGeometry = context.elementGridGeometry(); + const auto& elemVolVars = context.elementVariables().elemVolVars(); + + // add contributions from volume flux sources + source += problem.source(element, fvGeometry, elemVolVars, scv); + + // add contribution from possible point sources + source += problem.scvPointSources(element, fvGeometry, elemVolVars, scv); + + // multiply with scv volume + source *= Extrusion::volume(scv)*elemVolVars[scv].extrusionFactor(); + + return source; + } +}; + +} // end namespace Dumux::Experimental + +#endif -- GitLab From 3e48c0e0259462497f909d6d75aa5f0e9989e33f Mon Sep 17 00:00:00 2001 From: "Dennis.Glaeser" Date: Thu, 22 Oct 2020 17:07:19 +0200 Subject: [PATCH 27/32] [nonisothermal] add fv-type operators --- .../nonisothermal/CMakeLists.txt | 1 + .../nonisothermal/operators.hh | 154 ++++++++++++++++++ 2 files changed, 155 insertions(+) create mode 100644 dumux/porousmediumflow/nonisothermal/operators.hh diff --git a/dumux/porousmediumflow/nonisothermal/CMakeLists.txt b/dumux/porousmediumflow/nonisothermal/CMakeLists.txt index 24dc71bd45..5db30722a0 100644 --- a/dumux/porousmediumflow/nonisothermal/CMakeLists.txt +++ b/dumux/porousmediumflow/nonisothermal/CMakeLists.txt @@ -3,5 +3,6 @@ indices.hh iofields.hh localresidual.hh model.hh +operators.hh volumevariables.hh DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/dumux/porousmediumflow/nonisothermal) diff --git a/dumux/porousmediumflow/nonisothermal/operators.hh b/dumux/porousmediumflow/nonisothermal/operators.hh new file mode 100644 index 0000000000..1949e41397 --- /dev/null +++ b/dumux/porousmediumflow/nonisothermal/operators.hh @@ -0,0 +1,154 @@ +// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- +// vi: set et ts=4 sw=4 sts=4: +/***************************************************************************** + * 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 NIModel + * \brief Sub-control entity-local evaluation of the operators + * involved in the system of equations non-isothermal models + * that assume thermal equilibrium between all phases. + */ +#ifndef DUMUX_FV_NON_ISOTHERMAL_OPERATORS_HH +#define DUMUX_FV_NON_ISOTHERMAL_OPERATORS_HH + +namespace Dumux::Experimental { + +/*! + * \ingroup NIModel + * \brief Sub-control entity-local evaluation of the operators + * involved in the system of equations of non-isothermal models + * that assume thermal equilibrium between all phases. + * \tparam ModelTraits Model-specific traits. + * \tparam LocalContext Element-local context (geometry & primary/secondary variables) + */ +template +class FVNonIsothermalOperators +{ + // The variables required for the evaluation of the equation + using ElementVariables = typename LocalContext::ElementVariables; + + // The grid geometry on which the scheme operates + using GridGeometry = typename LocalContext::ElementGridGeometry::GridGeometry; + using SubControlVolume = typename GridGeometry::SubControlVolume; + +public: + + /*! + * \brief The energy storage in the fluid phase with index phaseIdx + * + * \param storage The mass of the component within the sub-control volume + * \param scv The sub-control volume + * \param context The element-local context (primary/secondary variables) + * \param phaseIdx The phase index + */ + template + static void addFluidPhaseStorage(NumEqVector& storage, + const SubControlVolume& scv, + const LocalContext& context, + int phaseIdx) + { + // do no-op in case energy balance is not active + if constexpr (ModelTraits::enableEnergyBalance()) + { + const auto& volVars = context.elementVariables().elemVolVars()[scv]; + storage[ModelTraits::Indices::energyEqIdx] += volVars.porosity() + *volVars.density(phaseIdx) + *volVars.internalEnergy(phaseIdx) + *volVars.saturation(phaseIdx); + } + } + + /*! + * \brief The energy storage in the solid matrix + * + * \param storage The mass of the component within the sub-control volume + * \param scv The sub-control volume + * \param context The element-local context (primary/secondary variables) + */ + template + static void addSolidPhaseStorage(NumEqVector& storage, + const SubControlVolume& scv, + const LocalContext& context) + { + // do no-op in case energy balance is not active + if constexpr (ModelTraits::enableEnergyBalance()) + { + const auto& volVars = context.elementVariables().elemVolVars()[scv]; + storage[ModelTraits::Indices::energyEqIdx] += volVars.temperature() + *volVars.solidHeatCapacity() + *volVars.solidDensity() + *(1.0 - volVars.porosity()); + } + } + + /*! + * \brief The advective phase energy fluxes + * + * \param flux The flux + * \param fluxVars The flux variables. + * \param phaseIdx The phase index + */ + template + static void addHeatConvectionFlux(NumEqVector& flux, + FluxVariables& fluxVars, + int phaseIdx) + { + // do no-op in case energy balance is not active + if constexpr (ModelTraits::enableEnergyBalance()) + { + auto upwindTerm = [phaseIdx](const auto& volVars) + { return volVars.density(phaseIdx)*volVars.mobility(phaseIdx)*volVars.enthalpy(phaseIdx); }; + + flux[ModelTraits::Indices::energyEqIdx] += fluxVars.advectiveFlux(phaseIdx, upwindTerm); + } + } + + /*! + * \brief The diffusive energy fluxes + * + * \param flux The flux + * \param fluxVars The flux variables. + */ + template + static void addHeatConductionFlux(NumEqVector& flux, + FluxVariables& fluxVars) + { + // do no-op in case energy balance is not active + if constexpr (ModelTraits::enableEnergyBalance()) + flux[ModelTraits::Indices::energyEqIdx] += fluxVars.heatConductionFlux(); + } + + /*! + * \brief heat transfer between the phases for nonequilibrium models + * + * \param source The source which ought to be simulated + * \param element An element which contains part of the control volume + * \param fvGeometry The finite-volume geometry + * \param context The element-local context (primary/secondary variables) + * \param scv The sub-control volume over which we integrate the source term + */ + template + static void addSourceEnergy(NumEqVector& source, + const LocalContext& context, + const SubControlVolume &scv) + {} +}; + +} // end namespace Dumux::Experimental + +#endif -- GitLab From e9cc4cecb4843e56e10c674b337b1af6ee9b06f8 Mon Sep 17 00:00:00 2001 From: "Dennis.Glaeser" Date: Thu, 22 Oct 2020 17:09:16 +0200 Subject: [PATCH 28/32] [immiscible] add fv-type operators --- .../immiscible/CMakeLists.txt | 1 + .../porousmediumflow/immiscible/operators.hh | 169 ++++++++++++++++++ 2 files changed, 170 insertions(+) create mode 100644 dumux/porousmediumflow/immiscible/operators.hh diff --git a/dumux/porousmediumflow/immiscible/CMakeLists.txt b/dumux/porousmediumflow/immiscible/CMakeLists.txt index 192e7955c6..ccb0268dc6 100644 --- a/dumux/porousmediumflow/immiscible/CMakeLists.txt +++ b/dumux/porousmediumflow/immiscible/CMakeLists.txt @@ -1,3 +1,4 @@ install(FILES localresidual.hh +operators.hh DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/dumux/porousmediumflow/immiscible) diff --git a/dumux/porousmediumflow/immiscible/operators.hh b/dumux/porousmediumflow/immiscible/operators.hh new file mode 100644 index 0000000000..33337e79ad --- /dev/null +++ b/dumux/porousmediumflow/immiscible/operators.hh @@ -0,0 +1,169 @@ +// -*- 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 PorousmediumflowModels + * \brief Sub-control entity-local evaluation of the operators + * within in the systems of equations of n-phase immiscible models. + */ +#ifndef DUMUX_FV_IMMISCIBLE_OPERATORS_HH +#define DUMUX_FV_IMMISCIBLE_OPERATORS_HH + +#include +#include +#include + +namespace Dumux::Experimental { + +/*! + * \ingroup PorousmediumflowModels + * \brief Sub-control entity-local evaluation of the operators + * within in the systems of equations of n-phase immiscible models. + * \tparam ModelTraits defines model-related types and variables (e.g. number of phases) + * \tparam FluxVariables the type that is responsible for computing the individual + * flux contributions, i.e., advective, diffusive, convective... + * \tparam LocalContext the element-stencil-local data required to evaluate the terms + * \tparam EnergyOperators optional template argument, specifying the class that + * handles the operators related to non-isothermal effects. + * The default energy operators are compatible with isothermal + * simulations. + */ +template> +class FVImmiscibleOperators +: public FVOperators +{ + using ParentType = FVOperators; + + // The variables required for the evaluation of the equation + using + using GridVariables = typename ElementVariables::GridVariables; + using VolumeVariables = typename GridVariables::VolumeVariables; + using ElementVolumeVariables = typename ElementVariables::ElementVolumeVariables; + using ElementFluxVariablesCache = typename ElementVariables::ElementFluxVariablesCache; + + // The grid geometry on which the scheme operates + using GridGeometry = typename GridVariables::GridGeometry; + using FVElementGeometry = typename GridGeometry::LocalView; + using SubControlVolume = typename GridGeometry::SubControlVolume; + using SubControlVolumeFace = typename GridGeometry::SubControlVolumeFace; + using Extrusion = Extrusion_t; + + using GridView = typename GridGeometry::GridView; + using Element = typename GridView::template Codim<0>::Entity; + + using Problem = std::decay_t().gridVolVars().problem())>; + using Indices = typename ModelTraits::Indices; + + static constexpr int numPhases = ModelTraits::numFluidPhases(); + static constexpr int conti0EqIdx = ModelTraits::Indices::conti0EqIdx; + static constexpr bool enableEnergyBalance = ModelTraits::enableEnergyBalance();; + +public: + //! export the type used to store scalar values for all equations + using typename ParentType::NumEqVector; + + // export the types of the flux/source/storage terms + using typename ParentType::FluxTerm; + using typename ParentType::SourceTerm; + using typename ParentType::StorageTerm; + + /*! + * \brief Compute the storage term of the equations for the given sub-control volume + * \param problem The problem to be solved (could store additionally required quantities) + * \param context The element-stencil-local required data + * \param scv The sub-control volume + */ + template + static StorageTerm storage(const Problem& problem, + const LocalContext& context, + const SubControlVolume& scv) + { + const auto& volVars = context.elementVariables().elemVolVars()[scv]; + + StorageTerm storage; + for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) + { + auto eqIdx = conti0EqIdx + phaseIdx; + storage[eqIdx] = volVars.porosity() + * volVars.density(phaseIdx) + * volVars.saturation(phaseIdx); + + // The energy storage in the fluid phase with index phaseIdx + if constexpr (enableEnergyBalance) + EnergyOperators::fluidPhaseStorage(storage, scv, volVars, phaseIdx); + } + + // The energy storage in the solid matrix + if constexpr (enableEnergyBalance) + EnergyOperators::solidPhaseStorage(storage, scv, volVars); + + // multiply with volume + storage *= Extrusion::volume(scv)*volVars.extrusionFactor(); + + return storage; + } + + /*! + * \brief Compute the flux term of the equations for the given sub-control volume face + * \param problem The problem to be solved (could store additionally required quantities) + * \param context The element-stencil-local required data + * \param scvf The sub-control volume face for which the flux term is to be computed + * \note This must be overloaded by the implementation + */ + template + static FluxTerm flux(const Problem& problem, + const LocalContext& context, + const SubControlVolumeFace& scvf) + { + const auto& fvGeometry = context.elementGridGeometry(); + const auto& elemVolVars = context.elementVariables().elemVolVars(); + const auto& elemFluxVarsCache = context.elementVariables().elemFluxVarsCache(); + + FluxVariables fluxVars; + fluxVars.init(problem, fvGeometry.element(), fvGeometry, elemVolVars, scvf, elemFluxVarsCache); + + NumEqVector flux; + for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) + { + // the physical quantities for which we perform upwinding + auto upwindTerm = [phaseIdx](const auto& volVars) + { return volVars.density(phaseIdx)*volVars.mobility(phaseIdx); }; + + auto eqIdx = conti0EqIdx + phaseIdx; + flux[eqIdx] = fluxVars.advectiveFlux(phaseIdx, upwindTerm); + + // Add advective phase energy fluxes. For isothermal model the contribution is zero. + if constexpr (enableEnergyBalance) + EnergyOperators::heatConvectionFlux(flux, fluxVars, phaseIdx); + } + + // Add diffusive energy fluxes. For isothermal model the contribution is zero. + if constexpr (enableEnergyBalance) + EnergyOperators::heatConductionFlux(flux, fluxVars); + + return flux; + } +}; + +} // end namespace Dumux::Experimental + +#endif -- GitLab From 0e640c91b94290395e060a6d5048f2c5a676ddea Mon Sep 17 00:00:00 2001 From: "Dennis.Glaeser" Date: Thu, 22 Oct 2020 17:11:54 +0200 Subject: [PATCH 29/32] [assembly] add fv local operator --- dumux/assembly/fv/CMakeLists.txt | 1 + dumux/assembly/fv/localoperator.hh | 231 +++++++++++++++++++++++++++++ 2 files changed, 232 insertions(+) create mode 100644 dumux/assembly/fv/localoperator.hh diff --git a/dumux/assembly/fv/CMakeLists.txt b/dumux/assembly/fv/CMakeLists.txt index 72883ac6ae..fb70764d7e 100644 --- a/dumux/assembly/fv/CMakeLists.txt +++ b/dumux/assembly/fv/CMakeLists.txt @@ -1,3 +1,4 @@ install(FILES +localoperator.hh operators.hh DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/dumux/assembly/fv) diff --git a/dumux/assembly/fv/localoperator.hh b/dumux/assembly/fv/localoperator.hh new file mode 100644 index 0000000000..33dd7e08c1 --- /dev/null +++ b/dumux/assembly/fv/localoperator.hh @@ -0,0 +1,231 @@ +// -*- 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 An element-wise local operator for finite-volume schemes. + */ +#ifndef DUMUX_FV_LOCAL_OPERATOR_HH +#define DUMUX_FV_LOCAL_OPERATOR_HH + +#include + +#include +#include + +#include +#include +#include + +namespace Dumux::Experimental { + +/*! + * \ingroup Assembly + * \brief The element-wise local operator for finite volume schemes. + * This allows for element-wise evaluation of individual terms + * of the equations to be solved. + * \tparam LC element local data (geometry & primary/secondary variables) + * \tparam OP The model-specific operators + */ +template +class FVLocalOperator +{ + using ElementVariables = typename LC::ElementVariables; + using GridVars = typename ElementVariables::GridVariables; + using PrimaryVariables = typename GridVars::PrimaryVariables; + + using FVElementGeometry = typename LC::ElementGridGeometry; + using SubControlVolumeFace = typename GridGeometry::SubControlVolumeFace; + using Extrusion = Extrusion_t; + + static constexpr int numEq = NumEqVectorTraits::numEq(); + static constexpr bool isBox = GridGeometry::discMethod == DiscretizationMethod::box; + +public: + //! export the expected local context type + using LocalContext = LC; + + //! export the underlying operators + using Operators = OP; + + //! vector type storing the operator values on all dofs of an element + //! TODO: Use ReservedBlockVector + using ElementResidualVector = Dune::BlockVector; + + //! Constructor from a local context + explicit FVLocalOperator(const LocalContext& context) + : context_(context) + {} + + /*! + * \name Main interface + * \note Methods used by the assembler to compute derivatives and residual + */ + // \{ + + /*! + * \brief Compute the terms of the local residual that do not appear in + * time derivatives. These are the sources and the fluxes. + */ + ElementResidualVector evalFluxesAndSources() const + { + auto result = getEmptyResidual(); + const auto& problem = elemVariables_.gridVariables().gridVolVars().problem(); + + // source term + for (const auto& scv : scvs(fvGeometry_)) + result[scv.localDofIndex()] -= Operators::source(problem, context_, scv); + + // flux term + for (const auto& scvf : scvfs(fvGeometry_)) + addFlux_(result, scvf); + + return result; + } + + /*! + * \brief Compute the storage term, i.e. the term appearing in the time derivative. + */ + ElementResidualVector evalStorage() const + { + const auto& problem = elemVariables_.gridVariables().gridVolVars().problem(); + + // TODO: Until now, FVLocalResidual defined storage as the entire + // time derivative. Now it means the term above the time derivative. + // We should think about the correct naming here... + // TODO: Should storage() NOT multiply with volume?? That was different until + // now but flux() also returns the flux multiplied with area so this should + // be more consistent + auto result = getEmptyResidual(); + for (const auto& scv : scvs(fvGeometry_)) + result[scv.localDofIndex()] += operators_.storage(problem, context_ scv); + + return result; + } + + ElementResidualVector getEmptyResidual() const + { + ElementResidualVector res(fvGeometry_.numScv()); + res = 0.0; + return res; + } + + // \} + + /*! + * \name Interfaces for analytic Jacobian computation + */ + // \{ + + //! \todo TODO: Add interfaces. Or, should this be here at all!? + + //\} + + // \} + +protected: + + //! compute and add the flux across the given face to the container (cc schemes) + template = 0> + void addFlux_(ElementResidualVector& r, const SubControlVolumeFace& scvf) const + { + const auto& evv = context_.elementVariables().elemVolVars(); + const auto& problem = evv.gridVolVars().problem(); + + const auto& insideScv = fvGeometry_.scv(scvf.insideScvIdx()); + const auto localDofIdx = insideScv.localDofIndex(); + + // TODO: Modify problem interfaces to receive context + const auto& fvGeometry = context_.elementGridGeometry(); + const auto& element = fvGeometry.element(); + const auto& efvc = context_.elementVariables().elemFluxVarsCache(); + + if (!scvf.boundary()) + r[localDofIdx] += Operators::flux(problem, context, scvf); + else + { + const auto& bcTypes = problem.boundaryTypes(element, scvf); + + // Dirichlet boundaries + if (bcTypes.hasDirichlet() && !bcTypes.hasNeumann()) + r[localDofIdx] += Operators::flux(problem, context, scvf); + + // Neumann and Robin ("solution dependent Neumann") boundary conditions + else if (bcTypes.hasNeumann() && !bcTypes.hasDirichlet()) + { + auto neumannFluxes = problem.neumann(element, fvGeometry, evv, efvc, scvf); + + // multiply neumann fluxes with the area and the extrusion factor + neumannFluxes *= Extrusion::area(scvf)*evv[insideScv].extrusionFactor(); + r[localDofIdx] += neumannFluxes; + } + + else + DUNE_THROW(Dune::NotImplemented, "Mixed boundary conditions for cell-centered schemes. " << + "Use pure boundary conditions by converting Dirichlet BCs to Robin BCs"); + } + } + + //! compute and add the flux across the given face to the container (box scheme) + template = 0> + void addFlux_(ElementResidualVector& r, const SubControlVolumeFace& scvf) const + { + const auto& evv = context_.elementVariables().elemVolVars(); + const auto& problem = evv.gridVolVars().problem(); + + // TODO: Modify problem interfaces to receive context + const auto& fvGeometry = context_.elementGridGeometry(); + const auto& element = fvGeometry.element(); + const auto& efvc = context_.elementVariables().elemFluxVarsCache(); + + // inner faces + if (!scvf.boundary()) + { + const auto flux = Operators::flux(problem, element, fvGeometry, evv, efvc, scvf); + r[fvGeometry_.scv(scvf.insideScvIdx()).localDofIndex()] += flux; + r[fvGeometry_.scv(scvf.outsideScvIdx()).localDofIndex()] -= flux; + } + + // boundary faces + else + { + const auto& scv = fvGeometry_.scv(scvf.insideScvIdx()); + const auto& bcTypes = problem.boundaryTypes(element, scv); + + // Treat Neumann and Robin ("solution dependent Neumann") boundary conditions. + if (bcTypes.hasNeumann()) + { + const auto neumannFluxes = problem.neumann(element, fvGeometry, evv, efvc, scvf); + const auto area = Extrusion::area(scvf)*evv[scv].extrusionFactor(); + + // only add fluxes to equations for which Neumann is set + for (int eqIdx = 0; eqIdx < NumEqVector::size(); ++eqIdx) + if (bcTypes.isNeumann(eqIdx)) + r[scv.localDofIndex()][eqIdx] += neumannFluxes[eqIdx]*area; + } + } + } + +private: + const LocalContext& context_; +}; + +} // end namespace Dumux::Experimental + +#endif -- GitLab From 4268a4c51a43a0430d4d7f8c402d0a0f45ffd67e Mon Sep 17 00:00:00 2001 From: "Dennis.Glaeser" Date: Thu, 22 Oct 2020 17:12:51 +0200 Subject: [PATCH 30/32] [assembly] add new box local assembler --- dumux/assembly/fv/CMakeLists.txt | 1 + dumux/assembly/fv/boxlocalassembler.hh | 455 +++++++++++++++++++++++++ 2 files changed, 456 insertions(+) create mode 100644 dumux/assembly/fv/boxlocalassembler.hh diff --git a/dumux/assembly/fv/CMakeLists.txt b/dumux/assembly/fv/CMakeLists.txt index fb70764d7e..40d6910f65 100644 --- a/dumux/assembly/fv/CMakeLists.txt +++ b/dumux/assembly/fv/CMakeLists.txt @@ -1,4 +1,5 @@ install(FILES +boxlocalassembler.hh localoperator.hh operators.hh DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/dumux/assembly/fv) diff --git a/dumux/assembly/fv/boxlocalassembler.hh b/dumux/assembly/fv/boxlocalassembler.hh new file mode 100644 index 0000000000..92b13ec36d --- /dev/null +++ b/dumux/assembly/fv/boxlocalassembler.hh @@ -0,0 +1,455 @@ +// -*- 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 An assembler for Jacobian and residual contribution + * per element for the box scheme. + */ +#ifndef DUMUX_BOX_LOCAL_ASSEMBLER_HH +#define DUMUX_BOX_LOCAL_ASSEMBLER_HH + +#include +#include + +#include + +#include +#include + +#include +#include +#include + + +namespace Dumux::Experimental { + +/*! + * \ingroup Assembly + * \brief An assembler for Jacobian and residual contribution + * per element for the box scheme. + * \tparam Assembler The grid-wide assembler type + */ +template +class BoxLocalAssembler +{ + using LocalContext = typename LocalOperator::LocalContext; + using ElementVariables = typename LocalContext::ElementVariables; + using ElementGridGeometry = typename LocalContext::ElementGridGeometry; + + using GG = typename ElementGridGeometry::GridGeometry; + using GV = typename ElementVariables::GridVariables; + + using Element = typename GG::GridView::template Codim<0>::Entity; + using PrimaryVariables = typename GV::PrimaryVariables; + using Scalar = typename GV::Scalar; + + static constexpr int numEq = NumEqVectorTraits::numEq(); + +public: + + //! the parameters of a stage in time integration + using StageParams = MultiStageParams; + + //! export the grid variables type on which to operate + using GridVariables = typename ElementVariables::GridVariables; + + /*! + * \brief Constructor for stationary problems. + */ + explicit BoxLocalAssembler(const Element& element, + GridVariables& gridVariables, + DiffMethod dm = DiffMethod::numeric) + : diffMethod_(dm) + , gridVariables_(gridVariables) + , fvGeometry_(localView(gridVariables.gridGeometry())) + , elemVariables_(localView(gridVariables)) + , prevElemVariables_(0) + , elementIsGhost_(element.partitionType() == Dune::GhostEntity) + , stageParams_(nullptr) + { + if (diffMethod_ != DiffMethod::numeric) + DUNE_THROW(Dune::NotImplemented, "Provided differentiation method"); + + fvGeometry_.bind(element); + elemVariables_.bind(element, fvGeometry_); + } + + /*! + * \brief Constructor for instationary problems. + * \note Using this constructor, we assemble one stage within + * a time integration step using multi-stage methods. + */ + explicit BoxLocalAssembler(const Element& element, + GridVariables& gridVariables, + std::vector& prevGridVars, + std::shared_ptr stageParams, + DiffMethod dm = DiffMethod::numeric) + : diffMethod_(dm) + , gridVariables_(gridVariables) + , fvGeometry_(localView(gridVariables.gridGeometry())) + , elemVariables_(localView(gridVariables)) + , prevElemVariables_() + , elementIsGhost_(element.partitionType() == Dune::GhostEntity) + , stageParams_(stageParams) + { + if (diffMethod_ != DiffMethod::numeric) + DUNE_THROW(Dune::NotImplemented, "Provided differentiation method"); + + if (prevGridVars.size() != stageParams.size() - 1) + DUNE_THROW(Dune::InvalidStateException, "Grid Variables for all stages needed"); + + fvGeometry_.bind(element); + for (const auto& prevGV : prevGridVars) + { + prevElemVariables_.emplace_back(localView(prevGV)); + prevElemVariables_.bind(element, fvGeometry_); + } + + elemVariables_.bind(element, fvGeometry_); + } + + /*! + * \brief Computes the derivatives with respect to the given element and adds + * them to the global matrix. The element residual is written into the + * right hand side. + */ + template + void assembleJacobianAndResidual(JacobianMatrix& jac, + ResidualVector& res, + const PartialReassembler* partialReassembler = nullptr) + { + const auto eIdxGlobal = fvGeometry_().gridGeometry().elementMapper().index(element()); + + if (partialReassembler && partialReassembler->elementColor(eIdxGlobal) == EntityColor::green) + { + const auto residual = evalLocalResidual(); + for (const auto& scv : scvs(fvGeometry_())) + res[scv.dofIndex()] += residual[scv.localDofIndex()]; + } + else if (!elementIsGhost_) + { + const auto residual = assembleJacobianAndResidual_(jac, partialReassembler); + for (const auto& scv : scvs(fvGeometry_())) + res[scv.dofIndex()] += residual[scv.localDofIndex()]; + } + else + { + static constexpr int dim = GG::GridView::dimension; + + for (const auto& scv : scvs(fvGeometry_())) + { + const auto vIdxLocal = scv.indexInElement(); + const auto& v = fvGeometry_().element().template subEntity(vIdxLocal); + + // do not change the non-ghost vertices + if (v.partitionType() == Dune::InteriorEntity || + v.partitionType() == Dune::BorderEntity) + continue; + + const auto dofIdx = scv.dofIndex(); + res[dofIdx] = 0; + for (int pvIdx = 0; pvIdx < numEq; ++pvIdx) + jac[dofIdx][dofIdx][pvIdx][pvIdx] = 1.0; + } + } + + auto applyDirichlet = [&] (const auto& scvI, + const auto& dirichletValues, + const auto eqIdx, + const auto pvIdx) + { + res[scvI.dofIndex()][eqIdx] = elemVariables_().elemVolVars()[scvI].priVars()[pvIdx] - dirichletValues[pvIdx]; + + auto& row = jac[scvI.dofIndex()]; + for (auto col = row.begin(); col != row.end(); ++col) + row[col.index()][eqIdx] = 0.0; + + jac[scvI.dofIndex()][scvI.dofIndex()][eqIdx][pvIdx] = 1.0; + + // if a periodic dof has Dirichlet values also apply the same Dirichlet values to the other dof + if (fvGeometry_().gridGeometry().dofOnPeriodicBoundary(scvI.dofIndex())) + { + const auto periodicDof = fvGeometry_().gridGeometry().periodicallyMappedDof(scvI.dofIndex()); + res[periodicDof][eqIdx] = elemVariables_().elemVolVars()[scvI].priVars()[pvIdx] - dirichletValues[pvIdx]; + const auto end = jac[periodicDof].end(); + for (auto it = jac[periodicDof].begin(); it != end; ++it) + (*it) = periodicDof != it.index() ? 0.0 : 1.0; + } + }; + + enforceDirichletConstraints(applyDirichlet); + } + + /*! + * \brief Computes the derivatives with respect to the given element and adds them + * to the global matrix. + */ + template + void assembleJacobian(JacobianMatrix& jac) + { + assembleJacobianAndResidual_(jac); + + auto applyDirichlet = [&] (const auto& scvI, + const auto& dirichletValues, + const auto eqIdx, + const auto pvIdx) + { + auto& row = jac[scvI.dofIndex()]; + for (auto col = row.begin(); col != row.end(); ++col) + row[col.index()][eqIdx] = 0.0; + + jac[scvI.dofIndex()][scvI.dofIndex()][eqIdx][pvIdx] = 1.0; + }; + + enforceDirichletConstraints(applyDirichlet); + } + + /*! + * \brief Assemble the residual only + */ + template + void assembleResidual(ResidualVector& res) + { + const auto residual = evalLocalResidual(); + for (const auto& scv : scvs(fvGeometry_())) + res[scv.dofIndex()] += residual[scv.localDofIndex()]; + + auto applyDirichlet = [&] (const auto& scvI, + const auto& dirichletValues, + const auto eqIdx, + const auto pvIdx) + { + res[scvI.dofIndex()][eqIdx] = elemVariables_().elemVolVars()[scvI].priVars()[pvIdx] - dirichletValues[pvIdx]; + }; + + enforceDirichletConstraints(applyDirichlet); + } + + //! Enforce Dirichlet constraints + template + void enforceDirichletConstraints(const ApplyFunction& applyDirichlet) + { + // enforce Dirichlet boundary conditions + evalDirichletBoundaries(applyDirichlet); + // TODO: take care of internal Dirichlet constraints (if enabled) + // enforceInternalDirichletConstraints(applyDirichlet); + } + + /*! + * \brief Evaluates Dirichlet boundaries + */ + template< typename ApplyDirichletFunctionType > + void evalDirichletBoundaries(ApplyDirichletFunctionType applyDirichlet) + { + for (const auto& scvI : scvs(fvGeometry())) + { + const auto bcTypes = problem().boundaryTypes(element(), scvI); + if (bcTypes.hasDirichlet()) + { + const auto dirichletValues = problem().dirichlet(element(), scvI); + + // set the Dirichlet conditions in residual and jacobian + for (int eqIdx = 0; eqIdx < numEq; ++eqIdx) + { + if (bcTypes.isDirichlet(eqIdx)) + { + const auto pvIdx = bcTypes.eqToDirichletIndex(eqIdx); + assert(0 <= pvIdx && pvIdx < numEq); + applyDirichlet(scvI, dirichletValues, eqIdx, pvIdx); + } + } + } + } + } + + /*! + * \brief Evaluate the complete local residual for the current element. + */ + ElementResidualVector evalLocalResidual() const + { + if (isStationary()) + { + const auto context = makeLocalContext(fvGeometry_, elemVariables_); + LocalOperator localOperator(context); + return elementIsGhost_ ? localOperator.getEmptyResidual() + : localOperator.evalFluxesAndSources(); + } + else + { + ElementResidualVector residual(fvGeometry_().numScv()); + residual = 0.0; + + if (!elementIsGhost_) + { + for (std::size_t k = 0; k < stageParams_->size()-1; ++k) + addStageTerms_(residual, k, + makeLocalContext(fvGeometry_, prevElemVariables_[k])); + addStageTerms_(residual, stageParams_->size()-1, + makeLocalContext(fvGeometry_, elemVariables_)); + } + + return residual; + } + } + +protected: + + //! add the terms of a stage to the current element residual + void addStageTerms_(ElementResidualVector& r, + std::size_t stageIdx, + const LocalContext& context) + { + LocalOperator localOperator(context); + if (!stageParams_->skipTemporal(stageIdx)) + residual.axpy(stageParams_->temporalWeight(stageIdx), + localOperator.evalStorage()); + if (!stageParams_->skipSpatial(stageIdx)) + residual.axpy(stageParams_->spatialWeight(stageIdx), + localOperator.evalFluxesAndSources()); + } + + /*! + * \brief Computes the derivatives with respect to the dofs of the given + * element and adds them to the global matrix. + * \return The element residual at the current solution. + */ + template + ElementResidualVector assembleJacobianAndResidual_(JacobianMatrix& A, + const PartialReassembler* partialReassembler = nullptr) + { + if constexpr (diffMethod == DiffMethod::numeric) + return assembleJacobianAndResidualNumeric_(A, partialReassembler); + else + DUNE_THROW(Dune::NotImplemented, "Analytic assembler for box"); + } + + /*! + * \brief Computes the derivatives by means of numeric differentiation + * and adds them to the global matrix. + * \return The element residual at the current solution. + * \note This specialization is for the box scheme with numeric differentiation + */ + template< class PartialReassembler = DefaultPartialReassembler > + ElementResidualVector assembleJacobianAndResidualNumeric_(JacobianMatrix& A, + const PartialReassembler* partialReassembler = nullptr) + { + // get the variables of the current stage + auto& curVariables = elemVariables(); + auto& curElemVolVars = curVariables.elemVolVars(); + const auto& x = curVariables.gridVariables().dofs(); + + const auto origResiduals = evalLocalResidual(); + const auto origElemSol = elementSolution(element(), x, fvGeometry().gridGeometry()); + auto elemSol = origElemSol; + + ////////////////////////////////////////////////////////////////////////////////////////////// + // Calculate derivatives of the residual of all dofs in element with respect to themselves. // + ////////////////////////////////////////////////////////////////////////////////////////////// + + ElementResidualVector partialDerivs(fvGeometry().numScv()); + for (const auto& scvI : scvs(fvGeometry())) + { + // dof index and corresponding actual pri vars + const auto globalI = scvI.dofIndex(); + const auto localI = scvI.localDofIndex(); + + const auto origCurVolVars = curElemVolVars[scvI]; + auto& curVolVars = curElemVolVars[scvI]; + + // calculate derivatives w.r.t to the privars at the dof at hand + for (int pvIdx = 0; pvIdx < numEq; pvIdx++) + { + partialDerivs = 0.0; + auto evalResiduals = [&](Scalar priVar) + { + // update the volume variables and compute element residual + elemSol[scvI.localDofIndex()][pvIdx] = priVar; + curVolVars.update(elemSol, problem(), element(), scvI); + return evalLocalResidual(); + }; + + // derive the residuals numerically + static const NumericEpsilon eps_{problem().paramGroup()}; + static const int numDiffMethod = getParamFromGroup(problem().paramGroup(), "Assembly.NumericDifferenceMethod"); + NumericDifferentiation::partialDerivative(evalResiduals, elemSol[localI][pvIdx], partialDerivs, + origResiduals, eps_(elemSol[localI][pvIdx], pvIdx), + numDiffMethod); + + // TODO: Distinguish between implicit/explicit here. For explicit schemes, + // no entries between different scvs of an element are reserved. Thus, + // we currently get an error when using explicit schemes. + // TODO: Doesn't this mean we only have to evaluate the residual for a single + // scv instead of calling evalLocalResidual()? That computes the residuals + // and derivs for all other scvs of the element, too, which are never used. + // Note: this is the same in the current implementation of master. + // Should we try to optimize this for explicit schemes? Or adjust the Jacobian pattern? + // update the global stiffness matrix with the current partial derivatives + for (const auto& scvJ : scvs(fvGeometry())) + { + const auto globalJ = scvJ.dofIndex(); + const auto localJ = scvJ.localDofIndex(); + + // don't add derivatives for green entities + if (!partialReassembler || partialReassembler->dofColor(globalJ) != EntityColor::green) + { + for (int eqIdx = 0; eqIdx < numEq; eqIdx++) + { + // A[i][col][eqIdx][pvIdx] is the rate of change of the + // the residual of equation 'eqIdx' at dof 'i' + // depending on the primary variable 'pvIdx' at dof 'col' + A[globalJ][globalI][eqIdx][pvIdx] += partialDerivs[localJ][eqIdx]; + } + } + } + + // restore the original element solution & volume variables + elemSol[localI][pvIdx] = origElemSol[localI][pvIdx]; + curVolVars = origCurVolVars; + + // TODO additional dof dependencies + } + } + + return origResiduals; + } + + //! Returns if a stationary problem is assembled + bool isStationary() const { return !stageParams_; } + + //! Return a reference to the underlying problem + //! TODO: Should grid vars return problem directly!? + const auto& problem_() const + { return elemVariables().gridVariables().gridVolVars().problem(); } + + DiffMethod diffMethod_; //!< the type of differentiation method + GridVariables& gridVariables_; //!< reference to the grid variables + FVElementGeometry fvGeometry_; //!< element-local finite volume geometry + ElementVariables elemVariables_; //!< element variables of the current stage + std::vector prevElemVariables_; //!< element variables of prior stages + + bool elementIsGhost_; + std::shared_ptr stageParams_; +}; + +} // end namespace Dumux::Experimental + +#endif -- GitLab From d1b61cd221711d9c26576ab62d41ae84ee661b1d Mon Sep 17 00:00:00 2001 From: "Dennis.Glaeser" Date: Thu, 22 Oct 2020 17:14:08 +0200 Subject: [PATCH 31/32] [assembly] add scheme-agnostic assembler --- dumux/assembly/CMakeLists.txt | 2 + dumux/assembly/assembler.hh | 418 +++++++++++++++++++++++++++++++ dumux/assembly/localassembler.hh | 56 +++++ 3 files changed, 476 insertions(+) create mode 100644 dumux/assembly/assembler.hh create mode 100644 dumux/assembly/localassembler.hh diff --git a/dumux/assembly/CMakeLists.txt b/dumux/assembly/CMakeLists.txt index 2543a6631f..63a13d1e43 100644 --- a/dumux/assembly/CMakeLists.txt +++ b/dumux/assembly/CMakeLists.txt @@ -1,4 +1,5 @@ install(FILES +assembler.hh boxlocalassembler.hh boxlocalresidual.hh cclocalassembler.hh @@ -10,6 +11,7 @@ fvlocalassemblerbase.hh fvlocalresidual.hh initialsolution.hh jacobianpattern.hh +localassembler.hh numericepsilon.hh partialreassembler.hh staggeredfvassembler.hh diff --git a/dumux/assembly/assembler.hh b/dumux/assembly/assembler.hh new file mode 100644 index 0000000000..4702c043ee --- /dev/null +++ b/dumux/assembly/assembler.hh @@ -0,0 +1,418 @@ +// -*- 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 Assembler class for residuals and jacobian matrices + * for grid-based numerical schemes. + */ +#ifndef DUMUX_ASSEMBLER_HH +#define DUMUX_ASSEMBLER_HH + +#include + +#include +#include +#include + +#include +#include +#include +#include + +#include "localassembler.hh" +#include "jacobianpattern.hh" + +namespace Dumux { + +//! Default types used for the linear system +template struct DefaultLinearSystemTraits; + +//! Default linear system types for Dune::BlockVector +template +struct DefaultLinearSystemTraits> +{ +private: + static constexpr int numEq = PrimaryVariables::size(); + using Scalar = typename PrimaryVariables::value_type; + using MatrixBlockType = Dune::FieldMatrix; + +public: + using ResidualVector = Dune::BlockVector; + using JacobianMatrix = Dune::BCRSMatrix; +}; + +/*! + * \ingroup Assembly + * \brief A linear system assembler (residual and Jacobian) for grid-based numerical schemes + * \tparam LO The local operator (evaluation of the terms of the equation) + * \tparam diffMethod The differentiation method to compute derivatives + * \tparam LST The linear system traits (types used for jacobians and residuals) + */ +template< class LO, DiffMethod diffMethod, + class LST = DefaultLinearSystemTraits > +class Assembler +{ + using ThisType = Assembler; + using GG = typename LO::GridVariables::GridGeometry; + using GridView = typename GG::GridView; + using Element = typename GridView::template Codim<0>::Entity; + +public: + //! export the types used for the linear system + using Scalar = typename LO::GridVariables::Scalar; + using JacobianMatrix = typename LST::JacobianMatrix; + using ResidualVector = typename LST::ResidualVector; + using ResidualType = ResidualVector; + + //! export the local operator type + using LocalOperator = LO; + + //! the local operator states the type of variables needed for evaluation + using GridVariables = typename LO::GridVariables; + + //! export the underlying grid geometry type + using GridGeometry = GG; + + //! export a grid-independent alias of the variables + using Variables = GridVariables; + + //! export type used for solution vectors + using SolutionVector = typename GridVariables::SolutionVector; + + /*! + * \brief The Constructor from a grid geometry. + * \param gridGeometry A grid geometry instance + * \note This assembler class is, after construction, defined for a specific equation + * (given by the template argument of the LocalOperator) and a specific grid + * geometry - which defines the connectivity of the degrees of freedom of the + * underlying discretization scheme on a particular grid. The evaluation point, + * consisting of a particular solution/variables/parameters may vary, and therefore, + * an instance of the grid variables is passed to the assembly functions. + */ + Assembler(std::shared_ptr gridGeometry) + : gridGeometry_(gridGeometry) + {} + + /*! + * \brief The Constructor from a grid geometry. + * \param gridGeometry A grid geometry instance + * \param isImplicit boolean to indicate whether an implicit or explicit pattern should be used + * \todo TODO replace bool with time stepping scheme once available + */ + Assembler(std::shared_ptr gridGeometry, + bool isImplicit) + : gridGeometry_(gridGeometry) + , isImplicit_(isImplicit) + {} + + /*! + * \brief Assembles the Jacobian matrix and the residual around the given evaluation point + * which is determined by the grid variables, containing all quantities required + * to evaluate the equations to be assembled. + * \param gridVariables The variables corresponding to the given solution state + * \note We assume the grid geometry on which the grid variables are defined + * to be the same as the one used to instantiate this class + */ + template + void assembleJacobianAndResidual(const GridVariables& gridVariables, + const PartialReassembler* partialReassembler = nullptr) + { + resetJacobian_(partialReassembler); + resetResidual_(); + + // TODO: Remove this assert? + assert(gridVariables.gridGeometry().numDofs() == gridGeometry().numDofs()); + + using LocalAssembler = Dumux::LocalAssembler; + assemble_([&](const Element& element) + { + auto fvGeometry = localView(gridGeometry()); + auto elemVars = localView(gridVariables); + + fvGeometry.bind(element); + elemVars.bind(element, fvGeometry); + + LocalAssembler localAssembler(element, fvGeometry, elemVars); + localAssembler.assembleJacobianAndResidual(*jacobian_, *residual_, partialReassembler); + }); + + // TODO: Put these into discretization-specific helpers? + enforceDirichletConstraints_(gridVariables, *jacobian_, *residual_); + enforceInternalConstraints_(gridVariables, *jacobian_, *residual_); + enforcePeriodicConstraints_(gridVariables, *jacobian_, *residual_); + } + + /*! + * \brief Assembles the Jacobian matrix of the discrete system of equations + * around a given state represented by the grid variables object. + */ + void assembleJacobian(const GridVariables& gridVariables) + { + resetJacobian_(); + + // TODO: Remove this assert? + assert(gridVariables.gridGeometry().numDofs() == gridGeometry().numDofs()); + + using LocalAssembler = Dumux::LocalAssembler; + assemble_([&](const Element& element) + { + auto fvGeometry = localView(gridGeometry()); + auto elemVars = localView(gridVariables); + + fvGeometry.bind(element); + elemVars.bind(element, fvGeometry); + + LocalAssembler localAssembler(element, fvGeometry, elemVars); + localAssembler.assembleJacobianAndResidual(*jacobian_); + }); + + // TODO: Put these into discretization-specific helpers? + enforceDirichletConstraints_(gridVariables, *jacobian_, *residual_); + enforceInternalConstraints_(gridVariables, *jacobian_, *residual_); + enforcePeriodicConstraints_(gridVariables, *jacobian_, *residual_); + } + + /*! + * \brief Assembles the residual for a given state represented by the provided + * grid variables object, using the internal residual vector to store the result. + */ + void assembleResidual(const GridVariables& gridVariables) + { + resetResidual_(); + assembleResidual(*residual_, gridVariables); + } + + /*! + * \brief Assembles the residual for a given state represented by the provided + * grid variables object, using the provided residual vector to store the result. + */ + void assembleResidual(ResidualVector& r, const GridVariables& gridVariables) const + { + using LocalAssembler = Dumux::LocalAssembler; + assemble_([&](const Element& element) + { + auto fvGeometry = localView(gridGeometry()); + auto elemVars = localView(gridVariables); + + fvGeometry.bind(element); + elemVars.bind(element, fvGeometry); + + LocalAssembler localAssembler(element, fvGeometry, elemVars); + localAssembler.assembleResidual(r); + }); + } + + //! Will become obsolete once the new linear solvers are available + Scalar residualNorm(const GridVariables& gridVars) const + { + ResidualType residual(numDofs()); + assembleResidual(residual, gridVars); + + // for box communicate the residual with the neighboring processes + if (GridGeometry::discMethod == DiscretizationMethod::box && gridView().comm().size() > 1) + { + using VertexMapper = typename GridGeometry::VertexMapper; + VectorCommDataHandleSum + sumResidualHandle(gridGeometry_->vertexMapper(), residual); + gridView().communicate(sumResidualHandle, + Dune::InteriorBorder_InteriorBorder_Interface, + Dune::ForwardCommunication); + } + + // calculate the square norm of the residual + Scalar result2 = residual.two_norm2(); + if (gridView().comm().size() > 1) + result2 = gridView().comm().sum(result2); + + using std::sqrt; + return sqrt(result2); + } + + + /*! + * \brief Tells the assembler which jacobian and residual to use. + * This also resizes the containers to the required sizes and sets the + * sparsity pattern of the jacobian matrix. + */ + void setLinearSystem(std::shared_ptr A, + std::shared_ptr r) + { + jacobian_ = A; + residual_ = r; + + // check and/or set the BCRS matrix's build mode + if (jacobian_->buildMode() == JacobianMatrix::BuildMode::unknown) + jacobian_->setBuildMode(JacobianMatrix::random); + 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(); + } + + /*! + * \brief The version without arguments uses the default constructor to create + * the jacobian and residual objects in this assembler if you don't need them outside this class + */ + void setLinearSystem() + { + jacobian_ = std::make_shared(); + jacobian_->setBuildMode(JacobianMatrix::random); + residual_ = std::make_shared(); + + setJacobianPattern(); + setResidualSize(); + } + + /*! + * \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 + // TODO: Does the bool need to be at compile time in getJacobianPattern? + const auto occupationPattern = isImplicit_ ? getJacobianPattern(gridGeometry()) + : getJacobianPattern(gridGeometry()); + + // export pattern to jacobian + occupationPattern.exportIdx(*jacobian_); + } + + //! Resizes the residual + void setResidualSize() + { residual_->resize(numDofs()); } + + //! Returns the number of degrees of freedom + std::size_t numDofs() const + { return gridGeometry().numDofs(); } + + //! The global finite volume geometry + const GridGeometry& gridGeometry() const + { return *gridGeometry_; } + + //! The gridview + const GridView& gridView() const + { return gridGeometry().gridView(); } + + //! The jacobian matrix + JacobianMatrix& jacobian() + { return *jacobian_; } + + //! The residual vector (rhs) + ResidualVector& residual() + { return *residual_; } + +protected: + // reset the residual vector to 0.0 + void resetResidual_() + { + if (!residual_) + { + residual_ = std::make_shared(); + setResidualSize(); + } + + (*residual_) = 0.0; + } + + // reset the Jacobian matrix to 0.0 + template + void resetJacobian_(const PartialReassembler *partialReassembler = nullptr) + { + if (!jacobian_) + { + jacobian_ = std::make_shared(); + jacobian_->setBuildMode(JacobianMatrix::random); + setJacobianPattern(); + } + + if (partialReassembler) + partialReassembler->resetJacobian(*this); + else + *jacobian_ = 0.0; + } + + /*! + * \brief A method assembling something per element + * \note Handles exceptions for parallel runs + * \throws NumericalProblem on all processes if something throwed during assembly + */ + template + void assemble_(AssembleElementFunc&& assembleElement) const + { + // a state that will be checked on all processes + bool succeeded = false; + + // 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 we get here, everything worked well on this process + succeeded = true; + } + // throw exception if a problem ocurred + catch (NumericalProblem &e) + { + std::cout << "rank " << gridView().comm().rank() + << " caught an exception while assembling:" << e.what() + << "\n"; + succeeded = false; + } + + // make sure everything worked well on all processes + if (gridView().comm().size() > 1) + succeeded = gridView().comm().min(succeeded); + + // if not succeeded rethrow the error on all processes + if (!succeeded) + DUNE_THROW(NumericalProblem, "A process did not succeed in linearizing the system"); + } + + void enforceDirichletConstraints_(const GridVariables& gridVars, JacobianMatrix& jac, SolutionVector& res) + { /*TODO: Implement / where to put this? Currently, local assemblers do this*/ } + + void enforceInternalConstraints_(const GridVariables& gridVars, JacobianMatrix& jac, SolutionVector& res) + { /*TODO: Implement / where to put this? Currently, local assemblers do this*/ } + + void enforcePeriodicConstraints_(const GridVariables& gridVars, JacobianMatrix& jac, SolutionVector& res) + { /*TODO: Implement / where to put this? Currently, local assemblers do this*/ } + +private: + //! the grid geometry on which it is assembled + std::shared_ptr gridGeometry_; + + //! states if an implicit of explicit scheme is used (affects jacobian pattern) + bool isImplicit_ = true; + + //! shared pointers to the jacobian matrix and residual + std::shared_ptr jacobian_; + std::shared_ptr residual_; +}; + +} // namespace Dumux + +#endif diff --git a/dumux/assembly/localassembler.hh b/dumux/assembly/localassembler.hh new file mode 100644 index 0000000000..b5fbe3da62 --- /dev/null +++ b/dumux/assembly/localassembler.hh @@ -0,0 +1,56 @@ +// -*- 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 Helper alias to select a local assembler based on the discretization scheme. + */ +#ifndef DUMUX_LOCAL_ASSEMBLER_HH +#define DUMUX_LOCAL_ASSEMBLER_HH + +#include +#include "fv/boxlocalassembler.hh" + +namespace Dumux { +namespace Impl { + + template + struct LocalAssemblerChooser; + + template + struct LocalAssemblerChooser + { using type = BoxLocalAssembler; }; + + template + using LocalAssemblerType = typename LocalAssemblerChooser::type; + +} // end namespace Immpl + +/*! + * \ingroup Assembly + * \brief Helper alias to select the local assembler type from an assembler. + */ +template +using LocalAssembler = Impl::LocalAssemblerType; + +} // namespace Dumux + +#endif -- GitLab From 31257106aa9721455759374daa77fea07599a848 Mon Sep 17 00:00:00 2001 From: "Dennis.Glaeser" Date: Thu, 22 Oct 2020 17:17:22 +0200 Subject: [PATCH 32/32] TEMP [test][1p][box] use new assembler --- .../1p/incompressible/main.cc | 33 +++++++++++-------- 1 file changed, 20 insertions(+), 13 deletions(-) diff --git a/test/porousmediumflow/1p/incompressible/main.cc b/test/porousmediumflow/1p/incompressible/main.cc index 3f8e191a6c..55434745fd 100644 --- a/test/porousmediumflow/1p/incompressible/main.cc +++ b/test/porousmediumflow/1p/incompressible/main.cc @@ -46,7 +46,9 @@ #include #include -#include +#include +#include +#include #include "properties.hh" #include "problem.hh" @@ -118,17 +120,14 @@ int main(int argc, char** argv) using Problem = GetPropType; auto problem = std::make_shared(gridGeometry); - // the solution vector - using SolutionVector = GetPropType; - SolutionVector x(gridGeometry->numDofs()); - // the grid variables using GridVariables = GetPropType; - auto gridVariables = std::make_shared(problem, gridGeometry); - gridVariables->init(x); + auto initX = [&] (auto& x) { x.resize(gridGeometry->numDofs()); x = 0.0; }; + auto gridVariables = std::make_shared(problem, gridGeometry, initX); // intialize the vtk output module - VtkOutputModule vtkWriter(*gridVariables, x, problem->name()); + using SolutionVector = GetPropType; + VtkOutputModule vtkWriter(*gridVariables, gridVariables->dofs(), problem->name()); using VelocityOutput = GetPropType; vtkWriter.addVelocityOutput(std::make_shared(*gridVariables)); using IOFields = GetPropType; @@ -136,15 +135,23 @@ int main(int argc, char** argv) vtkWriter.write(0.0); // create assembler & linear solver - using Assembler = FVAssembler; - auto assembler = std::make_shared(problem, gridGeometry, gridVariables); + // TODO: The operators or local operator could be a property, just as is LocalResidual currently + using ModelTraits = GetPropType; + using FluxVariables = GetPropType; + using ElemVariables = typename GridVariables::LocalView; + + using ImmiscibleOperators = FVImmiscibleOperators; + using LocalOperator = FVLocalOperator; + + using Assembler = Assembler; + auto assembler = std::make_shared(gridGeometry); using LinearSolver = SSORCGBackend; auto linearSolver = std::make_shared(); // solver the linear problem LinearPDESolver solver(assembler, linearSolver); - solver.solve(x); + solver.solve(*gridVariables); // output result to vtk vtkWriter.write(1.0); @@ -178,7 +185,7 @@ int main(int argc, char** argv) auto elemFluxVarsCache = localView(gridVariables->gridFluxVarsCache()); fvGeometry.bind(element); - elemVolVars.bind(element, fvGeometry, x); + elemVolVars.bind(element, fvGeometry, gridVariables->dofs()); elemFluxVarsCache.bind(element, fvGeometry, elemVolVars); velocityOutput.calculateVelocity(velocity, element, fvGeometry, elemVolVars, elemFluxVarsCache, 0); @@ -201,7 +208,7 @@ int main(int argc, char** argv) // For the mpfa test, write out the gradients in the scv centers if (getParam("IO.WriteMpfaVelocities", false)) - writeMpfaVelocities(*gridGeometry, *gridVariables, x); + writeMpfaVelocities(*gridGeometry, *gridVariables, gridVariables->dofs()); if (mpiHelper.rank() == 0) Parameters::print(); -- GitLab