diff --git a/dumux/nonlinear/newtoncontroller.hh b/dumux/nonlinear/newtoncontroller.hh
index c2bf0118a458f77e0208cbffb658f79dee93ce0e..a611e62a4d321898fcef76fc432322deed9c136e 100644
--- a/dumux/nonlinear/newtoncontroller.hh
+++ b/dumux/nonlinear/newtoncontroller.hh
@@ -56,12 +56,6 @@ NEW_PROP_TAG(SolutionVector);
//! Specifies the type of a global Jacobian matrix
NEW_PROP_TAG(JacobianMatrix);
-//! Specifies the type of the Vertex mapper
-NEW_PROP_TAG(VertexMapper);
-
-//! specifies the type of the time manager
-NEW_PROP_TAG(TimeManager);
-
//! specifies whether the convergence rate and the global residual
//! gets written out to disk for every Newton iteration (default is false)
NEW_PROP_TAG(NewtonWriteConvergence);
@@ -155,39 +149,30 @@ SET_INT_PROP(NewtonMethod, NewtonMaxSteps, 18);
template <class TypeTag>
class NewtonController
{
- typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
- typedef typename GET_PROP_TYPE(TypeTag, NewtonController) Implementation;
- typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
-
- typedef typename GET_PROP_TYPE(TypeTag, Problem) Problem;
- typedef typename GET_PROP_TYPE(TypeTag, Model) Model;
- typedef typename GET_PROP_TYPE(TypeTag, NewtonMethod) NewtonMethod;
- typedef typename GET_PROP_TYPE(TypeTag, JacobianMatrix) JacobianMatrix;
- typedef typename GET_PROP_TYPE(TypeTag, JacobianAssembler) JacobianAssembler;
- typedef typename GET_PROP_TYPE(TypeTag, TimeManager) TimeManager;
- typedef typename GET_PROP_TYPE(TypeTag, VertexMapper) VertexMapper;
+ using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
+ using Implementation = typename GET_PROP_TYPE(TypeTag, NewtonController);
+ using ConvergenceWriter = typename GET_PROP_TYPE(TypeTag, NewtonConvergenceWriter);
+ using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
+ using Communicator = typename GridView::CollectiveCommunication;
using NumEqVector = typename GET_PROP_TYPE(TypeTag, NumEqVector);
-
- typedef typename GET_PROP_TYPE(TypeTag, SolutionVector) SolutionVector;
-
- typedef typename GET_PROP_TYPE(TypeTag, NewtonConvergenceWriter) ConvergenceWriter;
-
- typedef typename GET_PROP_TYPE(TypeTag, LinearSolver) LinearSolver;
+ using SolutionVector = typename GET_PROP_TYPE(TypeTag, SolutionVector);
+ using JacobianMatrix = typename GET_PROP_TYPE(TypeTag, JacobianMatrix);
+ using JacobianAssembler = typename GET_PROP_TYPE(TypeTag, JacobianAssembler);
+ using LinearSolver = typename GET_PROP_TYPE(TypeTag, LinearSolver);
static constexpr int numEq = GET_PROP_VALUE(TypeTag, NumEq);
public:
/*!
- * \brief Constructor
+ * \brief Constructor without convergence writer
*/
- NewtonController(const Problem &problem)
- : endIterMsgStream_(std::ostringstream::out),
- linearSolver_(problem)
+ NewtonController(const Communicator& comm) : comm_(comm), endIterMsgStream_(std::ostringstream::out)
{
enablePartialReassemble_ = GET_PARAM_FROM_GROUP(TypeTag, bool, Implicit, EnablePartialReassemble);
enableJacobianRecycling_ = GET_PARAM_FROM_GROUP(TypeTag, bool, Implicit, EnableJacobianRecycling);
+ writeConvergence_ = GET_PARAM_FROM_GROUP(TypeTag, bool, Newton, WriteConvergence);
useLineSearch_ = GET_PARAM_FROM_GROUP(TypeTag, bool, Newton, UseLineSearch);
enableAbsoluteResidualCriterion_ = GET_PARAM_FROM_GROUP(TypeTag, bool, Newton, EnableAbsoluteResidualCriterion);
enableShiftCriterion_ = GET_PARAM_FROM_GROUP(TypeTag, bool, Newton, EnableShiftCriterion);
@@ -211,6 +196,13 @@ public:
numSteps_ = 0;
}
+ /*!
+ * \brief Constructor with a given convergence writer
+ */
+ NewtonController(const Communicator& comm, ConvergenceWriter&& writer)
+ : NewtonController(comm), convergenceWriter_(std::move(writer))
+ {}
+
/*!
* \brief Set the maximum acceptable difference of any primary variable
* between two iterations for declaring convergence.
@@ -326,20 +318,14 @@ public:
* \brief Called before the Newton method is applied to an
* non-linear system of equations.
*
- * \param method The object where the NewtonMethod is executed
* \param u The initial solution
*/
- void newtonBegin(NewtonMethod &method, const SolutionVector &u)
+ void newtonBegin(const SolutionVector &u)
{
- method_ = &method;
numSteps_ = 0;
- if (GET_PARAM_FROM_GROUP(TypeTag, bool, Newton, WriteConvergence))
- {
- if (!convergenceWriter_)
- convergenceWriter_ = std::make_unique<ConvergenceWriter>(asImp_(), problem_().gridView());
+ if (writeConvergence_)
convergenceWriter_->advanceTimeStep();
- }
}
/*!
@@ -381,14 +367,23 @@ public:
typename SolutionVector::block_type uNewI = uLastIter[i];
uNewI -= deltaU[i];
- Scalar shiftAtDof = model_().relativeShiftAtDof(uLastIter[i],
- uNewI);
+ Scalar shiftAtDof = relativeShiftAtDof_(uLastIter[i], uNewI);
using std::max;
shift_ = max(shift_, shiftAtDof);
}
- if (gridView_().comm().size() > 1)
- shift_ = gridView_().comm().max(shift_);
+ if (comm().size() > 1)
+ shift_ = comm().max(shift_);
+ }
+
+ /*!
+ * \brief Assemble the linear system of equations \f$\mathbf{A}x - b = 0\f$.
+ *
+ * \param assembler The jacobian assembler
+ */
+ void assembleLinearSystem(JacobianAssembler& assembler)
+ {
+ assembler.assembleJacobianAndResidual();
}
/*!
@@ -397,20 +392,22 @@ public:
* Throws NumericalProblem if the linear solver didn't
* converge.
*
+ * \param ls The linear solver to be used
* \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
*/
- void newtonSolveLinear(JacobianMatrix &A,
- SolutionVector &x,
- SolutionVector &b)
+ void solveLinearSystem(LinearSolver& ls,
+ JacobianMatrix& A,
+ SolutionVector& x,
+ SolutionVector& b)
{
try {
if (numSteps_ == 0)
{
Scalar norm2 = b.two_norm2();
- if (gridView_().comm().size() > 1)
- norm2 = gridView_().comm().sum(norm2);
+ if (comm().size() > 1)
+ norm2 = comm().sum(norm2);
using std::sqrt;
initialResidual_ = sqrt(norm2);
@@ -427,7 +424,7 @@ public:
for (int j = 0; j < numEq; ++j)
bTmp[i][j] = b[i][j];
- int converged = linearSolver_.solve(A, xTmp, bTmp);
+ int converged = ls.solve(A, xTmp, bTmp);
for (int i = 0; i < x.size(); ++i)
for (int j = 0; j < numEq; ++j)
@@ -435,8 +432,8 @@ public:
// make sure all processes converged
int convergedRemote = converged;
- if (gridView_().comm().size() > 1)
- convergedRemote = gridView_().comm().min(converged);
+ if (comm().size() > 1)
+ convergedRemote = comm().min(converged);
if (!converged) {
DUNE_THROW(NumericalProblem,
@@ -450,8 +447,8 @@ public:
catch (Dune::MatrixBlockError e) {
// make sure all processes converged
int converged = 0;
- if (gridView_().comm().size() > 1)
- converged = gridView_().comm().min(converged);
+ if (comm().size() > 1)
+ converged = comm().min(converged);
NumericalProblem p;
std::string msg;
@@ -463,8 +460,8 @@ public:
catch (const Dune::Exception &e) {
// make sure all processes converged
int converged = 0;
- if (gridView_().comm().size() > 1)
- converged = gridView_().comm().min(converged);
+ if (comm().size() > 1)
+ converged = comm().min(converged);
NumericalProblem p;
p.message(e.what());
@@ -483,13 +480,15 @@ public:
* subtract deltaU from uLastIter, i.e.
* \f[ u^{k+1} = u^k - \Delta u^k \f]
*
+ * \param assembler The assembler (needed for global residual evaluation)
* \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.
*/
- void newtonUpdate(SolutionVector &uCurrentIter,
+ void newtonUpdate(const JacobianAssembler& assembler,
+ SolutionVector &uCurrentIter,
const SolutionVector &uLastIter,
const SolutionVector &deltaU)
{
@@ -500,7 +499,7 @@ public:
if (useLineSearch_)
{
- lineSearchUpdate_(uCurrentIter, uLastIter, deltaU);
+ lineSearchUpdate_(assembler, uCurrentIter, uLastIter, deltaU);
}
else {
for (unsigned int i = 0; i < uLastIter.size(); ++i) {
@@ -510,8 +509,15 @@ public:
if (enableResidualCriterion_)
{
- SolutionVector tmp(uLastIter);
- residual_ = this->method().model().globalResidual(tmp, uCurrentIter);
+ // Originally here we handed in uCurrentIter into
+ // the global residual function to evaluate the
+ // residual for the given solution. Currently, the
+ // assembler evaluates the residual for the curSol
+ // container only. However, uCurrentIter is a ref
+ // to this container, so we don't need to do the copying
+ // etc.
+ // TODO: Should we re-include this??????
+ residual_ = assembler.globalResidual();
reduction_ = residual_;
reduction_ /= initialResidual_;
}
@@ -521,14 +527,24 @@ public:
/*!
* \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
*/
- void newtonEndStep(SolutionVector &uCurrentIter,
+ void newtonEndStep(JacobianAssembler& assembler,
+ SolutionVector &uCurrentIter,
const SolutionVector &uLastIter)
{
- // Update the volume variables
- this->model_().newtonEndStep();
+ // TODO
+ // This is a hack in order to be able to realize an update
+ // of the gridVolumeVariables etc... How could we achieve
+ // this differently??
+ // Originally here we had model.newtonEndStep() which updated
+ // the variables... Should the NewtonMethod get access to these
+ // classes directly as well?
+ // -> The assembler has all the data w.r.t. the problem etc...
+ // -> also, how do we call this?
+ assembler.updateVariables()
++numSteps_;
@@ -546,7 +562,8 @@ public:
endIterMsgStream_.str("");
// When the Newton iterations are done: ask the model to check whether it makes sense
- model_().checkPlausibility();
+ // TODO: how do we realize this?
+ // model_().checkPlausibility();
}
/*!
@@ -598,20 +615,6 @@ public:
return oldTimeStep*(1.0 + percent/1.2);
}
- /*!
- * \brief Returns a reference to the current Newton method
- * which is controlled by this controller.
- */
- NewtonMethod &method()
- { return *method_; }
-
- /*!
- * \brief Returns a reference to the current Newton method
- * which is controlled by this controller.
- */
- const NewtonMethod &method() const
- { return *method_; }
-
std::ostringstream &endIterMsg()
{ return endIterMsgStream_; }
@@ -625,56 +628,9 @@ public:
* \brief Returns true if the Newton method ought to be chatty.
*/
bool verbose() const
- { return verbose_ && gridView_().comm().rank() == 0; }
+ { return verbose_ && comm().rank() == 0; }
protected:
- /*!
- * \brief Returns a reference to the grid view.
- */
- const GridView &gridView_() const
- { return problem_().gridView(); }
-
- /*!
- * \brief Returns a reference to the vertex mapper.
- */
- const VertexMapper &vertexMapper_() const
- { return model_().vertexMapper(); }
-
- /*!
- * \brief Returns a reference to the problem.
- */
- Problem &problem_()
- { return method_->problem(); }
-
- /*!
- * \brief Returns a reference to the problem.
- */
- const Problem &problem_() const
- { return method_->problem(); }
-
- /*!
- * \brief Returns a reference to the time manager.
- */
- TimeManager &timeManager_()
- { return problem_().timeManager(); }
-
- /*!
- * \brief Returns a reference to the time manager.
- */
- const TimeManager &timeManager_() const
- { return problem_().timeManager(); }
-
- /*!
- * \brief Returns a reference to the problem.
- */
- Model &model_()
- { return problem_().model(); }
-
- /*!
- * \brief Returns a reference to the problem.
- */
- const Model &model_() const
- { return problem_().model(); }
// returns the actual implementation for the controller we do
// it this way in order to allow "poor man's virtual methods",
@@ -688,14 +644,15 @@ protected:
void writeConvergence_(const SolutionVector &uLastIter,
const SolutionVector &deltaU)
{
- if (GET_PARAM_FROM_GROUP(TypeTag, bool, Newton, WriteConvergence))
+ if (writeConvergence_)
{
convergenceWriter_->advanceIteration();
convergenceWriter_->write(uLastIter, deltaU);
}
}
- void lineSearchUpdate_(SolutionVector &uCurrentIter,
+ void lineSearchUpdate_(const JacobianAssembler& assembler,
+ SolutionVector &uCurrentIter,
const SolutionVector &uLastIter,
const SolutionVector &deltaU)
{
@@ -708,12 +665,20 @@ protected:
uCurrentIter += uLastIter;
// calculate the residual of the current solution
- residual_ = this->method().model().globalResidual(tmp, uCurrentIter);
+ // Originally here we handed in uCurrentIter into
+ // the global residual function to evaluate the
+ // residual for the given solution. Currently, the
+ // assembler evaluates the residual for the curSol
+ // container only. However, uCurrentIter is a ref
+ // to this container, so we don't need to do the copying
+ // etc.
+ // TODO: Should we re-include this??????
+ residual_ = assembler.globalResidual();
reduction_ = residual_;
reduction_ /= initialResidual_;
if (reduction_ < lastReduction_ || lambda <= 0.125) {
- this->endIterMsg() << ", residual reduction " << lastReduction_ << "->" << reduction_ << "@lambda=" << lambda;
+ endIterMsg() << ", residual reduction " << lastReduction_ << "->" << reduction_ << "@lambda=" << lambda;
return;
}
@@ -722,12 +687,39 @@ protected:
}
}
+ /*!
+ * \brief Returns the maximum relative shift between two vectors of
+ * primary variables.
+ *
+ * \param priVars1 The first vector of primary variables
+ * \param priVars2 The second vector of primary variables
+ */
+ Scalar relativeShiftAtDof_(const PrimaryVariables &priVars1,
+ const PrimaryVariables &priVars2)
+ {
+ Scalar result = 0.0;
+ using std::abs;
+ using std::max;
+ for (int j = 0; j < numEq; ++j) {
+ Scalar eqErr = abs(priVars1[j] - priVars2[j]);
+ eqErr /= max<Scalar>(1.0,abs(priVars1[j] + priVars2[j])/2);
+
+ result = max(result, eqErr);
+ }
+ return result;
+ }
+
+ // The grid view's communicator
+ const Communicator& comm_;
+
+ // message stream to be displayed at the end of iterations
std::ostringstream endIterMsgStream_;
- NewtonMethod *method_;
- bool verbose_;
+ // writes an output file for each iteration
+ ConvergenceWriter convergenceWriter_;
- std::unique_ptr<ConvergenceWriter> convergenceWriter_;
+ // switches on/off verbosity
+ bool verbose_;
// shift criterion variables
Scalar shift_;
@@ -749,9 +741,7 @@ protected:
// actual number of steps done so far
int numSteps_;
- // the linear solver
- LinearSolver linearSolver_;
-
+ // further parameters
bool enablePartialReassemble_;
bool enableJacobianRecycling_;
bool useLineSearch_;
@@ -759,6 +749,9 @@ protected:
bool enableShiftCriterion_;
bool enableResidualCriterion_;
bool satisfyResidualAndShiftCriterion_;
+
+ // specifies if we want to write down the intermediate solutions
+ bool writeConvergence_;
};
} // namespace Dumux
diff --git a/dumux/nonlinear/newtonconvergencewriter.hh b/dumux/nonlinear/newtonconvergencewriter.hh
index d0dcd8241fcfeb4d5b432b219b8db24066ffce2d..2ef85f6a5e6e1eba3f81427cf2e716e02a8d885c 100644
--- a/dumux/nonlinear/newtonconvergencewriter.hh
+++ b/dumux/nonlinear/newtonconvergencewriter.hh
@@ -29,17 +29,9 @@
#include <dumux/common/basicproperties.hh>
-#include "newtoncontroller.hh"
-
namespace Dumux
{
-namespace Properties
-{
-NEW_PROP_TAG(NewtonController);
-NEW_PROP_TAG(SolutionVector);
-}
-
/*!
* \ingroup Newton
* \brief Writes the intermediate solutions during
@@ -50,21 +42,19 @@ class NewtonConvergenceWriter
{
using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
- using NewtonController = typename GET_PROP_TYPE(TypeTag, NewtonController);
+ using LocalResidual = typename GET_PROP_TYPE(TypeTag, LocalResidual);
using SolutionVector = typename GET_PROP_TYPE(TypeTag, SolutionVector);
static constexpr int numEq = GET_PROP_VALUE(TypeTag, NumEq);
static constexpr bool isBox = GET_PROP_VALUE(TypeTag, ImplicitIsBox);
public:
- NewtonConvergenceWriter(NewtonController &ctl, const GridView& gridView)
- : ctl_(ctl),
- writer_(gridView, "convergence", "", "")
+ NewtonConvergenceWriter(const GridView& gridView, const std::size_t numDofs)
+ : writer_(gridView, "convergence", "", "")
{
timeStepIndex_ = 0;
iteration_ = 0;
- const auto numDofs = ctl_.method().model().numDofs();
for (int eqIdx = 0; eqIdx < numEq; ++eqIdx)
{
def_[eqIdx].resize(numDofs);
@@ -106,9 +96,7 @@ public:
void write(const SolutionVector &uLastIter,
const SolutionVector &deltaU)
{
- SolutionVector residual(uLastIter);
- ctl_.method().model().globalResidual(residual, uLastIter);
-
+ const auto residual = LocalResidual::globalResidual();
for (unsigned int dofIdxGlobal = 0; dofIdxGlobal < deltaU.size(); dofIdxGlobal++)
{
for (int eqIdx = 0; eqIdx < numEq; ++eqIdx)
@@ -130,7 +118,6 @@ private:
std::array<std::vector<Scalar>, numEq> delta_;
std::array<std::vector<Scalar>, numEq> x_;
- NewtonController &ctl_;
Dune::VTKSequenceWriter<GridView> writer_;
};
diff --git a/dumux/nonlinear/newtonmethod.hh b/dumux/nonlinear/newtonmethod.hh
index 55b5531036f4b9613929dd76505af9d57bc827b9..0b142216c5f5613043ca94ad76eca74190f204c1 100644
--- a/dumux/nonlinear/newtonmethod.hh
+++ b/dumux/nonlinear/newtonmethod.hh
@@ -43,8 +43,6 @@ namespace Properties
NEW_TYPE_TAG(NewtonMethod);
NEW_PROP_TAG(Scalar);
-NEW_PROP_TAG(Problem);
-NEW_PROP_TAG(Model);
NEW_PROP_TAG(NewtonController);
NEW_PROP_TAG(SolutionVector);
NEW_PROP_TAG(JacobianAssembler);
@@ -59,51 +57,42 @@ NEW_PROP_TAG(JacobianAssembler);
template <class TypeTag>
class NewtonMethod
{
- typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
- typedef typename GET_PROP_TYPE(TypeTag, Problem) Problem;
- typedef typename GET_PROP_TYPE(TypeTag, Model) Model;
- typedef typename GET_PROP_TYPE(TypeTag, NewtonController) NewtonController;
+ using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
+ using SolutionVector = typename GET_PROP_TYPE(TypeTag, SolutionVector);
+ using NewtonController = typename GET_PROP_TYPE(TypeTag, NewtonController);
+ using ConvergenceWriter = typename GET_PROP_TYPE(TypeTag, NewtonConvergenceWriter);
+ using JacobianAssembler = typename GET_PROP_TYPE(TypeTag, JacobianAssembler);
+ using LinearSolver = typename GET_PROP_TYPE(TypeTag, LinearSolver);
- typedef typename GET_PROP_TYPE(TypeTag, SolutionVector) SolutionVector;
- typedef typename GET_PROP_TYPE(TypeTag, JacobianAssembler) JacobianAssembler;
public:
- NewtonMethod(Problem &problem)
- : problem_(problem)
- { }
-
- /*!
- * \brief Returns a reference to the current numeric problem.
- */
- Problem &problem()
- { return problem_; }
-
- /*!
- * \brief Returns a reference to the current numeric problem.
- */
- const Problem &problem() const
- { return problem_; }
-
- /*!
- * \brief Returns a reference to the numeric model.
- */
- Model &model()
- { return problem().model(); }
-
- /*!
- * \brief Returns a reference to the numeric model.
- */
- const Model &model() const
- { return problem().model(); }
+ NewtonMethod(std::shared_ptr<JacobianAssembler> jacobianAssembler,
+ std::shared_ptr<LinearSolver> linearSolver)
+ : jacobianAssembler_(jacobianAssembler)
+ , linearSolver_(linearSolver)
+ , matrix_ = std::make_shared<JacobianMatrix>()
+ , residual_ = std::make_shared<SolutionVector>()
+ {
+ // construct the newton controller with or without convergence writer
+ if (GET_PARAM_FROM_GROUP(TypeTag, bool, Newton, WriteConvergence))
+ {
+ newtonController_ = std::make_unique<NewtonController>(assembler().gridView().comm(),
+ NewtonConvergenceWriter(gridView, assembler().numDofs()));
+ }
+ else
+ newtonController_ = std::make_unique<NewtonController>(assembler().gridView().comm());
+ // set the linear system (matrix & residual) in the assembler
+ assembler().setLinearSystem(matrix(), residual());
+ }
/*!
* \brief Run the newton method. The controller is responsible
* for all the strategic decisions.
*/
- bool execute(NewtonController &ctl)
+ bool execute()
{
try {
- return execute_(ctl);
+ return execute_();
}
catch (const NumericalProblem &e) {
if (ctl.verbose())
@@ -113,32 +102,47 @@ public:
}
}
-protected:
- bool execute_(NewtonController &ctl)
+ NewtonController& controller()
+ { return *newtonController_; }
+
+ JacobianAssembler& assembler()
+ { return *jacobianAssembler_; }
+
+ LinearSolver& linearSolver()
+ { return *linearSolver_; }
+
+ JacobianMatrix& matrix()
+ { return *matrix_; }
+
+ SolutionVector& residual()
+ { return residual_; }
+
+private:
+ bool execute_()
{
- SolutionVector &uCurrentIter = model().curSol();
+ // the current solution is the initial guess
+ SolutionVector& uCurrentIter = assembler.curSol();
SolutionVector uLastIter(uCurrentIter);
SolutionVector deltaU(uCurrentIter);
- JacobianAssembler &jacobianAsm = model().jacobianAssembler();
-
Dune::Timer assembleTimer(false);
Dune::Timer solveTimer(false);
Dune::Timer updateTimer(false);
// tell the controller that we begin solving
- ctl.newtonBegin(*this, uCurrentIter);
+ controller().newtonBegin(uCurrentIter);
// execute the method as long as the controller thinks
// that we should do another iteration
- while (ctl.newtonProceed(uCurrentIter))
+ while (controller().newtonProceed(uCurrentIter))
{
// notify the controller that we're about to start
// a new timestep
- ctl.newtonBeginStep();
+ controller().newtonBeginStep();
// make the current solution to the old one
- uLastIter = uCurrentIter;
+ if (controller().newtonNumSteps() > 0)
+ uLastIter = uCurrentIter;
if (ctl.verbose()) {
std::cout << "Assemble: r(x^k) = dS/dt + div F - q; M = grad r";
@@ -151,7 +155,7 @@ protected:
// linearize the problem at the current solution
assembleTimer.start();
- jacobianAsm.assemble();
+ controller().assembleLinearSystem(assembler());
assembleTimer.stop();
///////////////
@@ -175,9 +179,10 @@ protected:
// set the delta vector to zero before solving the linear system!
deltaU = 0;
// ask the controller to solve the linearized system
- ctl.newtonSolveLinear(jacobianAsm.matrix(),
- deltaU,
- jacobianAsm.residual());
+ controller().solveLinearSystem(linearSolver(),
+ jacobianAsm.matrix(),
+ deltaU,
+ jacobianAsm.residual());
solveTimer.stop();
///////////////
@@ -192,11 +197,11 @@ protected:
updateTimer.start();
// update the current solution (i.e. uOld) with the delta
// (i.e. u). The result is stored in u
- ctl.newtonUpdate(uCurrentIter, uLastIter, deltaU);
+ ctl.newtonUpdate(assembler(), uCurrentIter, uLastIter, deltaU);
updateTimer.stop();
// tell the controller that we're done with this iteration
- ctl.newtonEndStep(uCurrentIter, uLastIter);
+ ctl.newtonEndStep(assembler(), uCurrentIter, uLastIter);
}
// tell the controller that we're done
@@ -220,8 +225,13 @@ protected:
return true;
}
-private:
- Problem &problem_;
+ std::shared_ptr<JacobianAssembler> jacobianAssembler_;
+ std::shared_ptr<LinearSolver> linearSolver_;
+
+ std::unique_ptr<NewtonController> newtonController_;
+
+ std::shared_ptr<JacobianMatrix> matrix_;
+ std::shared_ptr<SolutionVector> residual_;
};
}