Skip to content
Snippets Groups Projects
Commit d99d5076 authored by Kilian Weishaupt's avatar Kilian Weishaupt
Browse files

Delete obsolete NewtonController and NewtonMethod headers

parent 76995aed
No related branches found
No related tags found
1 merge request!898Use NewtonSolver
Expand all Hide whitespace changes
Inline Side-by-side
Showing
with 4 additions and 1586 deletions
dumux/nonlinear/CMakeLists.txt
+ 1
− 2
View file @ d99d5076
#install headers #install headers
install(FILES install(FILES
newtoncontroller.hh
newtonconvergencewriter.hh newtonconvergencewriter.hh
newtonmethod.hh newtonsolver.hh
DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/dumux/nonlinear) DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/dumux/nonlinear)
dumux/nonlinear/newtoncontroller.hh deleted 100644 → 0
+ 0
− 651
View file @ 76995aed
This diff is collapsed.
dumux/nonlinear/newtonmethod.hh deleted 100644 → 0
+ 0
− 222
View file @ 76995aed
// -*- 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 2 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 The algorithmic part of the multi dimensional newton method.
*
* In order to use the method you need a Newtoncontroller
*/
#ifndef DUMUX_NEWTONMETHOD_HH
#define DUMUX_NEWTONMETHOD_HH
#include <memory>
#include <iostream>
#include <dune/common/timer.hh>
#include <dune/istl/istlexception.hh>
#include <dumux/common/exceptions.hh>
#include <dumux/common/properties.hh>
#include <dumux/common/parameters.hh>
#include <dune/common/deprecated.hh>
#warning "This file is deprecated. Use NewtonSolver instead."
namespace Dumux {
/*!
* \ingroup Newton
* \ingroup Nonlinear
* \brief The algorithmic part of the multi dimensional newton method.
*
* \tparam NewtonController the controller class is the driver implementation
* \tparam JacobianAssembler an assembler for the Jacobian and the residual
* \tparam LinearSolver the linear solver used to solve one iteration
*/
template <class NewtonController, class JacobianAssembler, class LinearSolver>
class DUNE_DEPRECATED_MSG("Use NewtonSolver instead.") NewtonMethod
{
//! provide an interface as a form of type erasure
//! this is the minimal requirements a convergence write passed to a newton method has to fulfill
struct ConvergenceWriterInferface
{
template<class SolutionVector>
void write(const SolutionVector &uLastIter, const SolutionVector &deltaU, const SolutionVector &residual) {}
};
public:
NewtonMethod(std::shared_ptr<NewtonController> controller,
std::shared_ptr<JacobianAssembler> assembler,
std::shared_ptr<LinearSolver> linearSolver,
const std::string& modelParamGroup = "")
: controller_(controller)
, assembler_(assembler)
, linearSolver_(linearSolver)
{
// set the linear system (matrix & residual) in the assembler
assembler_->setLinearSystem();
// set a different default for the linear solver residual reduction
// within the Newton the linear solver doesn't need to solve too exact
using Scalar = typename LinearSolver::Scalar;
linearSolver_->setResidualReduction(getParamFromGroup<Scalar>(modelParamGroup, "LinearSolver.ResidualReduction", 1e-6));
}
/*!
* \brief Run the newton method to solve a non-linear system.
* The controller is responsible for all the strategic decisions.
*/
template<class SolutionVector, class ConvergenceWriter = ConvergenceWriterInferface>
bool solve(SolutionVector& uCurrentIter, const std::unique_ptr<ConvergenceWriter>& convWriter = nullptr)
{
try
{
// the given solution is the initial guess
SolutionVector uLastIter(uCurrentIter);
SolutionVector deltaU(uCurrentIter);
Dune::Timer assembleTimer(false);
Dune::Timer solveTimer(false);
Dune::Timer updateTimer(false);
// tell the controller that we begin solving
controller_->newtonBegin(uCurrentIter);
// execute the method as long as the controller thinks
// that we should do another iteration
while (controller_->newtonProceed(uCurrentIter, controller_->newtonConverged()))
{
// notify the controller that we're about to start
// a new timestep
controller_->newtonBeginStep(uCurrentIter);
// make the current solution to the old one
if (controller_->newtonNumSteps() > 0)
uLastIter = uCurrentIter;
if (controller_->verbose()) {
std::cout << "Assemble: r(x^k) = dS/dt + div F - q; M = grad r";
std::cout.flush();
}
///////////////
// assemble
///////////////
// linearize the problem at the current solution
assembleTimer.start();
controller_->assembleLinearSystem(*assembler_, uCurrentIter);
assembleTimer.stop();
///////////////
// linear solve
///////////////
// Clear the current line using an ansi escape
// sequence. for an explanation see
// http://en.wikipedia.org/wiki/ANSI_escape_code
const char clearRemainingLine[] = { 0x1b, '[', 'K', 0 };
if (controller_->verbose()) {
std::cout << "\rSolve: M deltax^k = r";
std::cout << clearRemainingLine;
std::cout.flush();
}
// solve the resulting linear equation system
solveTimer.start();
// set the delta vector to zero before solving the linear system!
deltaU = 0;
// ask the controller to solve the linearized system
controller_->solveLinearSystem(*linearSolver_,
assembler_->jacobian(),
deltaU,
assembler_->residual());
solveTimer.stop();
///////////////
// update
///////////////
if (controller_->verbose()) {
std::cout << "\rUpdate: x^(k+1) = x^k - deltax^k";
std::cout << clearRemainingLine;
std::cout.flush();
}
updateTimer.start();
// update the current solution (i.e. uOld) with the delta
// (i.e. u). The result is stored in u
controller_->newtonUpdate(*assembler_, uCurrentIter, uLastIter, deltaU);
updateTimer.stop();
// tell the controller that we're done with this iteration
controller_->newtonEndStep(*assembler_, uCurrentIter, uLastIter);
// if a convergence writer was specified compute residual and write output
if (convWriter)
{
assembler_->assembleResidual(uCurrentIter);
convWriter->write(uLastIter, deltaU, assembler_->residual());
}
}
// tell controller we are done
controller_->newtonEnd();
// reset state if newton failed
if (!controller_->newtonConverged())
{
controller_->newtonFail(*assembler_, uCurrentIter);
return false;
}
// tell controller we converged successfully
controller_->newtonSucceed();
if (controller_->verbose()) {
const auto elapsedTot = assembleTimer.elapsed() + solveTimer.elapsed() + updateTimer.elapsed();
std::cout << "Assemble/solve/update time: "
<< assembleTimer.elapsed() << "(" << 100*assembleTimer.elapsed()/elapsedTot << "%)/"
<< solveTimer.elapsed() << "(" << 100*solveTimer.elapsed()/elapsedTot << "%)/"
<< updateTimer.elapsed() << "(" << 100*updateTimer.elapsed()/elapsedTot << "%)"
<< "\n";
}
return true;
}
catch (const NumericalProblem &e)
{
if (controller_->verbose())
std::cout << "Newton: Caught exception: \"" << e.what() << "\"\n";
controller_->newtonFail(*assembler_, uCurrentIter);
return false;
}
}
private:
std::shared_ptr<NewtonController> controller_;
std::shared_ptr<JacobianAssembler> assembler_;
std::shared_ptr<LinearSolver> linearSolver_;
};
} // end namespace Dumux
#endif
dumux/nonlinear/staggerednewtoncontroller.hh deleted 100644 → 0
+ 0
− 299
View file @ 76995aed
// -*- 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 2 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
* \ingroup StaggeredDiscretization
* \brief A newton controller for staggered finite volume schemes
*/
#ifndef DUMUX_STAGGERED_NEWTON_CONTROLLER_HH
#define DUMUX_STAGGERED_NEWTON_CONTROLLER_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>
#include <dune/common/deprecated.hh>
#warning "This file is deprecated. Use NewtonSolver instead."
namespace Dumux {
/*!
* \ingroup Nonlinear
* \ingroup StaggeredDiscretization
* \brief A newton controller for staggered finite volume schemes
*/
template <class Scalar,
class Comm = Dune::CollectiveCommunication<Dune::MPIHelper::MPICommunicator> >
class DUNE_DEPRECATED_MSG("Use NewtonSolver instead.")
StaggeredNewtonController : public NewtonController<Scalar, Comm>
{
using ParentType = NewtonController<Scalar, Comm>;
public:
using ParentType::ParentType;
/*!
* \brief Solve the linear system of equations \f$\mathbf{A}x - b = 0\f$.
*
* Throws Dumux::NumericalProblem if the linear solver didn't
* converge.
*
* If the linear solver doesn't accept multitype matrices we copy the matrix
* into a 1x1 block BCRS matrix for solving.
*
* \param ls the linear solver
* \param A The matrix of the linear system of equations
* \param x The vector which solves the linear system
* \param b The right hand side of the linear system
*/
template<class LinearSolver, class JacobianMatrix, class SolutionVector>
void solveLinearSystem(LinearSolver& ls,
JacobianMatrix& A,
SolutionVector& x,
SolutionVector& b)
{
// check matrix sizes
assert(checkMatrix_(A) && "Sub blocks of MultiType matrix have wrong sizes!");
try
{
if (this->numSteps_ == 0)
{
Scalar norm2 = b.two_norm2();
if (this->comm().size() > 1)
norm2 = this->comm().sum(norm2);
using std::sqrt;
this->initialResidual_ = sqrt(norm2);
}
// solve by calling the appropriate implementation depending on whether the linear solver
// is capable of handling MultiType matrices or not
bool converged = solveLinearSystem_(ls, A, x, b,
std::integral_constant<bool, linearSolverAcceptsMultiTypeMatrix<LinearSolver>::value>());
// make sure all processes converged
int convergedRemote = converged;
if (this->comm().size() > 1)
convergedRemote = this->comm().min(converged);
if (!converged) {
DUNE_THROW(NumericalProblem,
"Linear solver did not converge");
}
else if (!convergedRemote) {
DUNE_THROW(NumericalProblem,
"Linear solver did not converge on a remote process");
}
}
catch (const Dune::Exception &e) {
// make sure all processes converged
int converged = 0;
if (this->comm().size() > 1)
converged = this->comm().min(converged);
NumericalProblem p;
p.message(e.what());
throw p;
}
}
/*!
* \brief Update the current solution with a delta vector.
*
* The error estimates required for the newtonConverged() and
* newtonProceed() methods should be updated inside this method.
*
* Different update strategies, such as line search and chopped
* updates can be implemented. The default behavior is just to
* subtract deltaU from uLastIter, i.e.
* \f[ u^{k+1} = u^k - \Delta u^k \f]
*
* \param assembler The assembler for Jacobian and residual
* \param uCurrentIter The solution vector 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.
*/
template<class JacobianAssembler, class SolutionVector>
void newtonUpdate(JacobianAssembler& assembler,
SolutionVector &uCurrentIter,
const SolutionVector &uLastIter,
const SolutionVector &deltaU)
{
if (this->enableShiftCriterion_)
this->newtonUpdateShift(uLastIter, deltaU);
if (this->useLineSearch_)
{
this->lineSearchUpdate_(assembler, uCurrentIter, uLastIter, deltaU);
}
else {
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_)
{
this->residualNorm_ = assembler.residualNorm(uCurrentIter);
this->reduction_ = this->residualNorm_;
this->reduction_ /= this->initialResidual_;
}
}
// update the variables class to the new solution
assembler.gridVariables().update(uCurrentIter);
}
/*!
* \brief Update the maximum relative shift of the solution compared to
* the previous iteration.
*
* \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;
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[dofTypeIdx][i], uNewI);
this->shift_ = std::max(this->shift_, shiftAtDof);
}
});
if (this->comm().size() > 1)
this->shift_ = this->comm().max(this->shift_);
}
private:
/*!
* \brief Solve the linear system of equations \f$\mathbf{A}x - b = 0\f$.
*
* Throws Dumux::NumericalProblem if the linear solver didn't
* converge.
*
* Specialization for linear solvers that can handle MultiType matrices.
*
*/
template<class LinearSolver, class JacobianMatrix, class SolutionVector>
bool solveLinearSystem_(LinearSolver& ls,
JacobianMatrix& A,
SolutionVector& x,
SolutionVector& b,
std::true_type)
{
// TODO: automatically derive the precondBlockLevel
return ls.template solve</*precondBlockLevel=*/2>(A, x, b);
}
/*!
* \brief Solve the linear system of equations \f$\mathbf{A}x - b = 0\f$.
*
* Throws Dumux::NumericalProblem if the linear solver didn't
* converge.
*
* Specialization for linear solvers that cannot handle MultiType matrices.
* We copy the matrix into a 1x1 block BCRS matrix before solving.
*
*/
template<class LinearSolver, class JacobianMatrix, class SolutionVector>
bool solveLinearSystem_(LinearSolver& ls,
JacobianMatrix& A,
SolutionVector& x,
SolutionVector& b,
std::false_type)
{
// create the bcrs matrix the IterativeSolver backend can handle
const auto M = MatrixConverter<JacobianMatrix>::multiTypeToBCRSMatrix(A);
// get the new matrix sizes
const std::size_t numRows = M.N();
assert(numRows == M.M());
// create the vector the IterativeSolver backend can handle
const auto bTmp = VectorConverter<SolutionVector>::multiTypeToBlockVector(b);
assert(bTmp.size() == numRows);
// create a blockvector to which the linear solver writes the solution
using VectorBlock = typename Dune::FieldVector<Scalar, 1>;
using BlockVector = typename Dune::BlockVector<VectorBlock>;
BlockVector y(numRows);
// solve
const bool converged = ls.solve(M, y, bTmp);
// copy back the result y into x
if(converged)
VectorConverter<SolutionVector>::retrieveValues(x, y);
return converged;
}
//! helper method to assure the MultiType matrix's sub blocks have the correct sizes
template<class JacobianMatrix>
bool checkMatrix_(const JacobianMatrix& A)
{
bool matrixHasCorrectSize = true;
using namespace Dune::Hybrid;
using namespace Dune::Indices;
forEach(A, [&matrixHasCorrectSize](const auto& rowOfMultiTypeMatrix)
{
const auto numRowsLeftMostBlock = rowOfMultiTypeMatrix[_0].N();
forEach(rowOfMultiTypeMatrix, [&matrixHasCorrectSize, &numRowsLeftMostBlock](const auto& subBlock)
{
if (subBlock.N() != numRowsLeftMostBlock)
matrixHasCorrectSize = false;
});
});
return matrixHasCorrectSize;
}
};
} // end namespace Dumux
#endif
dumux/porousmediumflow/2pnc/CMakeLists.txt
+ 0
− 1
View file @ d99d5076
...@@ -3,7 +3,6 @@ ...@@ -3,7 +3,6 @@
install(FILES install(FILES
indices.hh indices.hh
model.hh model.hh
newtoncontroller.hh
primaryvariableswitch.hh primaryvariableswitch.hh
volumevariables.hh volumevariables.hh
vtkoutputfields.hh vtkoutputfields.hh
......
dumux/porousmediumflow/compositional/privarswitchnewtoncontroller.hh deleted 100644 → 0
+ 0
− 220
View file @ 76995aed
// -*- 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 2 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 PorousmediumCompositional
* \brief Reference implementation of a controller class for the Newton solver.
*
* Usually this controller should be sufficient.
*/
#ifndef DUMUX_PRIVARSWITCH_NEWTON_CONTROLLER_HH
#define DUMUX_PRIVARSWITCH_NEWTON_CONTROLLER_HH
#include <memory>
#include <dumux/common/properties.hh>
#include <dumux/common/parameters.hh>
#include <dumux/discretization/methods.hh>
#include <dumux/discretization/elementsolution.hh>
#include <dumux/nonlinear/newtoncontroller.hh>
#include <dune/common/deprecated.hh>
#warning "This file is deprecated. Use PriVarSwitchNewtonSolver instead."
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 DUNE_DEPRECATED_MSG("Use PriVarSwitchNewtonSolver instead.")
PriVarSwitchNewtonController : public NewtonController<typename GET_PROP_TYPE(TypeTag, Scalar)>
{
using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
using ParentType = NewtonController<Scalar>;
using PrimaryVariableSwitch = typename GET_PROP_TYPE(TypeTag, PrimaryVariableSwitch);
static constexpr bool isBox = GET_PROP_TYPE(TypeTag, FVGridGeometry)::discMethod == DiscretizationMethod::box;
public:
using ParentType::ParentType;
/*!
* \brief Returns true if the error of the solution is below the
* tolerance.
*/
bool newtonConverged() const
{
if (switchedInLastIteration_)
return false;
return ParentType::newtonConverged();
}
/*!
*
* \brief Called before the Newton method is applied to an
* non-linear system of equations.
*
* \param u The initial solution
*/
template<class SolutionVector>
void newtonBegin(const SolutionVector &u)
{
ParentType::newtonBegin(u);
priVarSwitch_ = std::make_unique<PrimaryVariableSwitch>(u.size());
}
/*!
* \brief Indicates that one Newton iteration was finished.
*
* \param assembler The jacobian assembler
* \param uCurrentIter The solution after the current Newton iteration
* \param uLastIter The solution at the beginning of the current Newton iteration
*/
template<class JacobianAssembler, class SolutionVector>
void newtonEndStep(JacobianAssembler& assembler,
SolutionVector &uCurrentIter,
const SolutionVector &uLastIter)
{
ParentType::newtonEndStep(assembler, uCurrentIter, uLastIter);
// update the variable switch (returns true if the pri vars at at least one dof were switched)
// for disabled grid variable caching
const auto& fvGridGeometry = assembler.fvGridGeometry();
const auto& problem = assembler.problem();
auto& gridVariables = assembler.gridVariables();
// invoke the primary variable switch
switchedInLastIteration_ = priVarSwitch_->update(uCurrentIter, gridVariables,
problem, fvGridGeometry);
if(switchedInLastIteration_)
{
for (const auto& element : elements(fvGridGeometry.gridView()))
{
// if the volume variables are cached globally, we need to update those where the primary variables have been switched
updateSwitchedVolVars_(std::integral_constant<bool, GET_PROP_VALUE(TypeTag, EnableGridVolumeVariablesCache)>(),
element, assembler, uCurrentIter, uLastIter);
// if the flux variables are cached globally, we need to update those where the primary variables have been switched
// (not needed for box discretization)
updateSwitchedFluxVarsCache_(std::integral_constant<bool, (GET_PROP_VALUE(TypeTag, EnableGridFluxVariablesCache) && !isBox)>(),
element, assembler, uCurrentIter, uLastIter);
}
}
}
/*!
* \brief Called if the Newton method ended
* (not known yet if we failed or succeeded)
*/
void newtonEnd()
{
ParentType::newtonEnd();
// in any way, we have to reset the switch flag
switchedInLastIteration_ = false;
// free some memory
priVarSwitch_.release();
}
private:
/*!
* \brief Update the volume variables whose primary variables were
switched. Required when volume variables are cached globally.
*/
template<class Element, class JacobianAssembler, class SolutionVector>
void updateSwitchedVolVars_(std::true_type,
const Element& element,
JacobianAssembler& assembler,
const SolutionVector &uCurrentIter,
const SolutionVector &uLastIter)
{
const auto& fvGridGeometry = assembler.fvGridGeometry();
const auto& problem = assembler.problem();
auto& gridVariables = assembler.gridVariables();
// make sure FVElementGeometry is bound to the element
auto fvGeometry = localView(fvGridGeometry);
fvGeometry.bindElement(element);
// update the secondary variables if global caching is enabled
for (auto&& scv : scvs(fvGeometry))
{
const auto dofIdxGlobal = scv.dofIndex();
if (priVarSwitch_->wasSwitched(dofIdxGlobal))
{
const auto elemSol = elementSolution(element, uCurrentIter, fvGridGeometry);
auto& volVars = gridVariables.curGridVolVars().volVars(scv);
volVars.update(elemSol, problem, element, scv);
}
}
}
/*!
* \brief Update the fluxVars cache for dof whose primary variables were
switched. Required when flux variables are cached globally.
*/
template<class Element, class JacobianAssembler, class SolutionVector>
void updateSwitchedFluxVarsCache_(std::true_type,
const Element& element,
JacobianAssembler& assembler,
const SolutionVector &uCurrentIter,
const SolutionVector &uLastIter)
{
const auto& fvGridGeometry = assembler.fvGridGeometry();
auto& gridVariables = assembler.gridVariables();
// update the flux variables if global caching is enabled
const auto dofIdxGlobal = fvGridGeometry.dofMapper().index(element);
if (priVarSwitch_->wasSwitched(dofIdxGlobal))
{
// make sure FVElementGeometry and the volume variables are bound
auto fvGeometry = localView(fvGridGeometry);
fvGeometry.bind(element);
auto curElemVolVars = localView(gridVariables.curGridVolVars());
curElemVolVars.bind(element, fvGeometry, uCurrentIter);
gridVariables.gridFluxVarsCache().updateElement(element, fvGeometry, curElemVolVars);
}
}
//! brief Do nothing when volume variables are not cached globally.
template <typename... Args>
void updateSwitchedVolVars_(std::false_type, Args&&... args) const {}
//! brief Do nothing when flux variables are not cached globally.
template <typename... Args>
void updateSwitchedFluxVarsCache_(std::false_type, Args&&... args) const {}
//! the class handling the primary variable switch
std::unique_ptr<PrimaryVariableSwitch> priVarSwitch_;
//! if we switched primary variables in the last iteration
bool switchedInLastIteration_ = false;
};
} // end namespace Dumux
#endif
dumux/porousmediumflow/nonequilibrium/CMakeLists.txt
+ 1
− 1
View file @ d99d5076
...@@ -4,8 +4,8 @@ install(FILES ...@@ -4,8 +4,8 @@ install(FILES
localresidual.hh localresidual.hh
gridvariables.hh gridvariables.hh
indices.hh indices.hh
newtoncontroller.hh
volumevariables.hh volumevariables.hh
vtkoutputfields.hh vtkoutputfields.hh
model.hh model.hh
newtonsolver.hh
DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/dumux/porousmediumflow/nonequilibrium) DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/dumux/porousmediumflow/nonequilibrium)
dumux/porousmediumflow/nonequilibrium/newtoncontroller.hh deleted 100644 → 0
+ 0
− 66
View file @ 76995aed
// -*- 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 2 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 PorousmediumNonEquilibriumModel
* \brief A MpNc specific controller for the newton solver.
* This controller calls the velocity averaging in the model after each iteration.
*/
#ifndef DUMUX_NONEQUILIBRIUM_NEWTON_CONTROLLER_HH
#define DUMUX_NONEQUILIBRIUM_NEWTON_CONTROLLER_HH
#include <algorithm>
#include <dumux/nonlinear/newtoncontroller.hh>
namespace Dumux {
/*!
* \ingroup PorousmediumNonEquilibriumModel
* \brief A nonequilibrium specific controller for the newton solver.
* This controller calls the velocity averaging in the problem after each iteration.
*/
template <class Scalar,
class Comm = Dune::CollectiveCommunication<Dune::MPIHelper::MPICommunicator> >
class NonEquilibriumNewtonController : public NewtonController<Scalar, Comm>
{
using ParentType = NewtonController<Scalar, Comm>;
public:
using ParentType::ParentType;
template<class JacobianAssembler, class SolutionVector>
void newtonUpdate(JacobianAssembler& assembler,
SolutionVector &uCurrentIter,
const SolutionVector &uLastIter,
const SolutionVector &deltaU)
{
ParentType::newtonUpdate(assembler,
uCurrentIter,
uLastIter,
deltaU);
auto& gridVariables = assembler.gridVariables();
// Averages the face velocities of a vertex. Implemented in the model.
// The velocities are stored in the model.
gridVariables.calcVelocityAverage(uCurrentIter);
}
};
} // end namespace Dumux
#endif // DUMUX_VELO_PROB_AVERAGE_NEWTON_CONTROLLER_HH
dumux/porousmediumflow/richards/CMakeLists.txt
+ 2
− 1
View file @ d99d5076
...@@ -4,7 +4,8 @@ install(FILES ...@@ -4,7 +4,8 @@ install(FILES
indices.hh indices.hh
localresidual.hh localresidual.hh
model.hh model.hh
newtoncontroller.hh newtonsolver.hh
privarswitchnewtonsolver.hh
primaryvariableswitch.hh primaryvariableswitch.hh
volumevariables.hh volumevariables.hh
vtkoutputfields.hh vtkoutputfields.hh
......
dumux/porousmediumflow/richards/newtoncontroller.hh deleted 100644 → 0
+ 0
− 122
View file @ 76995aed
// -*- 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 2 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 RichardsModel
* \brief A Richards model specific controller for the newton solver.
*/
#ifndef DUMUX_RICHARDS_NEWTON_CONTROLLER_HH
#define DUMUX_RICHARDS_NEWTON_CONTROLLER_HH
#include <dumux/common/properties.hh>
#include <dumux/nonlinear/newtoncontroller.hh>
#include <dumux/discretization/elementsolution.hh>
#include <dune/common/deprecated.hh>
#warning "This file is deprecated. Use RichardsNewtonSolver instead."
namespace Dumux {
/*!
* \ingroup RichardsModel
* \brief A Richards model specific controller for the newton solver.
*
* 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 DUNE_DEPRECATED_MSG("Use RichardsNewtonSolver instead.")
RichardsNewtonController : public NewtonController<typename GET_PROP_TYPE(TypeTag, Scalar)>
{
using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
using ParentType = NewtonController<Scalar>;
using MaterialLaw = typename GET_PROP_TYPE(TypeTag, MaterialLaw);
using Indices = typename GET_PROP_TYPE(TypeTag, Indices);
enum { pressureIdx = Indices::pressureIdx };
public:
using ParentType::ParentType;
/*!
* \brief Update the current solution of the newton method
*
* This is basically the step
* \f[ u^{k+1} = u^k - \Delta u^k \f]
*
* \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
* iterative solution and the next one \f$ \Delta u^k \f$
*/
template<class JacobianAssembler, class SolutionVector>
void newtonUpdate(JacobianAssembler& assembler,
SolutionVector &uCurrentIter,
const SolutionVector &uLastIter,
const SolutionVector &deltaU)
{
ParentType::newtonUpdate(assembler, uCurrentIter, uLastIter, deltaU);
if (!this->useLineSearch_ && getParamFromGroup<bool>(this->paramGroup(), "Newton.EnableChop"))
{
// do not clamp anything after 5 iterations
if (this->numSteps_ > 4)
return;
// clamp saturation change to at most 20% per iteration
const auto& fvGridGeometry = assembler.fvGridGeometry();
for (const auto& element : elements(fvGridGeometry.gridView()))
{
auto fvGeometry = localView(fvGridGeometry);
fvGeometry.bindElement(element);
for (auto&& scv : scvs(fvGeometry))
{
auto dofIdxGlobal = scv.dofIndex();
// calculate the old wetting phase saturation
const auto& spatialParams = assembler.problem().spatialParams();
const auto elemSol = elementSolution(element, uCurrentIter, fvGridGeometry);
const auto& materialLawParams = spatialParams.materialLawParams(element, scv, elemSol);
const Scalar pcMin = MaterialLaw::pc(materialLawParams, 1.0);
const Scalar pw = uLastIter[dofIdxGlobal][pressureIdx];
using std::max;
const Scalar pn = max(assembler.problem().nonWettingReferencePressure(), pw + pcMin);
const Scalar pcOld = pn - pw;
const Scalar SwOld = max(0.0, MaterialLaw::sw(materialLawParams, pcOld));
// convert into minimum and maximum wetting phase
// pressures
const Scalar pwMin = pn - MaterialLaw::pc(materialLawParams, SwOld - 0.2);
const Scalar pwMax = pn - MaterialLaw::pc(materialLawParams, SwOld + 0.2);
// clamp the result
using std::min; using std::max;
uCurrentIter[dofIdxGlobal][pressureIdx] = max(pwMin, min(uCurrentIter[dofIdxGlobal][pressureIdx], pwMax));
}
}
}
}
};
} // end namespace Dumux
#endif
dumux/porousmediumflow/richardsnc/model.hh
+ 0
− 1
View file @ d99d5076
...@@ -70,7 +70,6 @@ ...@@ -70,7 +70,6 @@
#include <dumux/common/properties.hh> #include <dumux/common/properties.hh>
#include <dumux/porousmediumflow/compositional/localresidual.hh> #include <dumux/porousmediumflow/compositional/localresidual.hh>
#include <dumux/porousmediumflow/richards/newtoncontroller.hh>
#include <dumux/material/spatialparams/fv1p.hh> #include <dumux/material/spatialparams/fv1p.hh>
#include <dumux/material/fluidmatrixinteractions/diffusivitymillingtonquirk.hh> #include <dumux/material/fluidmatrixinteractions/diffusivitymillingtonquirk.hh>
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment