newtoncontroller.hh 20.2 KB
Newer Older
Andreas Lauser's avatar
Andreas Lauser committed
1
// $Id$
2
/****************************************************************************
Andreas Lauser's avatar
Andreas Lauser committed
3
4
5
6
7
8
9
10
11
12
13
14
15
16
 *   Copyright (C) 2008-2010 by Andreas Lauser                               *
 *   Copyright (C) 2008-2010 by Bernd Flemisch                               *
 *   Institute of Hydraulic Engineering                                      *
 *   University of Stuttgart, Germany                                        *
 *   email: <givenname>.<name>@iws.uni-stuttgart.de                          *
 *                                                                           *
 *   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, as long as this copyright notice    *
 *   is included in its original form.                                       *
 *                                                                           *
 *   This program is distributed WITHOUT ANY WARRANTY.                       *
 *****************************************************************************/
17
/*!
Andreas Lauser's avatar
Andreas Lauser committed
18
19
20
21
22
 * \file
 * \brief Reference implementation of a newton controller solver.
 *
 * Usually for most cases this controller should be sufficient.
 */
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
#ifndef DUMUX_NEWTON_CONTROLLER_HH
#define DUMUX_NEWTON_CONTROLLER_HH

#include <dumux/common/exceptions.hh>

#include <dune/istl/overlappingschwarz.hh>
#include <dune/istl/schwarz.hh>
#include <dune/istl/preconditioners.hh>
#include <dune/istl/solvers.hh>
#include "dune/istl/owneroverlapcopy.hh"

#include <dune/istl/io.hh>

#include <dune/common/mpihelper.hh>

#include <iostream>
#include <boost/format.hpp>

#include <dumux/common/pardiso.hh>

Andreas Lauser's avatar
Andreas Lauser committed
43
#include <dumux/common/pdelabpreconditioner.hh>
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
#include <dune/pdelab/backend/istlsolverbackend.hh>


namespace Dumux
{
namespace Properties
{
//! specifies the verbosity of the linear solver (by default it is 0,
//! i.e. it doesn't print anything)
NEW_PROP_TAG(NewtonLinearSolverVerbosity);

//! 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);

//! specifies whether the update should be done using the line search
//! method instead of the "raw" newton method. whether this property
//! has any effect depends on wether the line search method is
//! implemented for the actual model's newton controller's update()
//! method. By default we do not use line search.
NEW_PROP_TAG(NewtonUseLineSearch);

SET_PROP_DEFAULT(NewtonLinearSolverVerbosity)
{public:
    static const int value = 0;
};

SET_PROP_DEFAULT(NewtonWriteConvergence)
{public:
    static const bool value = false;
};

SET_PROP_DEFAULT(NewtonUseLineSearch)
{public:
    static const bool value = false;
};
};


template <class TypeTag, bool enable>
struct NewtonConvergenceWriter
{
    typedef typename GET_PROP_TYPE(TypeTag, PTAG(GridView)) GridView;
    typedef typename GET_PROP_TYPE(TypeTag, PTAG(NewtonController)) NewtonController;

Andreas Lauser's avatar
Andreas Lauser committed
89
    typedef typename GET_PROP_TYPE(TypeTag, PTAG(SolutionVector)) SolutionVector;
90
91
92
    typedef Dumux::VtkMultiWriter<GridView>  VtkMultiWriter;

    NewtonConvergenceWriter(NewtonController &ctl)
Andreas Lauser's avatar
Andreas Lauser committed
93
        : ctl_(ctl)
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
    {
        timeStepNum_ = 0;
        iteration_ = 0;
        vtkMultiWriter_ = new VtkMultiWriter("convergence");
    }

    ~NewtonConvergenceWriter()
    { delete vtkMultiWriter_; };

    void beginTimestep()
    {
        ++timeStepNum_;
        iteration_ = 0;
    };

    void beginIteration(const GridView &gv)
    {
        ++ iteration_;
        vtkMultiWriter_->beginTimestep(timeStepNum_ + iteration_ / 100.0,
Andreas Lauser's avatar
Andreas Lauser committed
113
                                       gv);
114
115
116
117
118
    };

