newtoncontroller.hh 25.6 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
 *   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                          *
 *                                                                           *
9
 *   This program is free software: you can redistribute it and/or modify    *
Andreas Lauser's avatar
Andreas Lauser committed
10
 *   it under the terms of the GNU General Public License as published by    *
11
12
 *   the Free Software Foundation, either version 2 of the License, or       *
 *   (at your option) any later version.                                     *
Andreas Lauser's avatar
Andreas Lauser committed
13
 *                                                                           *
14
15
16
17
18
19
20
 *   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/>.   *
Andreas Lauser's avatar
Andreas Lauser committed
21
 *****************************************************************************/
22
/*!
Andreas Lauser's avatar
Andreas Lauser committed
23
 * \file
24
 * \brief Reference implementation of a controller class for the Newton solver.
Andreas Lauser's avatar
Andreas Lauser committed
25
 *
26
 * Usually this controller should be sufficient.
Andreas Lauser's avatar
Andreas Lauser committed
27
 */
28
29
30
#ifndef DUMUX_NEWTON_CONTROLLER_HH
#define DUMUX_NEWTON_CONTROLLER_HH

31
#include <dumux/io/vtkmultiwriter.hh>
32
#include <dumux/common/exceptions.hh>
33
#include <dumux/common/math.hh>
34
35
36

#include <dumux/common/pardiso.hh>

Andreas Lauser's avatar
Andreas Lauser committed
37
#include <dumux/io/vtkmultiwriter.hh>
38
39
40

#if HAVE_DUNE_PDELAB

Andreas Lauser's avatar
Andreas Lauser committed
41
#include <dumux/common/pdelabpreconditioner.hh>
Andreas Lauser's avatar
Andreas Lauser committed
42

43
44
45
46
47
48
49
50
51
52
53
#else // ! HAVE_DUNE_PDELAB

#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>

#endif // HAVE_DUNE_PDELAB

54
55
56
57
58

namespace Dumux
{
namespace Properties
{
59
//! Specifies the implementation of the Newton controller
Andreas Lauser's avatar
Andreas Lauser committed
60
61
NEW_PROP_TAG(NewtonController);

62
//! Specifies the type of the actual Newton method
Andreas Lauser's avatar
Andreas Lauser committed
63
64
NEW_PROP_TAG(NewtonMethod);

65
//! Specifies the type of a solution
Andreas Lauser's avatar
Andreas Lauser committed
66
67
NEW_PROP_TAG(SolutionVector);

68
//! Specifies the type of a vector of primary variables at a degree of freedom
Andreas Lauser's avatar
Andreas Lauser committed
69
70
NEW_PROP_TAG(PrimaryVariables);

71
//! Specifies the type of a global Jacobian matrix
Andreas Lauser's avatar
Andreas Lauser committed
72
73
NEW_PROP_TAG(JacobianMatrix);

74
//! Specifies the type of the Jacobian matrix assembler
Andreas Lauser's avatar
Andreas Lauser committed
75
76
77
78
79
NEW_PROP_TAG(JacobianAssembler);

//! specifies the type of the time manager
NEW_PROP_TAG(TimeManager);

80
81
82
83
84
85
86
/*!
 * \brief Specifies the verbosity of the linear solver
 *
 * By default it is 0, i.e. it doesn't print anything. Setting this
 * property to 1 prints aggregated convergence rates, 2 prints the
 * convergence rate of every iteration of the scheme.
 */
87
88
89
90
91
92
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);

93
94
//! Specifies whether time step size should be increased during the
//! Newton methods first few iterations
Andreas Lauser's avatar
Andreas Lauser committed
95
96
NEW_PROP_TAG(EnableTimeStepRampUp);

97
//! Specifies whether the Jacobian matrix should only be reassembled
Andreas Lauser's avatar
Andreas Lauser committed
98
99
100
//! if the current solution deviates too much from the evaluation point
NEW_PROP_TAG(EnablePartialReassemble);

101
102
/*!
 * \brief Specifies whether the update should be done using the line search
103
 *        method instead of the plain Newton method.
104
105
106
107
108
 *
 * 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 line search is not used.
 */
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
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;
};
};

127
//! \cond INTERNAL
128
/*!
129
 * \brief Writes the intermediate solutions during
130
131
 *        the Newton scheme
 */
132
133
134
135
136
137
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
138
    typedef typename GET_PROP_TYPE(TypeTag, PTAG(SolutionVector)) SolutionVector;
139
140
141
    typedef Dumux::VtkMultiWriter<GridView>  VtkMultiWriter;

    NewtonConvergenceWriter(NewtonController &ctl)
Andreas Lauser's avatar
Andreas Lauser committed
142
        : ctl_(ctl)
143
    {
144
        timeStepIndex_ = 0;
145
146
147
148
149
150
151
152
153
        iteration_ = 0;
        vtkMultiWriter_ = new VtkMultiWriter("convergence");
    }

    ~NewtonConvergenceWriter()
    { delete vtkMultiWriter_; };

    void beginTimestep()
    {
154
        ++timeStepIndex_;
155
156
157
158
159
160
        iteration_ = 0;
    };

    void beginIteration(const GridView &gv)
    {
        ++ iteration_;
161
        vtkMultiWriter_->beginTimestep(timeStepIndex_ + iteration_ / 100.0,
Andreas Lauser's avatar
Andreas Lauser committed
162
                                       gv);
163
164
    };

165
    void writeFields(const SolutionVector &uLastIter,
166
167
                     const SolutionVector &deltaU)
    {
168
        ctl_.method().model().addConvergenceVtkFields(*vtkMultiWriter_, uLastIter, deltaU);
169
170
171
    };

    void endIteration()
Andreas Lauser's avatar
Andreas Lauser committed
172
    { vtkMultiWriter_->endTimestep(); };
173
174
175

    void endTimestep()
    {
176
        ++timeStepIndex_;
177
178
179
180
        iteration_ = 0;
    };

private:
181
    int timeStepIndex_;
182
183
184
185
186
    int iteration_;
    VtkMultiWriter *vtkMultiWriter_;
    NewtonController &ctl_;
};

187
/*!
188
189
 * \brief Writes the intermediate solutions during
 *        the Newton scheme.
190
191
192
193
 *
 * This is the dummy specialization for the case where we don't want
 * to do anything.
 */
194
195
196
197
198
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
199
    typedef typename GET_PROP_TYPE(TypeTag, PTAG(SolutionVector)) SolutionVector;
200
201
202
203
204
205
206
207
208
209
210
211

    typedef Dumux::VtkMultiWriter<GridView>  VtkMultiWriter;

    NewtonConvergenceWriter(NewtonController &ctl)
    {};

    void beginTimestep()
    { };

    void beginIteration(const GridView &gv)
    { };

212
    void writeFields(const SolutionVector &uLastIter,
Andreas Lauser's avatar
Andreas Lauser committed
213
                     const SolutionVector &deltaU)
214
215
216
217
218
219
220
221
    { };

    void endIteration()
    { };

    void endTimestep()
    { };
};
222
//! \endcond
223
224

/*!
225
226
 * \brief A reference implementation of a newton controller specific
 *        for the box scheme.
Andreas Lauser's avatar
Andreas Lauser committed
227
 *
228
229
230
 * 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.
Andreas Lauser's avatar
Andreas Lauser committed
231
 */
