newtoncontroller.hh 22.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
#ifndef DUMUX_NEWTON_CONTROLLER_HH
#define DUMUX_NEWTON_CONTROLLER_HH

#include <dumux/common/exceptions.hh>

Andreas Lauser's avatar
Andreas Lauser committed
28
#include <queue> // for std::priority_queue
29
30
31

#include <dumux/common/pardiso.hh>

Andreas Lauser's avatar
Andreas Lauser committed
32
#include <dumux/io/vtkmultiwriter.hh>
Andreas Lauser's avatar
Andreas Lauser committed
33
#include <dumux/common/pdelabpreconditioner.hh>
Andreas Lauser's avatar
Andreas Lauser committed
34

35
36
37
38
39
40


namespace Dumux
{
namespace Properties
{
Andreas Lauser's avatar
Andreas Lauser committed
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
//! specifies the implementation of the newton controller
NEW_PROP_TAG(NewtonController);

//! specifies the type of the actual newton method
NEW_PROP_TAG(NewtonMethod);

//! specifies the type of a solution
NEW_PROP_TAG(SolutionVector);

//! specifies the type of a vector of primary variables at an DOF
NEW_PROP_TAG(PrimaryVariables);

//! specifies the type of a global jacobian matrix
NEW_PROP_TAG(JacobianMatrix);

//! specifies the type of the jacobian matrix assembler
NEW_PROP_TAG(JacobianAssembler);

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

62
63
64
65
66
67
68
69
//! 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);

Andreas Lauser's avatar
Andreas Lauser committed
70
71
72
73
74
75
76
77
//! specifies whether time step size should be increased during the
//! newton methods first few iterations
NEW_PROP_TAG(EnableTimeStepRampUp);

//! specifies whether the jacobian matrix should only be reassembled
//! if the current solution deviates too much from the evaluation point
NEW_PROP_TAG(EnablePartialReassemble);

78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
//! 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
108
    typedef typename GET_PROP_TYPE(TypeTag, PTAG(SolutionVector)) SolutionVector;
109
110
111
    typedef Dumux::VtkMultiWriter<GridView>  VtkMultiWriter;

    NewtonConvergenceWriter(NewtonController &ctl)
Andreas Lauser's avatar
Andreas Lauser committed
112
        : ctl_(ctl)
113
    {
114
        timeStepIndex_ = 0;
115
116
117
118
119
120
121
122
123
        iteration_ = 0;
        vtkMultiWriter_ = new VtkMultiWriter("convergence");
    }

    ~NewtonConvergenceWriter()
    { delete vtkMultiWriter_; };

    void beginTimestep()
    {
124
        ++timeStepIndex_;
125
126
127
128
129
130
        iteration_ = 0;
    };

    void beginIteration(const GridView &gv)
    {
        ++ iteration_;
131
        vtkMultiWriter_->beginTimestep(timeStepIndex_ + iteration_ / 100.0,
Andreas Lauser's avatar
Andreas Lauser committed
132
                                       gv);
133
134
135
136
137
    };

    void writeFields(const SolutionVector &uOld,
                     const SolutionVector &deltaU)
    {
Andreas Lauser's avatar
Andreas Lauser committed
138
        ctl_.method().model().addConvergenceVtkFields(*vtkMultiWriter_, uOld, deltaU);
139
140
141
    };

    void endIteration()
Andreas Lauser's avatar
Andreas Lauser committed
142
    { vtkMultiWriter_->endTimestep(); };
143
144
145

    void endTimestep()
    {
146
        ++timeStepIndex_;
147
148
149
150
        iteration_ = 0;
    };

private:
151
    int timeStepIndex_;
152
153
154
155
156
157
158
159
160
161
    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
162
    typedef typename GET_PROP_TYPE(TypeTag, PTAG(SolutionVector)) SolutionVector;
163
164
165
166
167
168
169
170
171
172
173
174
175

    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
176
                     const SolutionVector &deltaU)
177
178
179
180
181
182
183
184
185
186
    { };

    void endIteration()
    { };

    void endTimestep()
    { };
};

/*!
Andreas Lauser's avatar
Andreas Lauser committed
187
188
189
190
191
192
193
 * \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.
 */