    void writeFields(const SolutionVector &uOld,
                     const SolutionVector &deltaU)
    {
Andreas Lauser's avatar
Andreas Lauser committed
119
        ctl_.method().model().addConvergenceVtkFields(*vtkMultiWriter_, uOld, deltaU);
120
121
122
    };

    void endIteration()
Andreas Lauser's avatar
Andreas Lauser committed
123
    { vtkMultiWriter_->endTimestep(); };
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142

    void endTimestep()
    {
        ++timeStepNum_;
        iteration_ = 0;
    };

private:
    int timeStepNum_;
    int iteration_;
    VtkMultiWriter *vtkMultiWriter_;
    NewtonController &ctl_;
};

template <class TypeTag>
struct NewtonConvergenceWriter<TypeTag, false>
{
    typedef typename GET_PROP_TYPE(TypeTag, PTAG(GridView)) GridView;
    typedef typename GET_PROP_TYPE(TypeTag, PTAG(NewtonController)) NewtonController;
Andreas Lauser's avatar
Andreas Lauser committed
143
    typedef typename GET_PROP_TYPE(TypeTag, PTAG(SolutionVector)) SolutionVector;
144
145
146
147
148
149
150
151
152
153
154
155
156

    typedef Dumux::VtkMultiWriter<GridView>  VtkMultiWriter;

    NewtonConvergenceWriter(NewtonController &ctl)
    {};

    void beginTimestep()
    { };

    void beginIteration(const GridView &gv)
    { };

    void writeFields(const SolutionVector &uOld,
Andreas Lauser's avatar
Andreas Lauser committed
157
                     const SolutionVector &deltaU)
158
159
160
161
162
163
164
165
166
167
    { };

    void endIteration()
    { };

    void endTimestep()
    { };
};

/*!
Andreas Lauser's avatar
Andreas Lauser committed
168
169
170
171
172
173
174
 * \brief The reference implementation of a newton controller.
 *
 * 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.
 */
175
176
177
template <class TypeTag>
class NewtonController
{
Andreas Lauser's avatar
Andreas Lauser committed
178
179
180
    typedef typename GET_PROP_TYPE(TypeTag, PTAG(Scalar)) Scalar;
    typedef typename GET_PROP_TYPE(TypeTag, PTAG(NewtonController)) Implementation;
    typedef typename GET_PROP_TYPE(TypeTag, PTAG(GridView)) GridView;
181

Andreas Lauser's avatar
Andreas Lauser committed
182
183
    typedef typename GET_PROP_TYPE(TypeTag, PTAG(Problem)) Problem;
    typedef typename GET_PROP_TYPE(TypeTag, PTAG(Model)) Model;
184
185
    typedef typename GET_PROP_TYPE(TypeTag, PTAG(NewtonMethod)) NewtonMethod;

Andreas Lauser's avatar
Andreas Lauser committed
186
187
188
    typedef typename GET_PROP_TYPE(TypeTag, PTAG(GridFunctionSpace)) GridFunctionSpace;
    typedef typename GET_PROP_TYPE(TypeTag, PTAG(ConstraintsTrafo)) ConstraintsTrafo;
    typedef typename GET_PROP_TYPE(TypeTag, PTAG(SolutionVector)) SolutionVector;
189
190
191
192
193

    typedef NewtonConvergenceWriter<TypeTag, GET_PROP_VALUE(TypeTag, PTAG(NewtonWriteConvergence))>  ConvergenceWriter;

public:
    NewtonController()
Andreas Lauser's avatar
Andreas Lauser committed
194
195
        : endIterMsgStream_(std::ostringstream::out),
          convergenceWriter_(asImp_())
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
    {
        numSteps_ = 0;

        // maximum tolerated deflection between two iterations
        setRelTolerance(1e-7);
        setTargetSteps(8);
        setMaxSteps(12);

        curPhysicalness_ = 0;
        maxPhysicalness_ = 0;
    };

    /*!
     * \brief Set the relative deflection of any degree of freedom
     *        between two iterations which is tolerated.
     */
    void setRelTolerance(Scalar tolerance)
    { tolerance_ = tolerance; }

    /*!
     * \brief Set the number of iterations at which the Newton method
     *        should aim at.
     */
    void setTargetSteps(int targetSteps)
    { targetSteps_ = targetSteps; }

    /*!
     * \brief Set the number of iterations after which the Newton
     *        method gives up.
     */
    void setMaxSteps(int maxSteps)
    { maxSteps_ = maxSteps; }

    /*!
Andreas Lauser's avatar
Andreas Lauser committed
230
231
     * \brief Returns true if another iteration should be done.
     */
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
    bool newtonProceed(SolutionVector &u)
    {
        if (numSteps_ < 2)
            return true; // we always do at least two iterations
        else if (numSteps_ >= maxSteps_) {
            // we have exceeded the allowed number of steps.  if the
            // relative error was reduced by a factor of at least 5,
            // we proceed even if we are above the maximum number of
            // steps
            return error_*4.0 < lastError_ && !asImp_().newtonConverged();
        }
        else if (asImp_().newtonConverged())
            return false; // we are below the desired tolerance

        Scalar tmp = asImp_().physicalness_(u);
Andreas Lauser's avatar
Andreas Lauser committed
247
        curPhysicalness_ = gridView_().comm().min(tmp);
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
        curPhysicalness_ = std::min(curPhysicalness_, Scalar(1.0));


        // check for the physicalness of the solution
        if (curPhysicalness_ <= 0)
            // not physical enough even for a temporary
            // solution
            return false;
        else if (curPhysicalness_ < ((Scalar) numSteps_)/(maxSteps_ - 1)) {
            // we require that the solution gets more physical
            // with every step and at the last step the
            // solution must be completely physical.
            return false;
        }
        else if (curPhysicalness_ < maxPhysicalness_)
        {
            if (probationCount_ > 1) {
                // an iterative solution was more physical
                // than the current solution and at least 2
                // others.
                return false;
            }
            else {
                // we are physical enough, but some earlier
                // solution was more physical, so we let the
                // solver continue on probation.
                ++probationCount_;
                return true;
            }
        }
        else {
            // everything's fine: the solution is physical
            // enough for the number of iterations we did and
            // it is the most physical so far.
            maxPhysicalness_ = curPhysicalness_;
            probationCount_ = std::max(0, probationCount_ - 1);

            return true; // do another round
        };
    }

    /*!
Andreas Lauser's avatar
Andreas Lauser committed
290
291
292
     * \brief Returns true iff the error of the solution is below the
     *        tolerance.
     */
293
294
295
296
297
298
    bool newtonConverged() const
    {
        return (error_ <= tolerance_) && (curPhysicalness_ >= 1.0);
    }

    /*!
Andreas Lauser's avatar
Andreas Lauser committed
299
300
301
302
     * \brief Called before the newton method is applied to an
     *        non-linear system of equations.
     */
    void newtonBegin(NewtonMethod &method, SolutionVector &u)
303
    {
Andreas Lauser's avatar
Andreas Lauser committed
304
        method_ = &method;
305
306
307
308
309
310
311
312
313
        numSteps_ = 0;
        probationCount_ = 0;
        maxPhysicalness_ = 0;
        curPhysicalness_ = 0;

        convergenceWriter_.beginTimestep();
    }

    /*!
Andreas Lauser's avatar
Andreas Lauser committed
314
315
     * \brief Indidicates the beginning of a newton iteration.
     */
316
317
318
319
    void newtonBeginStep()
    { lastError_ = error_; }

    /*!
Andreas Lauser's avatar
Andreas Lauser committed
320
321
322
     * \brief Returns the number of steps done since newtonBegin() was
     *        called.
     */
323
324
325
326
327
    int newtonNumSteps()
    { return numSteps_; }