232
233
234
template <class TypeTag>
class NewtonController
{
Andreas Lauser's avatar
Andreas Lauser committed
235
236
237
    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;
238

Andreas Lauser's avatar
Andreas Lauser committed
239
240
    typedef typename GET_PROP_TYPE(TypeTag, PTAG(Problem)) Problem;
    typedef typename GET_PROP_TYPE(TypeTag, PTAG(Model)) Model;
241
    typedef typename GET_PROP_TYPE(TypeTag, PTAG(NewtonMethod)) NewtonMethod;
Bernd Flemisch's avatar
Bernd Flemisch committed
242
    typedef typename GET_PROP_TYPE(TypeTag, PTAG(JacobianMatrix)) JacobianMatrix;
243
    typedef typename GET_PROP_TYPE(TypeTag, PTAG(TimeManager)) TimeManager;
244

Andreas Lauser's avatar
Andreas Lauser committed
245
    typedef typename GET_PROP_TYPE(TypeTag, PTAG(JacobianAssembler)) JacobianAssembler;
Andreas Lauser's avatar
Andreas Lauser committed
246
    typedef typename GET_PROP_TYPE(TypeTag, PTAG(SolutionVector)) SolutionVector;
Andreas Lauser's avatar
Andreas Lauser committed
247
    typedef typename GET_PROP_TYPE(TypeTag, PTAG(PrimaryVariables)) PrimaryVariables;
248
249
250

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

251
    enum { enableTimeStepRampUp = GET_PROP_VALUE(TypeTag, PTAG(EnableTimeStepRampUp)) };
Andreas Lauser's avatar
Andreas Lauser committed
252
253
    enum { enablePartialReassemble = GET_PROP_VALUE(TypeTag, PTAG(EnablePartialReassemble)) };

254
public:
255
256
257
    /*!
     * \brief Constructor
     */
258
    NewtonController()
Andreas Lauser's avatar
Andreas Lauser committed
259
260
        : endIterMsgStream_(std::ostringstream::out),
          convergenceWriter_(asImp_())
261
    {
262
        verbose_ = true;
263
264
        numSteps_ = 0;

Andreas Lauser's avatar
Andreas Lauser committed
265
266
267
268
269
        this->setRelTolerance(1e-8);
        this->rampUpSteps_ = 0;

        if (enableTimeStepRampUp) {
            this->rampUpSteps_ = 9;
270

Andreas Lauser's avatar
Andreas Lauser committed
271
272
273
274
275
276
277
278
            // the ramp-up steps are not counting
            this->setTargetSteps(10);
            this->setMaxSteps(12);
        }
        else {
            this->setTargetSteps(10);
            this->setMaxSteps(18);
        }
279
280
281
    };

    /*!
282
     * \brief Set the maximum acceptable difference for convergence of
283
284
285
286
287
     *        any primary variable between two iterations.
     *
     * \param tolerance The maximum relative error between two Newton
     *                  iterations at which the scheme is considered
     *                  finished
288
289
290
291
292
293
294
     */
    void setRelTolerance(Scalar tolerance)
    { tolerance_ = tolerance; }

    /*!
     * \brief Set the number of iterations at which the Newton method
     *        should aim at.
295
296
297
298
299
300
     *
     * This is used to control the time step size. The heuristic used
     * is to scale the last time step size by the deviation of the
     * number of iterations used from the target steps.
     *
     * \param targetSteps Number of iterations which are considered "optimal"
301
302
303
304
305
306
307
     */
    void setTargetSteps(int targetSteps)
    { targetSteps_ = targetSteps; }

    /*!
     * \brief Set the number of iterations after which the Newton
     *        method gives up.
308
309
     *
     * \param maxSteps Number of iterations after we give up
310
311
312
     */
    void setMaxSteps(int maxSteps)
    { maxSteps_ = maxSteps; }
313

Andreas Lauser's avatar
Andreas Lauser committed
314
315
316
317
318
319
320
321
322
323
324
325
326
    /*!
     * \brief Returns the number of iterations used for the time step
     *        ramp-up.
     */
    Scalar rampUpSteps() const
    { return enableTimeStepRampUp?rampUpSteps_:0; }

    /*!
     * \brief Returns whether the time-step ramp-up is still happening
     */
    bool inRampUp() const
    { return numSteps_ < rampUpSteps(); }

327
    /*!
Andreas Lauser's avatar
Andreas Lauser committed
328
     * \brief Returns true if another iteration should be done.
329
     *
330
     * \param uCurrentIter The solution of the current newton iteration
Andreas Lauser's avatar
Andreas Lauser committed
331
     */
332
    bool newtonProceed(const SolutionVector &uCurrentIter)
333
    {
Andreas Lauser's avatar
Andreas Lauser committed
334
        if (numSteps_ < rampUpSteps() + 2)
335
            return true; // we always do at least two iterations
Andreas Lauser's avatar
Andreas Lauser committed
336
337
338
        else if (asImp_().newtonConverged())
            return false; // we are below the desired tolerance
        else if (numSteps_ >= rampUpSteps() + maxSteps_) {
339
            // we have exceeded the allowed number of steps.  if the
340
            // relative error was reduced by a factor of at least 4,
341
342
            // we proceed even if we are above the maximum number of
            // steps
Andreas Lauser's avatar
Andreas Lauser committed
343
            return error_*4.0 < lastError_;
344
345
        }

346
        return true;
347
348
349
    }

    /*!
Andreas Lauser's avatar
Andreas Lauser committed
350
351
352
     * \brief Returns true iff the error of the solution is below the
     *        tolerance.
     */
353
354
    bool newtonConverged() const
    {
355
        return error_ <= tolerance_;
356
357
358
    }

    /*!
Andreas Lauser's avatar
Andreas Lauser committed
359
360
     * \brief Called before the newton method is applied to an
     *        non-linear system of equations.
361
362
363
     *
     * \param method The object where the NewtonMethod is executed
     * \param u The initial solution
Andreas Lauser's avatar
Andreas Lauser committed
364
     */
365
    void newtonBegin(NewtonMethod &method, const SolutionVector &u)
366
    {
Andreas Lauser's avatar
Andreas Lauser committed
367
        method_ = &method;
368
369
        numSteps_ = 0;

Andreas Lauser's avatar
Andreas Lauser committed
370
371
372
        model_().jacobianAssembler().reassembleAll();
        dtInitial_ = timeManager_().timeStepSize();
        if (enableTimeStepRampUp) {
373
374
            rampUpDelta_ =
                timeManager_().timeStepSize()
Andreas Lauser's avatar
Andreas Lauser committed
375
376
377
378
379
380
381
                /
                rampUpSteps()
                *
                2;

            // reduce initial time step size for ramp-up.
            timeManager_().setTimeStepSize(rampUpDelta_);
382
383
        }

384
385
386
387
        convergenceWriter_.beginTimestep();
    }

    /*!
Andreas Lauser's avatar
Andreas Lauser committed
388
389
     * \brief Indidicates the beginning of a newton iteration.
     */
390
391
392
393
    void newtonBeginStep()
    { lastError_ = error_; }

    /*!
Andreas Lauser's avatar
Andreas Lauser committed
394
395
396
     * \brief Returns the number of steps done since newtonBegin() was
     *        called.
     */
397
398
399
400
    int newtonNumSteps()
    { return numSteps_; }

    /*!
401
402
403
404
405
406
     * \brief Update the relative error of the solution compared to
     *        the previous iteration.
     *
     * The relative error can be seen as a norm of the difference
     * between the current and the next iteration.
     *
407
     * \param uLastIter The current iterative solution
408
     * \param deltaU The difference between the current and the next solution
Andreas Lauser's avatar
Andreas Lauser committed
409
     */
410
    void newtonUpdateRelError(const SolutionVector &uLastIter,
411
412
413
414
415
416
417
418
                              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;
Andreas Lauser's avatar
Andreas Lauser committed
419
        int aboveTol = 0;
420
421
        for (int i = 0; i < int(uLastIter.size()); ++i) {
            PrimaryVariables uNewI = uLastIter[i];
Andreas Lauser's avatar
Andreas Lauser committed
422
            uNewI -= deltaU[i];
423
            Scalar vertErr =
Andreas Lauser's avatar
Andreas Lauser committed
424
                model_().relativeErrorVertex(i,
425
                                             uLastIter[i],
Andreas Lauser's avatar
Andreas Lauser committed
426
                                             uNewI);
427

Andreas Lauser's avatar
Andreas Lauser committed
428
429
430
431
432
433
            if (vertErr > tolerance_)
                ++aboveTol;
            if (vertErr > error_) {
                idxI = i;
                error_ = vertErr;
            }
434
        }
435

Bernd Flemisch's avatar
Bernd Flemisch committed
436
437
        if (gridView_().comm().size() > 1)
        	error_ = gridView_().comm().max(error_);
438
439
440
    }

    /*!
441
     * \brief Solve the linear system of equations \f$\mathbf{A}x - b = 0\f$.
Andreas Lauser's avatar
Andreas Lauser committed
442
443
444
     *
     * Throws Dumux::NumericalProblem if the linear solver didn't
     * converge.
445
446
     *
     * \param A The matrix of the linear system of equations
447
     * \param x The vector which solves the linear system
448
     * \param b The right hand side of the linear system
Andreas Lauser's avatar
Andreas Lauser committed
449
     */
Bernd Flemisch's avatar
Bernd Flemisch committed
450
451
    template <class Vector>
    void newtonSolveLinear(const JacobianMatrix &A,
452
                           Vector &x,
453
                           const Vector &b)
454
455
456
457
458
459
460
461
462
    {
        // 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 {
463
            solveLinear_(A, x, b, residReduction);
464

Andreas Lauser's avatar
Andreas Lauser committed
465
466
            // make sure all processes converged
            int converged = 1;
Bernd Flemisch's avatar
Bernd Flemisch committed
467
468
            if (gridView_().comm().size() > 1)
            	gridView_().comm().min(converged);
469

Andreas Lauser's avatar
Andreas Lauser committed
470
471
            if (!converged) {
                DUNE_THROW(NumericalProblem,
472
                           "A process threw NumericalProblem");
Andreas Lauser's avatar
Andreas Lauser committed
473
            }
474
475
476
477
        }
        catch (Dune::MatrixBlockError e) {
            // make sure all processes converged
            int converged = 0;
Bernd Flemisch's avatar
Bernd Flemisch committed
478
479
            if (gridView_().comm().size() > 1)
            	gridView_().comm().min(converged);
480
481
482
483
484
485
486
487

            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;
        }
488
489
490
        catch (const Dune::ISTLError &e) {
            // make sure all processes converged
            int converged = 0;
Bernd Flemisch's avatar
Bernd Flemisch committed
491
492
            if (gridView_().comm().size() > 1)
            	gridView_().comm().min(converged);
493
494
495
496
497

            Dumux::NumericalProblem p;
            p.message(e.what());
            throw p;
        }
498
499
500
    }

    /*!
501
     * \brief Update the current solution with a delta vector.
Andreas Lauser's avatar
Andreas Lauser committed
502
503
     *
     * The error estimates required for the newtonConverged() and
504
     * newtonProceed() methods should be updated inside this method.
Andreas Lauser's avatar
Andreas Lauser committed
505
506
507
     *
     * Different update strategies, such as line search and chopped
     * updates can be implemented. The default behaviour is just to
508
     * subtract deltaU from uLastIter, i.e.
509
     * \f[ u^{k+1} = u^k - \Delta u^k \f]
Andreas Lauser's avatar
Andreas Lauser committed
510
     *
511
512
     * \param uCurrentIter The solution vector after the current iteration
     * \param uLastIter The solution vector after the last iteration
Andreas Lauser's avatar
Andreas Lauser committed
513
514
515
516
     * \param deltaU The delta as calculated from solving the linear
     *               system of equations. This parameter also stores
     *               the updated solution.
     */
517
518
519
    void newtonUpdate(SolutionVector &uCurrentIter,
                      const SolutionVector &uLastIter,
                      const SolutionVector &deltaU)
520
    {
521
        writeConvergence_(uLastIter, deltaU);
522

523
        newtonUpdateRelError(uLastIter, deltaU);
Andreas Lauser's avatar
Andreas Lauser committed
524
525
526
527

        // compute the vertex and element colors for partial
        // reassembly
        if (enablePartialReassemble) {
528
529
            Scalar reassembleTol = Dumux::geometricMean(error_, 0.1*tolerance_);
            reassembleTol = std::max(reassembleTol, 0.1*tolerance_);
530
            this->model_().jacobianAssembler().updateDiscrepancy(uLastIter, deltaU);
Andreas Lauser's avatar
Andreas Lauser committed
531
            this->model_().jacobianAssembler().computeColors(reassembleTol);
532
        }
533

534
535
        uCurrentIter = uLastIter;
        uCurrentIter -= deltaU;
Andreas Lauser's avatar
Andreas Lauser committed
536
    }
537

538
    /*!
Andreas Lauser's avatar
Andreas Lauser committed
539
     * \brief Indicates that one newton iteration was finished.
540
     *
541
542
     * \param uCurrentIter The solution after the current Newton iteration
     * \param uLastIter The solution at the beginning of the current Newton iteration
Andreas Lauser's avatar
Andreas Lauser committed
543
     */
544
    void newtonEndStep(const SolutionVector &uCurrentIter,
545
                       const SolutionVector &uLastIter)
546
547
    {
        ++numSteps_;
Andreas Lauser's avatar
Andreas Lauser committed
548
549
550
551
552
553
554
555
556

        Scalar realError = error_;
        if (inRampUp() && error_ < 1.0) {
            // change time step size
            Scalar dt = timeManager_().timeStepSize();
            dt += rampUpDelta_;
            timeManager_().setTimeStepSize(dt);

            endIterMsg() << ", dt=" << timeManager_().timeStepSize() << ", ddt=" << rampUpDelta_;
557
        }
558

Andreas Lauser's avatar
Andreas Lauser committed
559
        if (verbose())
560
            std::cout << "\rNewton iteration " << numSteps_ << " done: "
Andreas Lauser's avatar
Andreas Lauser committed
561
                      << "error=" << realError << endIterMsg().str() << "\n";
562
563
564
565
        endIterMsgStream_.str("");
    }

    /*!
566
567
     * \brief Indicates that we're done solving the non-linear system
     *        of equations.
Andreas Lauser's avatar
Andreas Lauser committed
568
     */
569
570
571
572
573
574
    void newtonEnd()
    {
        convergenceWriter_.endTimestep();
    }

    /*!
Andreas Lauser's avatar
Andreas Lauser committed
575
576
577
578
     * \brief Called if the newton method broke down.
     *
     * This method is called _after_ newtonEnd()
     */
579
580
    void newtonFail()
    {
Andreas Lauser's avatar
Andreas Lauser committed
581
        timeManager_().setTimeStepSize(dtInitial_);
582
583
584
585
        numSteps_ = targetSteps_*2;
    }

    /*!
Andreas Lauser's avatar
Andreas Lauser committed
586
587
588
589
     * \brief Called when the newton method was sucessful.
     *
     * This method is called _after_ newtonEnd()
     */
590
591
592
593
    void newtonSucceed()
    { }

    /*!
594
595
     * \brief Suggest a new time stepsize based on the old time step
     *        size.
Andreas Lauser's avatar
Andreas Lauser committed
596
597
598
599
600
     *
     * 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.
     */
601
602
    Scalar suggestTimeStepSize(Scalar oldTimeStep) const
    {
Andreas Lauser's avatar
Andreas Lauser committed
603
        if (enableTimeStepRampUp)
604
            return oldTimeStep;
Andreas Lauser's avatar
Andreas Lauser committed
605

606
        Scalar n = numSteps_;
Andreas Lauser's avatar
Andreas Lauser committed
607
608
        n -= rampUpSteps();

609
610
611
612
613
        // 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.
614
615
        if (n > targetSteps_) {
            Scalar percent = (n - targetSteps_)/targetSteps_;
616
617
618
619
620
            return oldTimeStep/(1.0 + percent);
        }
        else {
            /*Scalar percent = (Scalar(1))/targetSteps_;
              return oldTimeStep*(1 + percent);
Andreas Lauser's avatar
Andreas Lauser committed
621
            */
622
            Scalar percent = (targetSteps_ - n)/targetSteps_;
623
624
625
626
627
            return oldTimeStep*(1.0 + percent/1.2);
        }
    }

    /*!
Andreas Lauser's avatar
Andreas Lauser committed
628
629
630
     * \brief Returns a reference to the current newton method
     *        which is controlled by this controller.
     */
631
632
633
634
    NewtonMethod &method()
    { return *method_; }

    /*!
Andreas Lauser's avatar
Andreas Lauser committed
635
636
637
     * \brief Returns a reference to the current newton method
     *        which is controlled by this controller.
     */
638
639
640
    const NewtonMethod &method() const
    { return *method_; }

Andreas Lauser's avatar
Andreas Lauser committed
641
642
643
    std::ostringstream &endIterMsg()
    { return endIterMsgStream_; }

644
645
646
647
648
649
    /*!
     * \brief Specifies if the newton method ought to be chatty.
     */
    void setVerbose(bool val)
    { verbose_ = val; }

650
    /*!
Andreas Lauser's avatar
Andreas Lauser committed
651
652
     * \brief Returns true iff the newton method ought to be chatty.
     */
653
    bool verbose() const
654
    { return verbose_ && gridView_().comm().rank() == 0; }
655

Andreas Lauser's avatar
Andreas Lauser committed
656
protected:
657
    /*!
Andreas Lauser's avatar
Andreas Lauser committed
658
659
660
661
     * \brief Returns a reference to the grid view.
     */
    const GridView &gridView_() const
    { return problem_().gridView(); }
662

Andreas Lauser's avatar
Andreas Lauser committed
663
664
665
666
667
668
669
670
671
672
673
674
    /*!
     * \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(); }

675
676
677
678
679
680
    /*!
     * \brief Returns a reference to the time manager.
     */
    TimeManager &timeManager_()
    { return problem_().timeManager(); }

Andreas Lauser's avatar
Andreas Lauser committed
681
682
683
684
685
686
    /*!
     * \brief Returns a reference to the time manager.
     */
    const TimeManager &timeManager_() const
    { return problem_().timeManager(); }

Andreas Lauser's avatar
Andreas Lauser committed
687
688
689
690
691
692
693
694
695
696
697
    /*!
     * \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(); }
698
699
700
701
702
703
704
705
706
707

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

708
    void writeConvergence_(const SolutionVector &uLastIter,
709
710
                           const SolutionVector &deltaU)
    {
Andreas Lauser's avatar
Andreas Lauser committed
711
        convergenceWriter_.beginIteration(this->gridView_());
712
        convergenceWriter_.writeFields(uLastIter, deltaU);
713
714
715
        convergenceWriter_.endIteration();
    };

716
    /*!
717
     * \brief Actually invoke the linear solver
718
719
720
     *
     * Usually we use the solvers from DUNE-ISTL.
     */
