Commit e881081a authored by Dennis Gläser's avatar Dennis Gläser
Browse files

Merge branch 'feature/newton-variables-backend' into 'master'

Feature/newton variables backend

See merge request !2291
parents 6e58bfb8 0a98d2b0
......@@ -51,4 +51,5 @@ timemanager.hh
valgrind.hh
variablelengthspline_.hh
variables.hh
variablesbackend.hh
DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/dumux/common)
......@@ -28,6 +28,7 @@
#include <utility>
#include <dune/common/hybridutilities.hh>
#include <dune/common/std/type_traits.hh>
#include <dumux/common/timeloop.hh>
......@@ -35,9 +36,25 @@
namespace Dune {
template <class FirstRow, class ... Args>
class MultiTypeBlockMatrix;
}
} // end namespace Dune
namespace Dumux {
namespace Detail {
template<class Assembler>
using AssemblerVariablesType = typename Assembler::Variables;
template<class Assembler>
inline constexpr bool exportsVariables = Dune::Std::is_detected_v<AssemblerVariablesType, Assembler>;
template<class A, bool exports = exportsVariables<A>> struct VariablesChooser;
template<class A> struct VariablesChooser<A, true> { using Type = AssemblerVariablesType<A>; };
template<class A> struct VariablesChooser<A, false> { using Type = typename A::ResidualType; };
template<class Assembler>
using AssemblerVariables = typename VariablesChooser<Assembler>::Type;
} // end namespace Detail
/*!
* \ingroup Common
......@@ -53,11 +70,19 @@ namespace Dumux {
template<class Assembler, class LinearSolver>
class PDESolver
{
using SolutionVector = typename Assembler::ResidualType;
using Scalar = typename Assembler::Scalar;
using TimeLoop = TimeLoopBase<Scalar>;
public:
//! export the type of variables that represent a numerical solution
using Variables = Detail::AssemblerVariables<Assembler>;
/*!
* \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> assembler,
std::shared_ptr<LinearSolver> 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
*/
......
// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
// vi: set et ts=4 sw=4 sts=4:
/****************************************************************************
* See the file COPYING for full copying permissions. *
* *
* This program is free software: you can redistribute it and/or modify *
* it under the terms of the GNU General Public License as published by *
* the Free Software Foundation, either version 3 of the License, or *
* (at your option) any later version. *
* *
* This program is distributed in the hope that it will be useful, *
* but WITHOUT ANY WARRANTY; without even the implied warranty of *
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
* GNU General Public License for more details. *
* *
* You should have received a copy of the GNU General Public License *
* along with this program. If not, see <http://www.gnu.org/licenses/>. *
*****************************************************************************/
/*!
* \file
* \ingroup Common
* \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_COMMON_VARIABLES_BACKEND_HH
#define DUMUX_COMMON_VARIABLES_BACKEND_HH
#include <array>
#include <utility>
#include <type_traits>
#include <dune/common/indices.hh>
#include <dune/common/typetraits.hh>
#include <dune/common/hybridutilities.hh>
#include <dune/common/std/type_traits.hh>
#include <dune/istl/bvector.hh>
// forward declaration
namespace Dune {
template<class... Args>
class MultiTypeBlockVector;
} // end namespace Dune
namespace Dumux {
/*!
* \ingroup Common
* \brief Class providing operations with primary variable vectors
*/
template<class DofVector, bool isScalar = Dune::IsNumber<DofVector>::value>
class DofBackend;
/*!
* \ingroup Common
* \brief Specialization providing operations for scalar/number types
*/
template<class Scalar>
class DofBackend<Scalar, true>
{
public:
using DofVector = Scalar; //!< the type of the dofs parametrizing the variables object
using SizeType = std::size_t;
//! Return the number of entries in the dof vector
static SizeType size(const DofVector& d)
{ return 1; }
//! Make a zero-initialized dof vector instance
static DofVector zeros(SizeType size)
{ return 0.0; }
//! Perform axpy operation (y += a * x)
static void axpy(Scalar a, const DofVector& x, DofVector& y)
{ y += a*x; }
};
/*!
* \ingroup Common
* \brief Specialization providing operations for block vectors
*/
template<class BT>
class DofBackend<Dune::BlockVector<BT>, false>
{
public:
using DofVector = Dune::BlockVector<BT>; //!< the type of the dofs parametrizing the variables object
using SizeType = std::size_t;
//! Return the number of entries in the dof vector
static SizeType size(const DofVector& d)
{ return d.size(); }
//! Make a zero-initialized dof vector instance
static DofVector zeros(SizeType size)
{ DofVector d; d.resize(size); return d; }
//! Perform axpy operation (y += a * x)
static void axpy(typename DofVector::field_type a, const DofVector& x, DofVector& y)
{ y.axpy(a, x); }
};
/*!
* \ingroup Common
* \brief Specialization providing operations for multitype block vectors
*/
template<class... Blocks>
class DofBackend<Dune::MultiTypeBlockVector<Blocks...>, false>
{
using DV = Dune::MultiTypeBlockVector<Blocks...>;
static constexpr auto numBlocks = DV::size();
using VectorSizeInfo = std::array<std::size_t, numBlocks>;
public:
using DofVector = DV; //!< the type of the dofs parametrizing the variables object
using SizeType = VectorSizeInfo;
//! Return the number of entries in the sub-dof-vectors
static SizeType size(const DofVector& d)
{
VectorSizeInfo result;
using namespace Dune::Hybrid;
forEach(std::make_index_sequence<numBlocks>{}, [&](auto i) {
result[i] = d[Dune::index_constant<i>{}].size();
});
return result;
}
//! Make a zero-initialized dof vector instance
static DofVector zeros(const SizeType& size)
{
DofVector result;
using namespace Dune::Hybrid;
forEach(std::make_index_sequence<numBlocks>{}, [&](auto i) {
result[Dune::index_constant<i>{}].resize(size[i]);
});
return result;
}
//! Perform axpy operation (y += a * x)
template<class Scalar, std::enable_if_t< Dune::IsNumber<Scalar>::value, int> = 0>
static void axpy(Scalar a, const DofVector& x, DofVector& y)
{ y.axpy(a, x); }
};
namespace Detail {
template<class Vars>
using SolutionVectorType = typename Vars::SolutionVector;
template<class Vars, bool varsExportSolution>
class VariablesBackend;
/*!
* \ingroup Common
* \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 Vars>
class VariablesBackend<Vars, false>
: public DofBackend<Vars>
{
using ParentType = DofBackend<Vars>;
public:
using Variables = Vars;
using typename ParentType::DofVector;
//! update to new solution vector
static void update(Variables& v, const DofVector& dofs)
{ v = dofs; }
//! return const reference to dof vector
static const DofVector& dofs(const Variables& v)
{ return v; }
//! return reference to dof vector
static DofVector& dofs(Variables& v)
{ return v; }
};
/*!
* \file
* \ingroup Common
* \brief Class providing operations for generic variable classes,
* containing primary and possibly also secondary variables.
*/
template<class Vars>
class VariablesBackend<Vars, true>
: public DofBackend<typename Vars::SolutionVector>
{
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); }
//! return const reference to dof vector
static const DofVector& dofs(const Variables& v)
{ return v.dofs(); }
//! return reference to dof vector
static DofVector& dofs(Variables& v)
{ return v.dofs(); }
};
} // end namespace Detail
/*!
* \ingroup Common
* \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<class Vars>
using VariablesBackend = Detail::VariablesBackend<Vars, Dune::Std::is_detected_v<Detail::SolutionVectorType, Vars>>;
} // end namespace Dumux
#endif
......@@ -50,7 +50,10 @@ template <class Assembler, class LinearSolver, class CouplingManager,
class MultiDomainNewtonSolver: public NewtonSolver<Assembler, LinearSolver, Reassembler, Comm>
{
using ParentType = NewtonSolver<Assembler, LinearSolver, Reassembler, Comm>;
using SolutionVector = typename Assembler::ResidualType;
using typename ParentType::Backend;
using typename ParentType::SolutionVector;
static constexpr bool assemblerExportsVariables = Detail::exportsVariables<Assembler>;
template<std::size_t i>
using PrimaryVariableSwitch = typename Detail::GetPVSwitchMultiDomain<Assembler, i>::type;
......@@ -63,6 +66,7 @@ class MultiDomainNewtonSolver: public NewtonSolver<Assembler, LinearSolver, Reas
using PriVarSwitchPtrTuple = typename Assembler::Traits::template Tuple<PrivarSwitchPtr>;
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::dofs(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<Assembler::Traits::numSubDomains>{}, [&](auto&& id)
{
this->initPriVarSwitch_(u, id, HasPriVarsSwitch<std::decay_t<decltype(id)>::value>{});
this->initPriVarSwitch_(vars, id, HasPriVarsSwitch<std::decay_t<decltype(id)>::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<Assembler::Traits::numSubDomains>{}, [&](auto&& id)
{
this->invokePriVarSwitch_(uCurrentIter[id], id, HasPriVarsSwitch<std::decay_t<decltype(id)>::value>{});
auto& uCurrentIter = Backend::dofs(varsCurrentIter)[id];
if constexpr (!assemblerExportsVariables)
this->invokePriVarSwitch_(this->assembler().gridVariables(id),
uCurrentIter, id, HasPriVarsSwitch<std::decay_t<decltype(id)>::value>{});
else
this->invokePriVarSwitch_(varsCurrentIter[id], uCurrentIter, id, HasPriVarsSwitch<std::decay_t<decltype(id)>::value>{});
});
ParentType::newtonEndStep(uCurrentIter, uLastIter);
couplingManager_->updateSolution(uCurrentIter);
ParentType::newtonEndStep(varsCurrentIter, uLastIter);
couplingManager_->updateSolution(Backend::dofs(varsCurrentIter));
}
private:
......@@ -150,60 +159,61 @@ private:
* \brief Reset the privar switch state, noop if there is no priVarSwitch
*/
template<std::size_t i>
void initPriVarSwitch_(SolutionVector&, Dune::index_constant<i> id, std::false_type) {}
void initPriVarSwitch_(Variables&, Dune::index_constant<i> id, std::false_type) {}
/*!
* \brief Switch primary variables if necessary
*/
template<std::size_t i>
void initPriVarSwitch_(SolutionVector& sol, Dune::index_constant<i> id, std::true_type)
void initPriVarSwitch_(Variables& vars, Dune::index_constant<i> id, std::true_type)
{
using namespace Dune::Hybrid;
auto& priVarSwitch = *elementAt(priVarSwitches_, id);
auto& sol = Backend::dofs(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<class SubSol, std::size_t i>
void invokePriVarSwitch_(SubSol&, Dune::index_constant<i> id, std::false_type) {}
template<class SubVars, class SubSol, std::size_t i>
void invokePriVarSwitch_(SubVars&, SubSol&, Dune::index_constant<i> id, std::false_type) {}
/*!
* \brief Switch primary variables if necessary
*/
template<class SubSol, std::size_t i>
void invokePriVarSwitch_(SubSol& uCurrentIter, Dune::index_constant<i> id, std::true_type)
template<class SubVars, class SubSol, std::size_t i>
void invokePriVarSwitch_(SubVars& subVars, SubSol& uCurrentIter, Dune::index_constant<i> 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);
}
}
}
......
......@@ -2,5 +2,6 @@ install(FILES
findscalarroot.hh
newtonconvergencewriter.hh
newtonsolver.hh
primaryvariableswitchadapter.hh
staggerednewtonconvergencewriter.hh
DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/dumux/nonlinear)
This diff is collapsed.
// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
// vi: set et ts=4 sw=4 sts=4:
/****************************************************************************
* See the file COPYING for full copying permissions. *
* *
* This program is free software: you can redistribute it and/or modify *
* it under the terms of the GNU General Public License as published by *
* the Free Software Foundation, either version 3 of the License, or *
* (at your option) any later version. *
* *
* This program is distributed in the hope that it will be useful, *
* but WITHOUT ANY WARRANTY; without even the implied warranty of *
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
* GNU General Public License for more details. *
* *
* You should have received a copy of the GNU General Public License *
* along with this program. If not, see <http://www.gnu.org/licenses/>. *
*****************************************************************************/
/*!
* \file
* \ingroup 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 <memory>
#include <dune/common/std/type_traits.hh>
namespace Dumux {
namespace Detail {
//! helper aliases to extract a primary variable switch from the VolumeVariables (if defined, yields int otherwise)
template<class Variables>
using DetectPVSwitch = typename Variables::VolumeVariables::PrimaryVariableSwitch;
template<class Variables>
using GetPVSwitch = Dune::Std::detected_or<int, DetectPVSwitch, Variables>;
template<class Variables>
using PrimaryVariableSwitch = typename GetPVSwitch<Variables>::type;
} // end namespace Detail
/*!
* \ingroup Nonlinear
* \brief Helper boolean to check if the given variables involve primary variable switching.
*/
template<class Variables>
inline constexpr bool hasPriVarsSwitch = typename Detail::GetPVSwitch<Variables>::value_t();
/*!
* \ingroup Nonlinear
* \brief An adapter for the Newton to manage models with primary variable switch
*/
template <class Variables, bool isValid = hasPriVarsSwitch<Variables>>
class PrimaryVariableSwitchAdapter
{
using PrimaryVariableSwitch = typename Detail::PrimaryVariableSwitch<Variables>;
public:
PrimaryVariableSwitchAdapter(const std::string& paramGroup = "")
{
const int priVarSwitchVerbosity = getParamFromGroup<int>(paramGroup, "PrimaryVariableSwitch.Verbosity", 1);
priVarSwitch_ = std::make_unique<PrimaryVariableSwitch>(priVarSwitchVerbosity);
}
/*!
* \brief Initialize the privar switch
*/
template<class SolutionVector>
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<class SolutionVector>
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<PrimaryVariableSwitch> priVarSwitch_;