Commit 60841d89 authored by Timo Koch's avatar Timo Koch
Browse files

Make Newtoncontroller independent of TypeTag

* [test][richards] Fix solver include
  (Used to compile before due to indirect include over newtoncontroller.
parent d049dfcc
......@@ -26,59 +26,64 @@
#ifndef DUMUX_NEWTON_CONTROLLER_HH
#define DUMUX_NEWTON_CONTROLLER_HH
#include <cmath>
#include <dune/common/exceptions.hh>
#include <dune/common/parallel/mpicollectivecommunication.hh>
#include <dune/common/parallel/mpihelper.hh>
#include <dune/istl/bvector.hh>
#include <dumux/common/properties.hh>
#include <dumux/common/exceptions.hh>
#include <dumux/common/math.hh>
#include <dumux/common/timeloop.hh>
#include <dumux/linear/seqsolverbackend.hh>
namespace Dumux {
/*!
* \ingroup Nonlinear
* \brief An implementation of a Newton controller
*
* \tparam Scalar the scalar type
* \tparam Comm the communication object used to communicate with all processes
* \note If you want to specialize only some methods but are happy with the
* defaults of the reference controller, derive your controller from
* this class and simply overload the required methods.
*/
template <class TypeTag>
template <class Scalar,
class Comm = Dune::CollectiveCommunication<Dune::MPIHelper::MPICommunicator> >
class NewtonController
{
using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
using Communicator = typename GridView::CollectiveCommunication;
using NumEqVector = typename GET_PROP_TYPE(TypeTag, NumEqVector);
static constexpr int numEq = GET_PROP_VALUE(TypeTag, NumEq);
public:
//! the communication type used to communicate with all processes
using Communication = Comm;
/*!
* \brief Constructor for stationary problems
*/
NewtonController(const Communicator& comm)
NewtonController(const Communication& comm = Dune::MPIHelper::getCollectiveCommunication(),
const std::string& paramGroup = "")
: comm_(comm)
, endIterMsgStream_(std::ostringstream::out)
, paramGroup_(paramGroup)
{
initParams_();
initParams_(paramGroup);
}
/*!
* \brief Constructor for instationary problems
*/
NewtonController(const Communicator& comm, std::shared_ptr<TimeLoop<Scalar>> timeLoop)
NewtonController(std::shared_ptr<TimeLoop<Scalar>> timeLoop,
const Communication& comm = Dune::MPIHelper::getCollectiveCommunication(),
const std::string& paramGroup = "")
: comm_(comm)
, timeLoop_(timeLoop)
, endIterMsgStream_(std::ostringstream::out)
, paramGroup_(paramGroup)
{
initParams_();
initParams_(paramGroup);
}
//! the communicator for parallel runs
const Communicator& communicator() const
const Communication& comm() const
{ return comm_; }
/*!
......@@ -251,8 +256,8 @@ public:
shift_ = max(shift_, shiftAtDof);
}
if (communicator().size() > 1)
shift_ = communicator().max(shift_);
if (comm().size() > 1)
shift_ = comm().max(shift_);
}
/*!
......@@ -289,34 +294,37 @@ public:
if (numSteps_ == 0)
{
Scalar norm2 = b.two_norm2();
if (communicator().size() > 1)
norm2 = communicator().sum(norm2);
if (comm().size() > 1)
norm2 = comm().sum(norm2);
using std::sqrt;
initialResidual_ = sqrt(norm2);
}
//! Copy into a standard block vector. This is necessary for all model _not_ using a FieldVector<Scalar, numEq> as
//! Copy into a standard block vector.
//! This is necessary for all model _not_ using a FieldVector<Scalar, blockSize> as
//! primary variables vector in combination with UMFPack or SuperLU as their interfaces are hard coded
//! to this field vector type in Dune ISTL
//! Could be avoided for vectors that already have the right type using SFINAE
//! but it shouldn't impact performance too much
Dune::BlockVector<NumEqVector> xTmp; xTmp.resize(b.size());
Dune::BlockVector<NumEqVector> bTmp(xTmp);
constexpr auto blockSize = JacobianMatrix::block_type::rows;
using BlockType = Dune::FieldVector<Scalar, blockSize>;
Dune::BlockVector<BlockType> xTmp; xTmp.resize(b.size());
Dune::BlockVector<BlockType> bTmp(xTmp);
for (unsigned int i = 0; i < b.size(); ++i)
for (unsigned int j = 0; j < numEq; ++j)
for (unsigned int j = 0; j < blockSize; ++j)
bTmp[i][j] = b[i][j];
int converged = ls.solve(A, xTmp, bTmp);
for (unsigned int i = 0; i < x.size(); ++i)
for (unsigned int j = 0; j < numEq; ++j)
for (unsigned int j = 0; j < blockSize; ++j)
x[i][j] = xTmp[i][j];
// make sure all processes converged
int convergedRemote = converged;
if (communicator().size() > 1)
convergedRemote = communicator().min(converged);
if (comm().size() > 1)
convergedRemote = comm().min(converged);
if (!converged) {
DUNE_THROW(NumericalProblem,
......@@ -327,24 +335,11 @@ public:
"Linear solver did not converge on a remote process");
}
}
catch (Dune::MatrixBlockError e) {
// make sure all processes converged
int converged = 0;
if (communicator().size() > 1)
converged = communicator().min(converged);
NumericalProblem p;
std::string msg;
std::ostringstream ms(msg);
ms << e.what() << "M=" << A[e.r][e.c];
p.message(ms.str());
throw p;
}
catch (const Dune::Exception &e) {
// make sure all processes converged
int converged = 0;
if (communicator().size() > 1)
converged = communicator().min(converged);
if (comm().size() > 1)
converged = comm().min(converged);
NumericalProblem p;
p.message(e.what());
......@@ -515,15 +510,19 @@ public:
* \brief Returns true if the Newton method ought to be chatty.
*/
bool verbose() const
{ return verbose_ && communicator().rank() == 0; }
{ return verbose_ && comm().rank() == 0; }
/*!
* \brief Returns the parameter group
*/
const std::string& paramGroup() const
{ return paramGroup_; }
protected:
//! initialize the parameters by reading from the parameter tree
void initParams_()
void initParams_(const std::string& group = "")
{
const std::string group = GET_PROP_VALUE(TypeTag, ModelParameterGroup);
useLineSearch_ = getParamFromGroup<bool>(group, "Newton.UseLineSearch");
enableAbsoluteResidualCriterion_ = getParamFromGroup<bool>(group, "Newton.EnableAbsoluteResidualCriterion");
enableShiftCriterion_ = getParamFromGroup<bool>(group, "Newton.EnableShiftCriterion");
......@@ -590,8 +589,6 @@ protected:
using std::abs;
using std::max;
// iterate over all primary variables
// note: we use PrimaryVariables::dimension (== numEq)
// for compatibility with the staggered grid implementation
for (int j = 0; j < PrimaryVariables::dimension; ++j) {
Scalar eqErr = abs(priVars1[j] - priVars2[j]);
eqErr /= max<Scalar>(1.0,abs(priVars1[j] + priVars2[j])/2);
......@@ -601,8 +598,8 @@ protected:
return result;
}
//! The grid view's communicator
const Communicator& comm_;
//! The communication object
Communication comm_;
//! The time loop for stationary simulations
std::shared_ptr<TimeLoop<Scalar>> timeLoop_;
......@@ -640,6 +637,9 @@ protected:
bool enableShiftCriterion_;
bool enableResidualCriterion_;
bool satisfyResidualAndShiftCriterion_;
//! the parameter group for getting parameters from the parameter tree
std::string paramGroup_;
};
} // end namespace Dumux
......
......@@ -25,9 +25,12 @@
#ifndef DUMUX_STAGGERED_NEWTON_CONTROLLER_HH
#define DUMUX_STAGGERED_NEWTON_CONTROLLER_HH
#include <dumux/common/properties.hh>
#include <dumux/common/exceptions.hh>
#include <dune/common/indices.hh>
#include <dune/common/hybridutilities.hh>
#include <dune/common/fvector.hh>
#include <dune/istl/bvector.hh>
#include <dumux/common/exceptions.hh>
#include <dumux/nonlinear/newtoncontroller.hh>
#include <dumux/linear/linearsolveracceptsmultitypematrix.hh>
#include <dumux/linear/matrixconverter.hh>
......@@ -40,26 +43,11 @@ namespace Dumux {
* \brief A newton controller for staggered finite volume schemes
*/
template <class TypeTag>
class StaggeredNewtonController : public NewtonController<TypeTag>
template <class Scalar,
class Comm = Dune::CollectiveCommunication<Dune::MPIHelper::MPICommunicator> >
class StaggeredNewtonController : public NewtonController<Scalar, Comm>
{
using ParentType = NewtonController<TypeTag>;
using Problem = typename GET_PROP_TYPE(TypeTag, Problem);
using SolutionVector = typename GET_PROP_TYPE(TypeTag, SolutionVector);
using JacobianMatrix = typename GET_PROP_TYPE(TypeTag, JacobianMatrix);
using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
using Communicator = typename GridView::CollectiveCommunication;
using DofTypeIndices = typename GET_PROP(TypeTag, DofTypeIndices);
typename DofTypeIndices::CellCenterIdx cellCenterIdx;
typename DofTypeIndices::FaceIdx faceIdx;
enum {
numEqCellCenter = GET_PROP_VALUE(TypeTag, NumEqCellCenter),
numEqFace = GET_PROP_VALUE(TypeTag, NumEqFace)
};
using ParentType = NewtonController<Scalar, Comm>;
public:
using ParentType::ParentType;
......@@ -90,8 +78,9 @@ public:
this->initialResidual_ = b.two_norm();
// check matrix sizes
assert(A[cellCenterIdx][cellCenterIdx].N() == A[cellCenterIdx][faceIdx].N());
assert(A[faceIdx][cellCenterIdx].N() == A[faceIdx][faceIdx].N());
using namespace Dune::Indices;
assert(A[_0][_0].N() == A[_0][_1].N());
assert(A[_1][_0].N() == A[_1][_1].N());
// create the bcrs matrix the IterativeSolver backend can handle
const auto M = MatrixConverter<JacobianMatrix>::multiTypeToBCRSMatrix(A);
......@@ -158,14 +147,15 @@ public:
this->lineSearchUpdate_(assembler, uCurrentIter, uLastIter, deltaU);
}
else {
for (unsigned int i = 0; i < uLastIter[cellCenterIdx].size(); ++i) {
uCurrentIter[cellCenterIdx][i] = uLastIter[cellCenterIdx][i];
uCurrentIter[cellCenterIdx][i] -= deltaU[cellCenterIdx][i];
}
for (unsigned int i = 0; i < uLastIter[faceIdx].size(); ++i) {
uCurrentIter[faceIdx][i] = uLastIter[faceIdx][i];
uCurrentIter[faceIdx][i] -= deltaU[faceIdx][i];
}
using namespace Dune::Hybrid;
forEach(integralRange(Dune::Hybrid::size(uLastIter)), [&](const auto dofTypeIdx)
{
for (unsigned int i = 0; i < uLastIter[dofTypeIdx].size(); ++i)
{
uCurrentIter[dofTypeIdx][i] = uLastIter[dofTypeIdx][i];
uCurrentIter[dofTypeIdx][i] -= deltaU[dofTypeIdx][i];
}
});
if (this->enableResidualCriterion_)
{
......@@ -186,30 +176,27 @@ public:
* \param uLastIter The current iterative solution
* \param deltaU The difference between the current and the next solution
*/
template<class SolutionVector>
void newtonUpdateShift(const SolutionVector &uLastIter,
const SolutionVector &deltaU)
{
this->shift_ = 0;
for (int i = 0; i < int(uLastIter[cellCenterIdx].size()); ++i) {
auto uNewI = uLastIter[cellCenterIdx][i];
uNewI -= deltaU[cellCenterIdx][i];
Scalar shiftAtDof = this->relativeShiftAtDof_(uLastIter[cellCenterIdx][i],
uNewI);
this->shift_ = std::max(this->shift_, shiftAtDof);
}
for (int i = 0; i < int(uLastIter[faceIdx].size()); ++i) {
auto uNewI = uLastIter[faceIdx][i];
uNewI -= deltaU[faceIdx][i];
using namespace Dune::Hybrid;
forEach(integralRange(Dune::Hybrid::size(uLastIter)), [&](const auto dofTypeIdx)
{
for (int i = 0; i < int(uLastIter[dofTypeIdx].size()); ++i)
{
auto uNewI = uLastIter[dofTypeIdx][i];
uNewI -= deltaU[dofTypeIdx][i];
Scalar shiftAtDof = this->relativeShiftAtDof_(uLastIter[faceIdx][i],
uNewI);
this->shift_ = std::max(this->shift_, shiftAtDof);
}
Scalar shiftAtDof = this->relativeShiftAtDof_(uLastIter[dofTypeIdx][i], uNewI);
this->shift_ = std::max(this->shift_, shiftAtDof);
}
});
if (this->communicator().size() > 1)
this->shift_ = this->communicator().max(this->shift_);
if (this->comm().size() > 1)
this->shift_ = this->comm().max(this->shift_);
}
};
......
......@@ -40,39 +40,22 @@ namespace Dumux
/*!
* \ingroup PorousmediumCompositional
* \brief A newton controller that handles primary variable switches
* \todo make this independent of TypeTag by making PrimaryVariableSwitch a template argument
* and extracting everything model specific from there
* \todo Implement for volume variable caching enabled
*/
template <class TypeTag>
class PriVarSwitchNewtonController : public NewtonController<TypeTag>
class PriVarSwitchNewtonController : public NewtonController<typename GET_PROP_TYPE(TypeTag, Scalar)>
{
using ParentType = NewtonController<TypeTag>;
using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
using Communicator = typename GridView::CollectiveCommunication;
using NumEqVector = typename GET_PROP_TYPE(TypeTag, NumEqVector);
using ParentType = NewtonController<Scalar>;
using PrimaryVariableSwitch = typename GET_PROP_TYPE(TypeTag, PrimaryVariableSwitch);
using ElementSolution = typename GET_PROP_TYPE(TypeTag, ElementSolutionVector);
static constexpr int numEq = GET_PROP_VALUE(TypeTag, NumEq);
static constexpr bool isBox = GET_PROP_VALUE(TypeTag, DiscretizationMethod) == DiscretizationMethods::Box;
// using ElementSolution = typename GET_PROP_TYPE(TypeTag, ElementSolutionVector);
// static constexpr bool isBox = GET_PROP_VALUE(TypeTag, DiscretizationMethod) == DiscretizationMethods::Box;
public:
/*!
* \ingroup Newton
* \brief Constructor for stationary problems
*/
PriVarSwitchNewtonController(const Communicator& comm)
: ParentType(comm)
, switchedInLastIteration_(false)
{}
/*!
* \brief Constructor for stationary problems
*/
PriVarSwitchNewtonController(const Communicator& comm, std::shared_ptr<TimeLoop<Scalar>> timeLoop)
: ParentType(comm, timeLoop)
, switchedInLastIteration_(false)
{}
using ParentType::ParentType;
/*!
* \brief Returns true if the error of the solution is below the
......@@ -201,7 +184,7 @@ private:
//! the class handling the primary variable switch
std::unique_ptr<PrimaryVariableSwitch> priVarSwitch_;
//! if we switched primary variables in the last iteration
bool switchedInLastIteration_;
bool switchedInLastIteration_ = false;
};
} // end namespace Dumux
......
......@@ -35,13 +35,11 @@ namespace Dumux {
* \brief A nonequilibrium specific controller for the newton solver.
* This controller calls the velocity averaging in the problem after each iteration.
*/
template <class TypeTag>
class NonEquilibriumNewtonController : public NewtonController<TypeTag>
template <class Scalar,
class Comm = Dune::CollectiveCommunication<Dune::MPIHelper::MPICommunicator> >
class NonEquilibriumNewtonController : public NewtonController<Scalar, Comm>
{
using ParentType = NewtonController<TypeTag>;
using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
using Communicator = typename GridView::CollectiveCommunication;
using ParentType = NewtonController<Scalar, Comm>;
public:
using ParentType::ParentType;
......
......@@ -34,38 +34,23 @@ namespace Dumux {
*
* This controller 'knows' what a 'physically meaningful' solution is
* and can thus do update smarter than the plain Newton controller.
*
* \todo make this typetag independent by extracting anything model specific from assembler
* or from possible ModelTraits.
*/
template <class TypeTag>
class RichardsNewtonController : public NewtonController<TypeTag>
class RichardsNewtonController : public NewtonController<typename GET_PROP_TYPE(TypeTag, Scalar)>
{
using ParentType = NewtonController<TypeTag>;
using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
using SolutionVector = typename GET_PROP_TYPE(TypeTag, SolutionVector);
using SpatialParams = typename GET_PROP_TYPE(TypeTag, SpatialParams);
using FVElementGeometry = typename GET_PROP_TYPE(TypeTag, FVElementGeometry);
using ParentType = NewtonController<Scalar>;
using MaterialLaw = typename GET_PROP_TYPE(TypeTag, MaterialLaw);
using Problem = typename GET_PROP_TYPE(TypeTag, Problem);
using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
using Communicator = typename GridView::CollectiveCommunication;
using ElementSolution = typename GET_PROP_TYPE(TypeTag, ElementSolutionVector);
using Indices = typename GET_PROP_TYPE(TypeTag, Indices);
enum { pressureIdx = Indices::pressureIdx };
public:
/*!
* \brief Constructor for stationary problems
*/
RichardsNewtonController(const Communicator& comm)
: ParentType(comm)
{}
/*!
* \brief Constructor for stationary problems
*/
RichardsNewtonController(const Communicator& comm, std::shared_ptr<TimeLoop<Scalar>> timeLoop)
: ParentType(comm, timeLoop)
{}
using ParentType::ParentType;
/*!
* \brief Update the current solution of the newton method
......@@ -73,7 +58,7 @@ public:
* This is basically the step
* \f[ u^{k+1} = u^k - \Delta u^k \f]
*
* \param assembler TODO docme!
* \param assembler The Jacobian assembler
* \param uCurrentIter The solution 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
......@@ -86,8 +71,7 @@ public:
const SolutionVector &deltaU)
{
ParentType::newtonUpdate(assembler, uCurrentIter, uLastIter, deltaU);
const std::string group = GET_PROP_VALUE(TypeTag, ModelParameterGroup);
if (!this->useLineSearch_ && getParamFromGroup<bool>(group, "Newton.EnableChop"))
if (!this->useLineSearch_ && getParamFromGroup<bool>(this->paramGroup(), "Newton.EnableChop"))
{
// do not clamp anything after 5 iterations
if (this->numSteps_ > 4)
......
......@@ -176,9 +176,9 @@ int main(int argc, char** argv) try
auto linearSolver = std::make_shared<LinearSolver>();
// the non-linear solver
using NewtonController = StaggeredNewtonController<TypeTag>;
using NewtonController = StaggeredNewtonController<Scalar>;
using NewtonMethod = Dumux::NewtonMethod<NewtonController, Assembler, LinearSolver>;
auto newtonController = std::make_shared<NewtonController>(leafGridView.comm(), timeLoop);
auto newtonController = std::make_shared<NewtonController>(timeLoop);
NewtonMethod nonLinearSolver(newtonController, assembler, linearSolver);
// time loop
......
......@@ -174,9 +174,9 @@ int main(int argc, char** argv) try
auto linearSolver = std::make_shared<LinearSolver>();
// the non-linear solver
using NewtonController = StaggeredNewtonController<TypeTag>;
using NewtonController = StaggeredNewtonController<Scalar>;
using NewtonMethod = Dumux::NewtonMethod<NewtonController, Assembler, LinearSolver>;
auto newtonController = std::make_shared<NewtonController>(leafGridView.comm(), timeLoop);
auto newtonController = std::make_shared<NewtonController>(timeLoop);
NewtonMethod nonLinearSolver(newtonController, assembler, linearSolver);
// set up two planes over which fluxes are calculated
......
......@@ -171,9 +171,9 @@ int main(int argc, char** argv) try
auto linearSolver = std::make_shared<LinearSolver>();
// the non-linear solver
using NewtonController = StaggeredNewtonController<TypeTag>;
using NewtonController = StaggeredNewtonController<Scalar>;
using NewtonMethod = Dumux::NewtonMethod<NewtonController, Assembler, LinearSolver>;
auto newtonController = std::make_shared<NewtonController>(leafGridView.comm(), timeLoop);
auto newtonController = std::make_shared<NewtonController>(timeLoop);
NewtonMethod nonLinearSolver(newtonController, assembler, linearSolver);
// time loop
......
......@@ -158,9 +158,10 @@ int main(int argc, char** argv) try
auto linearSolver = std::make_shared<LinearSolver>();
// the non-linear solver
using NewtonController = StaggeredNewtonController<TypeTag>;
using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
using NewtonController = StaggeredNewtonController<Scalar>;
using NewtonMethod = Dumux::NewtonMethod<NewtonController, Assembler, LinearSolver>;
auto newtonController = std::make_shared<NewtonController>(leafGridView.comm());
auto newtonController = std::make_shared<NewtonController>();
NewtonMethod nonLinearSolver(newtonController, assembler, linearSolver);
// linearize & solve
......
......@@ -157,9 +157,10 @@ int main(int argc, char** argv) try
auto linearSolver = std::make_shared<LinearSolver>();
// the non-linear solver
using NewtonController = StaggeredNewtonController<TypeTag>;
using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
using NewtonController = StaggeredNewtonController<Scalar>;
using NewtonMethod = Dumux::NewtonMethod<NewtonController, Assembler, LinearSolver>;
auto newtonController = std::make_shared<NewtonController>(leafGridView.comm());
auto newtonController = std::make_shared<NewtonController>();
NewtonMethod nonLinearSolver(newtonController, assembler, linearSolver);
// linearize & solve
......
......@@ -173,9 +173,9 @@ int main(int argc, char** argv) try
auto linearSolver = std::make_shared<LinearSolver>();
// the non-linear solver
using NewtonController = StaggeredNewtonController<TypeTag>;
using NewtonController = StaggeredNewtonController<Scalar>;
using NewtonMethod = Dumux::NewtonMethod<NewtonController, Assembler, LinearSolver>;
auto newtonController = std::make_shared<NewtonController>(leafGridView.comm(), timeLoop);
auto newtonController = std::make_shared<NewtonController>(timeLoop);
NewtonMethod nonLinearSolver(newtonController, assembler, linearSolver);
// time loop
......
......@@ -172,9 +172,9 @@ int main(int argc, char** argv) try
auto linearSolver = std::make_shared<LinearSolver>();
// the non-linear solver
using NewtonController = StaggeredNewtonController<TypeTag>;
using NewtonController = StaggeredNewtonController<Scalar>;
using NewtonMethod = Dumux::NewtonMethod<NewtonController, Assembler, LinearSolver>;
auto newtonController = std::make_shared<NewtonController>(leafGridView.comm(), timeLoop);
auto newtonController = std::make_shared<NewtonController>(timeLoop);
NewtonMethod nonLinearSolver(newtonController, assembler, linearSolver);
// time loop
......
......@@ -171,9 +171,9 @@ int main(int argc, char** argv) try
auto linearSolver = std::make_shared<LinearSolver>();
// the non-linear solver
using NewtonController = StaggeredNewtonController<TypeTag>;
using NewtonController = StaggeredNewtonController<Scalar>;
using NewtonMethod = Dumux::NewtonMethod<NewtonController, Assembler, LinearSolver>;
auto newtonController = std::make_shared<NewtonController>(leafGridView