    /*!
Andreas Lauser's avatar
Andreas Lauser committed
328
329
330
     * \brief Update the error of the solution compared to the
     *        previous iteration.
     */
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
    void newtonUpdateRelError(const SolutionVector &uOld,
                              const SolutionVector &deltaU)
    {
        // calculate the relative error as the maximum relative
        // deflection in any degree of freedom.
        typedef typename SolutionVector::block_type FV;
        error_ = 0;

        int idxI = -1;
        int idxJ = -1;
        for (int i = 0; i < int(uOld.size()); ++i) {
            for (int j = 0; j < FV::size; ++j) {
                Scalar tmp
                    =
                    std::abs(deltaU[i][j])
                    / std::max(std::abs(uOld[i][j]), Scalar(1e-4));
                if (tmp > error_)
                {
Andreas Lauser's avatar
Andreas Lauser committed
349
350
351
                    idxI = i;
                    idxJ = j;
                    error_ = tmp;
352
353
354
                }
            }
        }
Andreas Lauser's avatar
Andreas Lauser committed
355
        error_ = gridView_().comm().max(error_);
356
357
358
359
    }


    /*!
Andreas Lauser's avatar
Andreas Lauser committed
360
361
362
363
364
365
     * \brief Solve the linear system of equations \f$ \mathbf{A}x - b
     *        = 0\f$.
     *
     * Throws Dumux::NumericalProblem if the linear solver didn't
     * converge.
     */
366
367
368
369
370
371
372
373
374
375
376
377
378
    template <class Matrix, class Vector>
    void newtonSolveLinear(Matrix &A,
                           Vector &u,
                           Vector &b)
    {
        // if the deflection of the newton method is large, we do not
        // need to solve the linear approximation accurately. Assuming
        // that the initial value for the delta vector u is quite
        // close to the final value, a reduction of 6 orders of
        // magnitude in the defect should be sufficient...
        Scalar residReduction = 1e-6;

        try {
Andreas Lauser's avatar
Andreas Lauser committed
379
            solveLinear_(A, u, b, residReduction);
380

Andreas Lauser's avatar
Andreas Lauser committed
381
382
383
            // make sure all processes converged
            int converged = 1;
            gridView_().comm().min(converged);
384

Andreas Lauser's avatar
Andreas Lauser committed
385
386
387
388
            if (!converged) {
                DUNE_THROW(NumericalProblem,
                           "A process threw MatrixBlockError");
            }
389
390
391
392
        }
        catch (Dune::MatrixBlockError e) {
            // make sure all processes converged
            int converged = 0;
Andreas Lauser's avatar
Andreas Lauser committed
393
            gridView_().comm().min(converged);
394
395
396
397
398
399
400
401
402
403
404

            Dumux::NumericalProblem p;
            std::string msg;
            std::ostringstream ms(msg);
            ms << e.what() << "M=" << A[e.r][e.c];
            p.message(ms.str());
            throw p;
        }
    }

    /*!
Andreas Lauser's avatar
Andreas Lauser committed
405
406
407
408
409
410
411
412
413
414
415
416
417
418
     * \brief Update the current solution function with a delta vector.
     *
     * The error estimates required for the newtonConverged() and
     * newtonProceed() methods should be updated here.
     *
     * Different update strategies, such as line search and chopped
     * updates can be implemented. The default behaviour is just to
     * subtract deltaU from uOld.
     *
     * \param deltaU The delta as calculated from solving the linear
     *               system of equations. This parameter also stores
     *               the updated solution.
     * \param uOld   The solution of the last iteration
     */
419
420
421
422
423
424
425
426
427
428
429
    void newtonUpdate(SolutionVector &deltaU, const SolutionVector &uOld)
    {
        writeConvergence_(uOld, deltaU);

        newtonUpdateRelError(uOld, deltaU);

        deltaU *= -1;
        deltaU += uOld;
    }

    /*!
Andreas Lauser's avatar
Andreas Lauser committed
430
431
     * \brief Indicates that one newton iteration was finished.
     */
432
433
434
435
436
    void newtonEndStep(SolutionVector &u, SolutionVector &uOld)
    {
        ++numSteps_;

        curPhysicalness_ = asImp_().physicalness_(u);
Andreas Lauser's avatar
Andreas Lauser committed
437
        if (verbose())
438
            std::cout << "\rNewton iteration " << numSteps_ << " done: "
Andreas Lauser's avatar
Andreas Lauser committed
439
                      << "error=" << error_ << endIterMsg().str() << "\n";
440
441
442
443
        endIterMsgStream_.str("");
    }