Bernd Flemisch's avatar
Bernd Flemisch committed
721
722
    template <class Vector>
    void solveLinear_(const JacobianMatrix &A,
723
724
                      Vector &x,
                      const Vector &b,
725
726
727
                      Scalar residReduction)
    {
        int verbosity = GET_PROP_VALUE(TypeTag, PTAG(NewtonLinearSolverVerbosity));
Andreas Lauser's avatar
Andreas Lauser committed
728
        if (gridView_().comm().rank() != 0)
729
730
            verbosity = 0;

731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
#if ! HAVE_DUNE_PDELAB
       typedef Dune::SeqILU0<JacobianMatrix, Vector, Vector> Preconditioner;
       Preconditioner precond(A, 1.0);

       typedef Dune::MatrixAdapter<JacobianMatrix,Vector,Vector> MatrixAdapter;
       MatrixAdapter operatorA(A);

       typedef Dune::BiCGSTABSolver<Vector> Solver;
       Solver solver(operatorA, precond, residReduction, 500, verbosity);
//        typedef Dune::RestartedGMResSolver<Vector> Solver;
//        Solver solver(operatorA, precond, residReduction, 50, 500, verbosity);

       Dune::InverseOperatorResult result;

        Vector bTmp(b);
       solver.apply(x, bTmp, result);

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

753
#if HAVE_PARDISO
754
        typedef Dumux::PDELab::ISTLBackend_NoOverlap_Loop_Pardiso<TypeTag> Solver;
Andreas Lauser's avatar
Andreas Lauser committed
755
        Solver solver(problem_(), 500, verbosity);
756
757
#else // !HAVE_PARDISO
#if HAVE_MPI
758
//        typedef Dune::PDELab::ISTLBackend_NOVLP_BCGS_NOPREC<GridFunctionSpace> Solver;
Andreas Lauser's avatar
Andreas Lauser committed
759
//        Solver solver(model_().jacobianAssembler().gridFunctionSpace(), 50000, verbosity);
760
        typedef Dumux::PDELab::ISTLBackend_NoOverlap_BCGS_ILU<TypeTag> Solver;
Andreas Lauser's avatar
Andreas Lauser committed
761
        Solver solver(problem_(), 500, verbosity);
762
#else
Andreas Lauser's avatar
Andreas Lauser committed
763
        typedef Dune::PDELab::ISTLBackend_SEQ_BCGS_SSOR Solver;
764
765
766
767
        Solver solver(500, verbosity);
#endif // HAVE_MPI
#endif // HAVE_PARDISO

Andreas Lauser's avatar
Andreas Lauser committed
768
        //    Solver solver(model_().jacobianAssembler().gridFunctionSpace(), 500, verbosity);
769
770
        Vector bTmp(b);
        solver.apply(A, x, bTmp, residReduction);
771
772
773
774

        if (!solver.result().converged)
            DUNE_THROW(Dumux::NumericalProblem,
                       "Solving the linear system of equations did not converge.");
775
#endif // HAVE_DUNE_PDELAB
776
777
778
779
780

        // 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();
Bernd Flemisch's avatar
Bernd Flemisch committed
781
782
        if (gridView_().comm().size() > 1)
        	gridView_().comm().sum(xNorm2);
783
784
785
786
787
788
789
790
        if (std::isnan(xNorm2) || !std::isfinite(xNorm2))
            DUNE_THROW(Dumux::NumericalProblem,
                       "The linear solver produced a NaN or inf somewhere.");
    }

    std::ostringstream endIterMsgStream_;

    NewtonMethod *method_;
791
    bool verbose_;
792
793
794
795
796

    ConvergenceWriter convergenceWriter_;

    Scalar error_;
    Scalar lastError_;
797
    Scalar tolerance_;
798

799
800
    // number of iterations for the time-step ramp-up
    Scalar rampUpSteps_;
Andreas Lauser's avatar
Andreas Lauser committed
801
802
803
804
    // the increase of the time step size during the rampup
    Scalar rampUpDelta_;

    Scalar dtInitial_; // initial time step size
805

806
    // optimal number of iterations we want to achive
807
    int targetSteps_;
808
    // maximum number of iterations we do before giving up
809
    int maxSteps_;
810
    // actual number of steps done so far
811
    int numSteps_;
812
};
Andreas Lauser's avatar
Andreas Lauser committed
813
} // namespace Dumux
814
815

#endif