multidomainnewtoncontroller.hh 11.1 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
// vi: set et ts=4 sw=4 sts=4:
/*****************************************************************************
 *   See the file COPYING for full copying permissions.                      *
 *                                                                           *
 *   This program is free software: you can redistribute it and/or modify    *
 *   it under the terms of the GNU General Public License as published by    *
 *   the Free Software Foundation, either version 2 of the License, or       *
 *   (at your option) any later version.                                     *
 *                                                                           *
 *   This program is distributed in the hope that it will be useful,         *
 *   but WITHOUT ANY WARRANTY; without even the implied warranty of          *
 *   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the            *
 *   GNU General Public License for more details.                            *
 *                                                                           *
 *   You should have received a copy of the GNU General Public License       *
 *   along with this program.  If not, see <http://www.gnu.org/licenses/>.   *
 *****************************************************************************/
/*!
* \file
21
* \brief Additional properties required for the coupled Newton controller.
22
*/
23
24
#ifndef DUMUX_MULTIDOMAIN_NEWTON_CONTROLLER_HH
#define DUMUX_MULTIDOMAIN_NEWTON_CONTROLLER_HH
25

26
#include <dumux/nonlinear/newtoncontroller.hh>
27
#include "multidomainconvergencewriter.hh"
28
29
30
31

namespace Dumux
{
template <class TypeTag>
32
class MultiDomainNewtonController;
33

34
template <class TypeTag>
35
struct MultiDomainConvergenceWriter;
36

37
38
namespace Properties
{
39
// set default values for Newton for multidomain problems
40
// they can be overwritten in the parameter file
41
42
43
SET_INT_PROP(MultiDomain, NewtonTargetSteps, 8);
SET_INT_PROP(MultiDomain, NewtonMaxSteps, 15);
SET_SCALAR_PROP(MultiDomain, NewtonRelTolerance, 1e-5);
44
45
46
47
48
49
50
51
52
53
54
55
}


/*!
 * \brief Reference implementation of a newton controller for coupled problems.
 *
 * 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.
 */
template <class TypeTag>
56
class MultiDomainNewtonController : public NewtonController<TypeTag>
57
{
58
    typedef NewtonController<TypeTag> ParentType;
59
60
61
62
63
64
65
66
    typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
    typedef typename GET_PROP_TYPE(TypeTag, NewtonController) Implementation;

    typedef typename GET_PROP_TYPE(TypeTag, NewtonMethod) NewtonMethod;
    typedef typename GET_PROP_TYPE(TypeTag, SolutionVector) SolutionVector;
    typedef typename GET_PROP_TYPE(TypeTag, Problem) Problem;
    typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;

67
68
    typedef typename GET_PROP_TYPE(TypeTag, SubDomain1TypeTag) SubDomain1TypeTag;
    typedef typename GET_PROP_TYPE(TypeTag, SubDomain2TypeTag) SubDomain2TypeTag;
69

70
71
    typedef typename GET_PROP_TYPE(SubDomain1TypeTag, GridView) GridView1;
    typedef typename GET_PROP_TYPE(SubDomain2TypeTag, GridView) GridView2;
72

73
    typedef MultiDomainConvergenceWriter<TypeTag>  ConvergenceWriter;
74
75
76
    typedef typename GET_PROP_TYPE(TypeTag, LinearSolver) LinearSolver;

public:
77
78
79
80
81
    /*!
     * \brief Constructor
     *
     * \param problem The problem
     */
82
    MultiDomainNewtonController(const Problem &problem)
83
84
        : ParentType(problem)
		, endIterMsgStream_(std::ostringstream::out)
85
86
87
88
        , linearSolver_(problem)
        , convergenceWriter_(asImp_())
    {
//		  Writes out, where the relative tolerance is defined
89
90
91
92
93
        std::cout << "ParameterNewtonRelTol= "
        		<< PROP_DIAGNOSTIC(TypeTag, NewtonRelTolerance)
        		<< ", "
        		<< GET_PROP_VALUE(TypeTag, NewtonRelTolerance)
        		<< std::endl;
94
95
96
97
98
    };

    /*!
     * \brief Destructor
     */
99
    ~MultiDomainNewtonController()
100
    { };
101
102
103
104
105
106


    /*!
    * \brief Update the error of the solution compared to the
    *        previous iteration.
    *
107
108
109
110
    * \param uLastIter The solution of the last iteration
    * \param deltaU The delta as calculated from solving the linear
    *               system of equations. This parameter also stores
    *               the updated solution.
111
112
113
114
115
116
    */
    void newtonUpdateRelError(const SolutionVector &uLastIter,
                              const SolutionVector &deltaU)
    {
        // calculate the relative error as the maximum relative
        // deflection in any degree of freedom.
117
        this->error_ = 0;
118
119
120
121

        SolutionVector uNewI = uLastIter;
        uNewI -= deltaU;

122
123
124
125
        for (unsigned int i = 0; i < uLastIter.base().size(); ++i) {
            for (unsigned int j = 0; j < uLastIter.base()[i].size(); ++j) {
                Scalar vertexError = std::abs(deltaU.base()[i][j]);
                vertexError /= std::max<Scalar>(1.0, std::abs(uLastIter.base()[i][j] + uNewI.base()[i][j])/2);
126

127
                this->error_ = std::max(this->error_, vertexError);
128
129
130
131
132
            }
        }
    }

    /*!
133
134
     * \brief Solve the linear system of equations
     *        \f$ \mathbf{A} x - b = 0\f$.
135
     *
136
137
138
     * \param A Coefficient matrix A
     * \param x Vector of unknowns
     * \param b Right hand side
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
     *
     * Throws Dumux::NumericalProblem if the linear solver didn't
     * converge.
     */
    template <class Matrix, class Vector>
    void newtonSolveLinear(Matrix &A,
                           Vector &x,
                           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...
        try {
154
            int converged = linearSolver_.solve(A.base(), x.base(), b.base());
155
156

#if HAVE_MPI
157
            // make sure all processes converged
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
            int convergedSend = 1;
            MPI_Allreduce(/*sendBuf=*/&convergedSend,
                          /*recvBuf=*/&converged,
                          /*count=*/1,
                          MPI_INT,
                          MPI_MIN,
                          MPI_COMM_WORLD);
#endif
            if (!converged) {
                DUNE_THROW(NumericalProblem,
                           "Linear solver did not converge");
            }
        }
        catch (const Dune::MatrixBlockError &e) {
#if HAVE_MPI
173
            // make sure all processes converged
174
175
176
177
178
179
180
181
182
183
184
185
186
187
            int convergedSend = 0;
            int converged;

            MPI_Allreduce(/*sendBuf=*/&convergedSend,
                          /*recvBuf=*/&converged,
                          /*count=*/1,
                          MPI_INT,
                          MPI_MIN,
                          MPI_COMM_WORLD);
#endif

            Dumux::NumericalProblem p;
            std::string msg;
            std::ostringstream ms(msg);
188
            ms << e.what() << "M=" << A.base()[e.r][e.c];
189
190
191
192
193
            p.message(ms.str());
            throw p;
        }
        catch (const Dune::Exception &e) {
#if HAVE_MPI
194
            // make sure all processes converged
195
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
            int convergedSend = 0;
            int converged;

            MPI_Allreduce(/*sendBuf=*/&convergedSend,
                          /*recvBuf=*/&converged,
                          /*count=*/1,
                          MPI_INT,
                          MPI_MIN,
                          MPI_COMM_WORLD);
#endif

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

    /*!
    * \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 uLastIter.
    *
222
223
    * \param uCurrentIter The solution of the current iteration
    * \param uLastIter The solution of the last iteration
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
    * \param deltaU The delta as calculated from solving the linear
    *               system of equations. This parameter also stores
    *               the updated solution.
    *
    */
    void newtonUpdate(SolutionVector &uCurrentIter,
                      const SolutionVector &uLastIter,
                      const SolutionVector &deltaU)
    {
        if (GET_PARAM_FROM_GROUP(TypeTag, bool, Newton, WriteConvergence)) {
            writeConvergence_(uLastIter, deltaU);
        }

        newtonUpdateRelError(uLastIter, deltaU);

        uCurrentIter = uLastIter;
        uCurrentIter -= deltaU;
    }

    /*!
    * \brief Indicates that one newton iteration was finished.
    *
246
    * \param uCurrentIter The solution of the current iteration
247
248
249
250
251
    * \param uLastIter The solution of the last iteration
    *
    */
    void newtonEndStep(SolutionVector &uCurrentIter, SolutionVector &uLastIter)
    {
252
        typedef Dumux::SplitAndMerge<TypeTag> Common;
253

254
        Common::splitSolVector(this->model_().curSol(),
255
256
                               this->model_().sdModel1().curSol(),
                               this->model_().sdModel2().curSol());
257

258
        ParentType::newtonEndStep(uCurrentIter, uLastIter);
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
    }

    /*!
    * \brief Called when the newton method was sucessful.
    *
    * This method is called _after_ newtonEnd()
    */
    void newtonSucceed()
    {
    }

    /*!
    * \brief the convergence writer produces the output
    *
    * \param uLastIter The solution of the last iteration
    * \param deltaU The delta as calculated from solving the linear
    *               system of equations. This parameter also stores
    *               the updated solution.
    *
    */
    void writeConvergence_(const SolutionVector &uLastIter,
                           const SolutionVector &deltaU)
    {
        if (GET_PARAM_FROM_GROUP(TypeTag, bool, Newton, WriteConvergence)) {
283
            convergenceWriter_.beginIteration(sdGridView1_(), sdGridView2_());
284
285
286
287
288
289
290
291
            convergenceWriter_.writeFields(uLastIter, deltaU);
            convergenceWriter_.endIteration();
        };
    };

    /*!
     * \brief the subdomain gridviews
     */
292
293
294
295
    const GridView1 sdGridView1_() const
    { return this->problem_().sdGridView1(); }
    const GridView2 sdGridView2_() const
    { return this->problem_().sdGridView2(); }
296

297

298
private:
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
    Implementation &asImp_()
    { return *static_cast<Implementation*>(this); }
    const Implementation &asImp_() const
    { return *static_cast<const Implementation*>(this); }

    bool verbose_;

    std::ostringstream endIterMsgStream_;
    NewtonMethod *method_;

    // optimal number of iterations we want to achieve
    int targetSteps_;
    // maximum number of iterations we do before giving up
    int maxSteps_;
    // actual number of steps done so far
    int numSteps_;

    // the linear solver
    LinearSolver linearSolver_;

    ConvergenceWriter convergenceWriter_;
};

} // namespace Dumux

324
#endif // DUMUX_MULTIDOMAIN_NEWTON_CONTROLLER_HH