    /*!
Andreas Lauser's avatar
Andreas Lauser committed
444
445
     * \brief Indicates that we're done solving the non-linear system of equations.
     */
446
447
448
449
450
451
    void newtonEnd()
    {
        convergenceWriter_.endTimestep();
    }

    /*!
Andreas Lauser's avatar
Andreas Lauser committed
452
453
454
455
     * \brief Called if the newton method broke down.
     *
     * This method is called _after_ newtonEnd()
     */
456
457
458
459
460
461
    void newtonFail()
    {
        numSteps_ = targetSteps_*2;
    }

    /*!
Andreas Lauser's avatar
Andreas Lauser committed
462
463
464
465
     * \brief Called when the newton method was sucessful.
     *
     * This method is called _after_ newtonEnd()
     */
466
467
468
469
    void newtonSucceed()
    { }

    /*!
Andreas Lauser's avatar
Andreas Lauser committed
470
471
472
473
474
475
     * \brief Suggest a new time stepsize based on the old time step size.
     *
     * The default behaviour is to suggest the old time step size
     * scaled by the ratio between the target iterations and the
     * iterations required to actually solve the last time step.
     */
476
477
478
479
480
481
482
483
484
485
486
487
488
489
    Scalar suggestTimeStepSize(Scalar oldTimeStep) const
    {
        // be agressive reducing the timestep size but
        // conservative when increasing it. the rationale is
        // that we want to avoid failing in the next newton
        // iteration which would require another linearization
        // of the problem.
        if (numSteps_ > targetSteps_) {
            Scalar percent = ((Scalar) numSteps_ - targetSteps_)/targetSteps_;
            return oldTimeStep/(1.0 + percent);
        }
        else {
            /*Scalar percent = (Scalar(1))/targetSteps_;
              return oldTimeStep*(1 + percent);
Andreas Lauser's avatar
Andreas Lauser committed
490
            */
491
492
493
494
495
496
            Scalar percent = ((Scalar) targetSteps_ - numSteps_)/targetSteps_;
            return oldTimeStep*(1.0 + percent/1.2);
        }
    }

    /*!
Andreas Lauser's avatar
Andreas Lauser committed
497
498
499
     * \brief Returns a reference to the current newton method
     *        which is controlled by this controller.
     */
500
501
502
503
    NewtonMethod &method()
    { return *method_; }

    /*!
Andreas Lauser's avatar
Andreas Lauser committed
504
505
506
     * \brief Returns a reference to the current newton method
     *        which is controlled by this controller.
     */
507
508
509
    const NewtonMethod &method() const
    { return *method_; }

Andreas Lauser's avatar
Andreas Lauser committed
510
511
512
    std::ostringstream &endIterMsg()
    { return endIterMsgStream_; }

513
    /*!
Andreas Lauser's avatar
Andreas Lauser committed
514
515
516
517
     * \brief Returns true iff the newton method ought to be chatty.
     */
    bool verbose() const 
    { return gridView_().comm().rank() == 0; }
518

Andreas Lauser's avatar
Andreas Lauser committed
519
protected:
520
    /*!
Andreas Lauser's avatar
Andreas Lauser committed
521
522
523
524
     * \brief Returns a reference to the grid view.
     */
    const GridView &gridView_() const
    { return problem_().gridView(); }
525

Andreas Lauser's avatar
Andreas Lauser committed
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
    /*!
     * \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 problem.
     */
    Model &model_()
    { return problem_().model(); }

    /*!
     * \brief Returns a reference to the problem.
     */
    const Model &model_() const
    { return problem_().model(); }
549
550
551
552
553
554
555
556
557
558
559
560
561

    // returns the actual implementation for the cotroller we do
    // it this way in order to allow "poor man's virtual methods",
    // i.e. methods of subclasses which can be called by the base
    // class.
    Implementation &asImp_()
    { return *static_cast<Implementation*>(this); }
    const Implementation &asImp_() const
    { return *static_cast<const Implementation*>(this); }