194
195
196
template <class TypeTag>
class NewtonController
{
Andreas Lauser's avatar
Andreas Lauser committed
197
198
199
    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;
200

Andreas Lauser's avatar
Andreas Lauser committed
201
202
    typedef typename GET_PROP_TYPE(TypeTag, PTAG(Problem)) Problem;
    typedef typename GET_PROP_TYPE(TypeTag, PTAG(Model)) Model;
203
    typedef typename GET_PROP_TYPE(TypeTag, PTAG(NewtonMethod)) NewtonMethod;
Bernd Flemisch's avatar
Bernd Flemisch committed
204
    typedef typename GET_PROP_TYPE(TypeTag, PTAG(JacobianMatrix)) JacobianMatrix;
205
    typedef typename GET_PROP_TYPE(TypeTag, PTAG(TimeManager)) TimeManager;
206

Andreas Lauser's avatar
Andreas Lauser committed
207
    typedef typename GET_PROP_TYPE(TypeTag, PTAG(JacobianAssembler)) JacobianAssembler;
Andreas Lauser's avatar
Andreas Lauser committed
208
    typedef typename GET_PROP_TYPE(TypeTag, PTAG(SolutionVector)) SolutionVector;
Andreas Lauser's avatar
Andreas Lauser committed
209
    typedef typename GET_PROP_TYPE(TypeTag, PTAG(PrimaryVariables)) PrimaryVariables;
210
211
212

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

213
    enum { enableTimeStepRampUp = GET_PROP_VALUE(TypeTag, PTAG(EnableTimeStepRampUp)) };
Andreas Lauser's avatar
Andreas Lauser committed
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
    enum { enablePartialReassemble = GET_PROP_VALUE(TypeTag, PTAG(EnablePartialReassemble)) };
    enum { Red = JacobianAssembler::Red };

    // class to keep track of the most offending vertices in a way
    // compatible with std::priority_queue
    class VertexError
    {
    public:
        VertexError(int idx, Scalar err)
        {
            idx_ = idx;
            err_ = err;
        }
        
        int index() const
        { return idx_; }
        
        bool operator<(const VertexError &a) const
        { return a.err_ < err_; }

    private:
        int idx_;
        Scalar err_;
    };
238

239
240
public:
    NewtonController()
Andreas Lauser's avatar
Andreas Lauser committed
241
242
        : endIterMsgStream_(std::ostringstream::out),
          convergenceWriter_(asImp_())
243
244
245
    {
        numSteps_ = 0;

Andreas Lauser's avatar
Andreas Lauser committed
246
247
248
249
250
251
252
253
254
255
256
257
258
259
        this->setRelTolerance(1e-8);
        this->rampUpSteps_ = 0;

        if (enableTimeStepRampUp) {
            this->rampUpSteps_ = 9;
            
            // the ramp-up steps are not counting
            this->setTargetSteps(10);
            this->setMaxSteps(12);
        }
        else {
            this->setTargetSteps(10);
            this->setMaxSteps(18);
        }
260
261
262
    };

    /*!
263
264
     * \brief Set the maximum acceptable difference for convergence of
     *        any primary variable between two iterations
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
     */
    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; }
282
    
Andreas Lauser's avatar
Andreas Lauser committed
283
284
285
286
287
288
289
290
291
292
293
294
295
    /*!
     * \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(); }

296
    /*!
Andreas Lauser's avatar
Andreas Lauser committed
297
298
     * \brief Returns true if another iteration should be done.
     */
299
    bool newtonProceed(const SolutionVector &u)
300
    {
Andreas Lauser's avatar
Andreas Lauser committed
301
        if (numSteps_ < rampUpSteps() + 2)
302
            return true; // we always do at least two iterations
Andreas Lauser's avatar
Andreas Lauser committed
303
304
305
        else if (asImp_().newtonConverged())
            return false; // we are below the desired tolerance
        else if (numSteps_ >= rampUpSteps() + maxSteps_) {
306
            // we have exceeded the allowed number of steps.  if the
307
            // relative error was reduced by a factor of at least 4,
308
309
            // we proceed even if we are above the maximum number of
            // steps
Andreas Lauser's avatar
Andreas Lauser committed
310
            return error_*4.0 < lastError_;
311
312
        }

313
        return true;
314
315
316
    }

    /*!
Andreas Lauser's avatar
Andreas Lauser committed
317
318
319
     * \brief Returns true iff the error of the solution is below the
     *        tolerance.
     */
320
321
    bool newtonConverged() const
    {
Andreas Lauser's avatar
Andreas Lauser committed
322
323
324
        return 
            error_ <= tolerance_ && 
            model_().jacobianAssembler().reassembleTolerance() <= tolerance_/2;
325
326
327
    }

    /*!
Andreas Lauser's avatar
Andreas Lauser committed
328
329
330
331
     * \brief Called before the newton method is applied to an
     *        non-linear system of equations.
     */
    void newtonBegin(NewtonMethod &method, SolutionVector &u)
332
    {
Andreas Lauser's avatar
Andreas Lauser committed
333
        method_ = &method;
334
335
        numSteps_ = 0;

Andreas Lauser's avatar
Andreas Lauser committed
336
        model_().jacobianAssembler().reassembleAll();
337

Andreas Lauser's avatar
Andreas Lauser committed
338
339
340
341
342
343
344
345
346
347
348
        dtInitial_ = timeManager_().timeStepSize();
        if (enableTimeStepRampUp) {
            rampUpDelta_ = 
                timeManager_().timeStepSize() 
                /
                rampUpSteps()
                *
                2;

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

351
352
353
354
        convergenceWriter_.beginTimestep();
    }

    /*!
Andreas Lauser's avatar
Andreas Lauser committed
355
356
     * \brief Indidicates the beginning of a newton iteration.
     */
357
358
359
360
    void newtonBeginStep()
    { lastError_ = error_; }

    /*!
Andreas Lauser's avatar
Andreas Lauser committed
361
362
363
     * \brief Returns the number of steps done since newtonBegin() was
     *        called.
     */
364
365
366
367
    int newtonNumSteps()
    { return numSteps_; }

    /*!
Andreas Lauser's avatar
Andreas Lauser committed
368
369
370
     * \brief Update the error of the solution compared to the
     *        previous iteration.
     */
371
372
373
374
375
376
377
378
379
    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;
Andreas Lauser's avatar
Andreas Lauser committed
380
        int aboveTol = 0;
381
        for (int i = 0; i < int(uOld.size()); ++i) {
Andreas Lauser's avatar
Andreas Lauser committed
382
383
384
385
386
387
            PrimaryVariables uNewI = uOld[i];
            uNewI -= deltaU[i];
            Scalar vertErr = 
                model_().relativeErrorVertex(i,
                                             uOld[i],
                                             uNewI);
388
            
Andreas Lauser's avatar
Andreas Lauser committed
389
390
391
392
393
394
            if (vertErr > tolerance_)
                ++aboveTol;
            if (vertErr > error_) {
                idxI = i;
                error_ = vertErr;
            }
395
        }
396

Andreas Lauser's avatar
Andreas Lauser committed
397
        error_ = gridView_().comm().max(error_);
398
399
400
    }

    /*!
Andreas Lauser's avatar
Andreas Lauser committed
401
402
403
404
405
406
     * \brief Solve the linear system of equations \f$ \mathbf{A}x - b
     *        = 0\f$.
     *
     * Throws Dumux::NumericalProblem if the linear solver didn't
     * converge.
     */
Bernd Flemisch's avatar
Bernd Flemisch committed
407
408
    template <class Vector>
    void newtonSolveLinear(const JacobianMatrix &A,
409
                           Vector &u,
410
                           const Vector &b)
411
412
413
414
415
416
417
418
419
    {
        // 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
420
            solveLinear_(A, u, b, residReduction);
421

Andreas Lauser's avatar
Andreas Lauser committed
422
423
424
            // make sure all processes converged
            int converged = 1;
            gridView_().comm().min(converged);
425

Andreas Lauser's avatar
Andreas Lauser committed
426
427
428
429
            if (!converged) {
                DUNE_THROW(NumericalProblem,
                           "A process threw MatrixBlockError");
            }
430
431
432
433
        }
        catch (Dune::MatrixBlockError e) {
            // make sure all processes converged
            int converged = 0;
Andreas Lauser's avatar
Andreas Lauser committed
434
            gridView_().comm().min(converged);
435
436
437
438
439
440
441
442
443
444
445

            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
446
447
448
449
450
451
452
453
454
455
456
457
458
459
     * \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
     */
460
461
462
463
464
    void newtonUpdate(SolutionVector &deltaU, const SolutionVector &uOld)
    {
        writeConvergence_(uOld, deltaU);

        newtonUpdateRelError(uOld, deltaU);
Andreas Lauser's avatar
Andreas Lauser committed
465

466
467
468
        deltaU *= -1;
        deltaU += uOld;

Andreas Lauser's avatar
Andreas Lauser committed
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
        // compute the vertex and element colors for partial
        // reassembly
        if (enablePartialReassemble) {
            Scalar maxDelta = 0;
            for (int i = 0; i < int(uOld.size()); ++i) {
                const PrimaryVariables &uEval = this->model_().jacobianAssembler().evalPoint()[i];
                const PrimaryVariables &uSol = this->model_().curSol()[i];
                Scalar tmp = 
                    model_().relativeErrorVertex(i,
                                                 uEval,
                                                 uSol);
                maxDelta = std::max(tmp, maxDelta);
            }
            
            Scalar reassembleTol = std::max(maxDelta/10, this->tolerance_/5);
            if (error_ < 10*tolerance_)
                reassembleTol = tolerance_/5;
            this->model_().jacobianAssembler().computeColors(reassembleTol);
        }               
    }
489

490
    /*!
Andreas Lauser's avatar
Andreas Lauser committed
491
492
     * \brief Indicates that one newton iteration was finished.
     */
493
494
495
    void newtonEndStep(SolutionVector &u, SolutionVector &uOld)
    {
        ++numSteps_;
Andreas Lauser's avatar
Andreas Lauser committed
496
497
498
499
500
501
502
503
504

        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_;
505
        }
506

Andreas Lauser's avatar
Andreas Lauser committed
507
        if (verbose())
508
            std::cout << "\rNewton iteration " << numSteps_ << " done: "
Andreas Lauser's avatar
Andreas Lauser committed
509
                      << "error=" << realError << endIterMsg().str() << "\n";
510
511
512
513
        endIterMsgStream_.str("");
    }

    /*!
Andreas Lauser's avatar
Andreas Lauser committed
514
515
     * \brief Indicates that we're done solving the non-linear system of equations.
     */
516
517
518
519
520
521
    void newtonEnd()
    {
        convergenceWriter_.endTimestep();
    }

    /*!
Andreas Lauser's avatar
Andreas Lauser committed
522
523
524
525
     * \brief Called if the newton method broke down.
     *
     * This method is called _after_ newtonEnd()
     */
526
527
    void newtonFail()
    {
Andreas Lauser's avatar
Andreas Lauser committed
528
        timeManager_().setTimeStepSize(dtInitial_);
529
530
531
532
        numSteps_ = targetSteps_*2;
    }

    /*!
Andreas Lauser's avatar
Andreas Lauser committed
533
534
535
536
     * \brief Called when the newton method was sucessful.
     *
     * This method is called _after_ newtonEnd()
     */
537
538
539
540
    void newtonSucceed()
    { }

    /*!
Andreas Lauser's avatar
Andreas Lauser committed
541
542
543
544
545
546
     * \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.
     */
547
548
    Scalar suggestTimeStepSize(Scalar oldTimeStep) const
    {
Andreas Lauser's avatar
Andreas Lauser committed
549
550
551
        if (enableTimeStepRampUp)
            return oldTimeStep; 

552
        Scalar n = numSteps_;
Andreas Lauser's avatar
Andreas Lauser committed
553
554
        n -= rampUpSteps();

555
556
557
558
559
        // 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.
560
561
        if (n > targetSteps_) {
            Scalar percent = (n - targetSteps_)/targetSteps_;
562
563
564
565
566
            return oldTimeStep/(1.0 + percent);
        }
        else {
            /*Scalar percent = (Scalar(1))/targetSteps_;
              return oldTimeStep*(1 + percent);
Andreas Lauser's avatar
Andreas Lauser committed
567
            */
568
            Scalar percent = (targetSteps_ - n)/targetSteps_;
569
570
571
572
573
            return oldTimeStep*(1.0 + percent/1.2);
        }
    }

    /*!
Andreas Lauser's avatar
Andreas Lauser committed
574
575
576
     * \brief Returns a reference to the current newton method
     *        which is controlled by this controller.
     */
577
578
579
580
    NewtonMethod &method()
    { return *method_; }

    /*!
Andreas Lauser's avatar
Andreas Lauser committed
581
582
583
     * \brief Returns a reference to the current newton method
     *        which is controlled by this controller.
     */
584
585
586
    const NewtonMethod &method() const
    { return *method_; }

Andreas Lauser's avatar
Andreas Lauser committed
587
588
589
    std::ostringstream &endIterMsg()
    { return endIterMsgStream_; }

590
    /*!
Andreas Lauser's avatar
Andreas Lauser committed
591
592
     * \brief Returns true iff the newton method ought to be chatty.
     */
593
    bool verbose() const
Andreas Lauser's avatar
Andreas Lauser committed
594
    { return gridView_().comm().rank() == 0; }
595

Andreas Lauser's avatar
Andreas Lauser committed
596
protected:
597
    /*!
Andreas Lauser's avatar
Andreas Lauser committed
598
599
600
601
     * \brief Returns a reference to the grid view.
     */
    const GridView &gridView_() const
    { return problem_().gridView(); }
602

Andreas Lauser's avatar
Andreas Lauser committed
603
604
605
606
607
608
609
610
611
612
613
614
    /*!
     * \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(); }

615
616
617
618
619
620
    /*!
     * \brief Returns a reference to the time manager.
     */
    TimeManager &timeManager_()
    { return problem_().timeManager(); }

Andreas Lauser's avatar
Andreas Lauser committed
621
622
623
624
625
626
    /*!
     * \brief Returns a reference to the time manager.
     */
    const TimeManager &timeManager_() const
    { return problem_().timeManager(); }

Andreas Lauser's avatar
Andreas Lauser committed
627
628
629
630
631
632
633
634
635
636
637
    /*!
     * \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(); }
638
639
640
641
642
643
644
645
646
647
648
649
650

    // 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
651
        convergenceWriter_.beginIteration(this->gridView_());
652
653
654
655
656
        convergenceWriter_.writeFields(uOld, deltaU);
        convergenceWriter_.endIteration();
    };


Bernd Flemisch's avatar
Bernd Flemisch committed
657
658
    template <class Vector>
    void solveLinear_(const JacobianMatrix &A,
659
660
                      Vector &x,
                      const Vector &b,
661
662
663
                      Scalar residReduction)
    {
        int verbosity = GET_PROP_VALUE(TypeTag, PTAG(NewtonLinearSolverVerbosity));
Andreas Lauser's avatar
Andreas Lauser committed
664
        if (gridView_().comm().rank() != 0)
665
666
667
            verbosity = 0;

#if HAVE_PARDISO
668
        typedef Dumux::PDELab::ISTLBackend_NoOverlap_Loop_Pardiso<TypeTag> Solver;
Andreas Lauser's avatar
Andreas Lauser committed
669
        Solver solver(problem_(), 500, verbosity);
670
671
#else // !HAVE_PARDISO
#if HAVE_MPI
672
//        typedef Dune::PDELab::ISTLBackend_NOVLP_BCGS_NOPREC<GridFunctionSpace> Solver;
Andreas Lauser's avatar
Andreas Lauser committed
673
//        Solver solver(model_().jacobianAssembler().gridFunctionSpace(), 50000, verbosity);
674
        typedef Dumux::PDELab::ISTLBackend_NoOverlap_BCGS_ILU<TypeTag> Solver;
Andreas Lauser's avatar
Andreas Lauser committed
675
        Solver solver(problem_(), 500, verbosity);
676
#else
Andreas Lauser's avatar
Andreas Lauser committed
677
        typedef Dune::PDELab::ISTLBackend_SEQ_BCGS_SSOR Solver;
678
679
680
681
        Solver solver(500, verbosity);
#endif // HAVE_MPI
#endif // HAVE_PARDISO

Andreas Lauser's avatar
Andreas Lauser committed
682
        //    Solver solver(model_().jacobianAssembler().gridFunctionSpace(), 500, verbosity);
683
684
        Vector bTmp(b);
        solver.apply(A, x, bTmp, residReduction);
685
686
687
688
689
690
691
692
693

        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
694
        gridView_().comm().sum(xNorm2);
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
        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_;

    ConvergenceWriter convergenceWriter_;

    Scalar tolerance_;

    Scalar error_;
    Scalar lastError_;

711
712
    // number of iterations for the time-step ramp-up
    Scalar rampUpSteps_;
Andreas Lauser's avatar
Andreas Lauser committed
713
714
715
716
    // the increase of the time step size during the rampup
    Scalar rampUpDelta_;

    Scalar dtInitial_; // initial time step size
717

718
    // optimal number of iterations we want to achive
719
    int targetSteps_;
720
    // maximum number of iterations we do before giving up
721
    int maxSteps_;
722
    // actual number of steps done so far
723
    int numSteps_;
724
};
Andreas Lauser's avatar
Andreas Lauser committed
725
} // namespace Dumux
726
727

#endif