    void writeConvergence_(const SolutionVector &uOld,
                           const SolutionVector &deltaU)
    {
Andreas Lauser's avatar
Andreas Lauser committed
562
        convergenceWriter_.beginIteration(this->gridView_());
563
564
565
566
567
568
569
570
571
572
573
574
575
576
        convergenceWriter_.writeFields(uOld, deltaU);
        convergenceWriter_.endIteration();
    };


    template <class Matrix, class Vector>
    void solveLinear_(Matrix &A,
                      Vector &u,
                      Vector &b,
                      Scalar residReduction)
    {
        Vector &x = u;

        int verbosity = GET_PROP_VALUE(TypeTag, PTAG(NewtonLinearSolverVerbosity));
Andreas Lauser's avatar
Andreas Lauser committed
577
        if (gridView_().comm().rank() != 0)
578
579
580
            verbosity = 0;

#if HAVE_PARDISO
Andreas Lauser's avatar
Andreas Lauser committed
581
582
        typedef  Dumux::PDELab::ISTLBackend_NoOverlap_Loop_Pardiso<TypeTag> Solver;
        Solver solver(problem_(), 500, verbosity);
583
584
585
#else // !HAVE_PARDISO
#if HAVE_MPI
//        typedef  Dune::PDELab::ISTLBackend_NOVLP_BCGS_NOPREC<GridFunctionSpace> Solver;
Andreas Lauser's avatar
Andreas Lauser committed
586
587
588
//        Solver solver(model_().jacobianAssembler().gridFunctionSpace(), 50000, verbosity);
        typedef  Dumux::PDELab::ISTLBackend_NoOverlap_BCGS_ILU<TypeTag> Solver;
        Solver solver(problem_(), 500, verbosity);
589
#else
Andreas Lauser's avatar
Andreas Lauser committed
590
        typedef  Dumux::PDELab::ISTLBackend_SEQ_BCGS_SSOR Solver;
591
592
593
594
        Solver solver(500, verbosity);
#endif // HAVE_MPI
#endif // HAVE_PARDISO

Andreas Lauser's avatar
Andreas Lauser committed
595
        //    Solver solver(model_().jacobianAssembler().gridFunctionSpace(), 500, verbosity);
596
597
598
599
600
601
602
603
604
605
        solver.apply(A, x, b, residReduction);

        if (!solver.result().converged)
            DUNE_THROW(Dumux::NumericalProblem,
                       "Solving the linear system of equations did not converge.");

        // make sure the solver didn't produce a nan or an inf
        // somewhere. this should never happen but for some strange
        // reason it happens anyway.
        Scalar xNorm2 = x.two_norm2();
Andreas Lauser's avatar
Andreas Lauser committed
606
        gridView_().comm().sum(xNorm2);
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
        if (std::isnan(xNorm2) || !std::isfinite(xNorm2))
            DUNE_THROW(Dumux::NumericalProblem,
                       "The linear solver produced a NaN or inf somewhere.");
    }

    //! this function is an indication of how "physically
    //! meaningful" a temporary solution is. 0 means it isn't
    //! meaningful at all (maybe because it has highly negative
    //! pressures, etc) and the newton method can be stopped
    //! immediately. Conversly 1 means that the solution is
    //! perfectly physically meaningful (although it doesn't need
    //! to be the solution in any way) and the method can to run.
    //! Values inbetween mean that the funtion is not meaningfull,
    //! but can be tolerated as temporary solution at some
    //! iteration. (The controller assumes that as the method
    //! progresses, the physicallness of the solution must
    //! increase.)
    Scalar physicalness_(SolutionVector &u)
    {
        return 1;
    }

    std::ostringstream endIterMsgStream_;

    NewtonMethod *method_;

    ConvergenceWriter convergenceWriter_;

    Scalar tolerance_;

    Scalar maxPhysicalness_;
    Scalar curPhysicalness_;
    Scalar error_;
    Scalar lastError_;
    int    probationCount_;

    // optimal number of iterations we want to achive
    int    targetSteps_;
    // maximum number of iterations we do before giving up
    int    maxSteps_;
    // actual number of steps done so far
    int    numSteps_;
};
Andreas Lauser's avatar
Andreas Lauser committed
650
} // namespace Dumux
651
652

#endif