nonlineargstest.cc 10.5 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
// -*- tab-width: 8; indent-tabs-mode: nil; c-basic-offset: 4 -*-
// vi: set et ts=8 sw=4 sts=4:

#include <config.h>

#include <iostream>
#include <sstream>

// dune-common includes
#include <dune/common/bitsetvector.hh>
#include <dune/common/stringutility.hh>

// dune-istl includes
#include <dune/istl/bcrsmatrix.hh>
15
16


17
18
19
20
21
22
23
24

// dune-grid includes
#include <dune/grid/yaspgrid.hh>
#include <dune/grid/io/file/vtk/vtkwriter.hh>

// dune-solver includes
#include <dune/solvers/solvers/loopsolver.hh>
#include <dune/solvers/norms/energynorm.hh>
25
#include <dune/solvers/common/defaultbitvector.hh>
26
27
#include <dune/solvers/transferoperators/compressedmultigridtransfer.hh>
#include <dune/solvers/iterationsteps/truncatedblockgsstep.hh>
28
29
30
31

// dune-tnnmg includes
#include <dune/tnnmg/functionals/quadraticfunctional.hh>
#include <dune/tnnmg/functionals/boxconstrainedquadraticfunctional.hh>
Carsten Gräser's avatar
Carsten Gräser committed
32
#include <dune/tnnmg/functionals/bcqfconstrainedlinearization.hh>
33
34
#include <dune/tnnmg/iterationsteps/nonlineargsstep.hh>
#include <dune/tnnmg/iterationsteps/acceleratednonlineargsstep.hh>
Carsten Gräser's avatar
Carsten Gräser committed
35
#include <dune/tnnmg/iterationsteps/tnnmgacceleration.hh>
36

37
38
#include <dune/tnnmg/localsolvers/scalarobstaclesolver.hh>

39
40
41
#include <dune/tnnmg/linearsolvers/ldlsolver.hh>
#include <dune/tnnmg/linearsolvers/ilu0cgsolver.hh>
#include <dune/tnnmg/linearsolvers/fastamgsolver.hh>
42
#include <dune/tnnmg/linearsolvers/mgsolver.hh>
43

44
45
#include <dune/tnnmg/projections/obstacledefectprojection.hh>

46
47
48
#include <dune/solvers/test/common.hh>


Carsten Gräser's avatar
Carsten Gräser committed
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
// template<class F>
// std::tuple<double, double, int> bisection(const F& f, double start, double guess, double tol)
// {
//     double a = start;
//     double x = guess;
//     double fx = f(x);
//     int evaluations = 1;
//     while (fx<0)
//     {
//         x = x + (guess-start);
//         fx = f(x);
//         ++evaluations;
//     }
//     double b = x;
//     while (std::fabs(b-a)>tol)
//     {
//         x = (a+b)*.5;
//         fx = f(x);
//         ++evaluations;
//         if (fx<0)
//             a = x;
//         else
//             b = x;
//     }
//     return std::make_tuple(x, fx, evaluations);
// };



78
79
template<class MatrixType, class VectorType, class BitVector, class TransferOperators>
void solveProblem(const MatrixType& mat, VectorType& x, const VectorType& rhs, const BitVector& ignore, const TransferOperators& transfer, int maxIterations=100, double tolerance=1.0e-8)
80
81
82
83
84
85
86
87
88
89
{
    using Vector = VectorType;
    using Matrix = MatrixType;

    auto lower = Vector(rhs.size());
    auto upper = Vector(rhs.size());

    lower = 0;
    upper = 1;

90
    using Functional = Dune::TNNMG::BoxConstrainedQuadraticFunctional<Matrix&, Vector&, Vector&, Vector&, double>;
91
    auto J = Functional(mat, rhs, lower, upper);
Carsten Gräser's avatar
Carsten Gräser committed
92

93
//    auto localSolver = [](auto& xi, const auto& restriction, const auto& ignoreI) {
94
//        gaussSeidelLoop(xi, restriction, ignoreI, Dune::TNNMG::ScalarObstacleSolver());
95
96
97
98
//    };

//    Dune::TNNMG::GaussSeidelSolver<Dune::TNNMG::ScalarObstacleSolver> gss(Dune::TNNMG::ScalarObstacleSolver());

99
    auto localSolver = gaussSeidelLocalSolver(std::make_tuple(Dune::TNNMG::ScalarObstacleSolver(), Dune::TNNMG::ScalarObstacleSolver()));
100
101

//    auto localSolver = gaussSeidelLocalSolver(Dune::TNNMG::ScalarObstacleSolver());
102

103
104

    auto smoother = TruncatedBlockGSStep<Matrix, Vector, BitVector>();
Carsten Gräser's avatar
Carsten Gräser committed
105
    auto linearSolver = Dune::TNNMG::MGSolver<Matrix, Vector>(transfer, smoother, smoother, 3, 3);
106
107

//    auto linearSolver = Dune::TNNMG::FastAMGSolver<Matrix, Vector>();
108
109
//    auto linearSolver = Dune::TNNMG::ILU0CGSolver(1e-5, 5);
//    auto linearSolver = Dune::TNNMG::LDLSolver();
Carsten Gräser's avatar
Carsten Gräser committed
110
111
112

    using Linearization = Dune::TNNMG::BoxConstrainedQuadraticFunctionalConstrainedLinearization<Functional, BitVector>;
    using LinearSolver = decltype(linearSolver);
113
114
    using DefectProjection = Dune::TNNMG::ObstacleDefectProjection;
    using LineSearchSolver = Dune::TNNMG::ScalarObstacleSolver;
Carsten Gräser's avatar
Carsten Gräser committed
115
116
    using Acceleration = Dune::TNNMG::TNNMGAcceleration<Functional, BitVector, Linearization, LinearSolver, DefectProjection, LineSearchSolver>;

117
    auto acceleration = Acceleration(linearSolver, DefectProjection(), LineSearchSolver());
Carsten Gräser's avatar
Carsten Gräser committed
118

119
120
121
122
    using Step = typename Dune::TNNMG::AcceleratedNonlinearGSStep<Vector, Functional, decltype(localSolver), decltype(acceleration)>;
    using Solver = LoopSolver<Vector>;
    using Norm =  EnergyNorm<Matrix, Vector>;

123
124
125
    auto step = std::make_shared<Step>(J, x, localSolver, acceleration);
    auto norm = std::make_shared<Norm>(mat);
    auto solver = Solver(step, 1e9, 0, norm, Solver::FULL);
126

127
128
129
    step->setIgnore(ignore);
    step->setPreSmoothingSteps(3);
    step->setPostSmoothingSteps(3);
130

Carsten Gräser's avatar
Carsten Gräser committed
131
132
133
134
135
136
    solver.addCriterion(
            [&](){
                return Dune::formatString("   % 12.5e", J(x));
            },
            "   energy      ");

137
138
139
140
141
142
143
144
145
146
147
    double initialEnergy = J(x);
    solver.addCriterion(
            [&](){
                static double oldEnergy=initialEnergy;
                double currentEnergy = J(x);
                double decrease = currentEnergy - oldEnergy;
                oldEnergy = currentEnergy;
                return Dune::formatString("   % 12.5e", decrease);
            },
            "   decrease    ");

Carsten Gräser's avatar
Carsten Gräser committed
148
149
    solver.addCriterion(
            [&](){
150
                return Dune::formatString("   % 12.5e", step->accelerationStep().lastDampingFactor());
Carsten Gräser's avatar
Carsten Gräser committed
151
152
153
154
155
156
            },
            "   damping     ");


    solver.addCriterion(
            [&](){
157
                return Dune::formatString("   % 12d", step->accelerationStep().linearization().truncated().count());
Carsten Gräser's avatar
Carsten Gräser committed
158
159
160
161
            },
            "   truncated   ");


162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
//    Vector temp(x.size());
//    auto gradientCriterion = Dune::Solvers::Criterion(
//            [&](){
//                temp = 0;
//                J.addGradient(x, temp);
//                auto err = temp.two_norm();
//                return std::make_tuple(err<tolerance, Dune::formatString("   % 12.5e", err));
//            },
//            "   |gradient|  ");

//    solver.addCriterion(
//            [&](){return Dune::formatString("   % 12.5e", J(x));},
//            "   energy      ");

//    double initialEnergy = J(x);
//    solver.addCriterion(
//            [&](){
//                static double oldEnergy=initialEnergy;
//                double currentEnergy = J(x);
//                double decrease = currentEnergy - oldEnergy;
//                oldEnergy = currentEnergy;
//                return Dune::formatString("   % 12.5e", decrease);
//            },
//            "   decrease    ");

//    solver.addCriterion(
//            [&](){return Dune::formatString("   % 12d", stepRule.evaluations_);},
//            "   evaluations ");

//    solver.addCriterion(
//            [&](){return Dune::formatString("   % 12.5e", stepRule.x_);},
//            "   step size   ");

//    solver.addCriterion(Dune::Solvers::maxIterCriterion(solver, maxIterations) | gradientCriterion);

    std::vector<double> correctionNorms;
198
    solver.addCriterion(Dune::Solvers::correctionNormCriterion(*step, *norm, 1e-10, correctionNorms));
199
200
201

//    solver.addCriterion(Dune::Solvers::correctionNormCriterion(step, norm, 1e-4));

202
    solver.preprocess();
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
    solver.solve();
    std::cout << correctionNorms.size() << std::endl;
}





template <class GridType>
bool checkWithGrid(const GridType& grid, const std::string fileName="")
{
    bool passed = true;

    static const int dim = GridType::dimension;

218
219
    typedef typename Dune::FieldMatrix<double, 2, 2> MatrixBlock;
    typedef typename Dune::FieldVector<double, 2> VectorBlock;
220
221
    typedef typename Dune::BCRSMatrix<MatrixBlock> Matrix;
    typedef typename Dune::BlockVector<VectorBlock> Vector;
222
    typedef typename Dune::Solvers::DefaultBitVector_t<Vector> BitVector;
223
224
225

    typedef typename GridType::LeafGridView GridView;
    typedef typename Dune::FieldVector<double, dim> DomainType;
226
    typedef typename Dune::FieldVector<double, 2> RangeType;
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243


    const GridView gridView = grid.leafGridView();

    Matrix A;
    constructPQ1Pattern(gridView, A);
    A = 0.0;
    assemblePQ1Stiffness(gridView, A);

    Matrix M;
    constructPQ1Pattern(gridView, M);
    M = 0.0;
    assemblePQ1Mass(gridView, M);
    A += M;

    Vector rhs(A.N());
    rhs = 0;
244
    auto f = [](const DomainType& x) -> RangeType { return (x.two_norm() < .7) ? RangeType(200.0) : RangeType(0.0);};
245
246
247
248
249
250
251
    assemblePQ1RHS(gridView, rhs, f);

    Vector u(A.N());
    u = 0;

    BitVector ignore(A.N());
    ignore.unsetAll();
252
    markBoundaryDOFs(gridView, ignore);
253
254


255
256
257
258
259
260
261
262
263
264
265
266
    using TransferOperator = CompressedMultigridTransfer<Vector>;
    using TransferOperators = std::vector<std::shared_ptr<TransferOperator>>;

    TransferOperators transfer(grid.maxLevel());
    for (size_t i = 0; i < transfer.size(); ++i)
    {
        // create transfer operator from level i to i+1
        transfer[i] = std::make_shared<TransferOperator>();
        transfer[i]->setup(grid, i, i+1);
    }

    solveProblem(A, u, rhs, ignore, transfer);
267

268
269
270
271
272
273
//    if (fileName!="")
//    {
//        typename Dune::VTKWriter<GridView> vtkWriter(gridView);
//        vtkWriter.addVertexData(u, "solution");
//        vtkWriter.write(fileName);
//    }
274
275
276
277
278
279
280
281
282
283
284
285
286

    return passed;
}


template <int dim>
bool checkWithYaspGrid(int refine, const std::string fileName="")
{
    bool passed = true;

    typedef Dune::YaspGrid<dim> GridType;

    typename Dune::FieldVector<typename GridType::ctype,dim> L(1.0);
lh1887's avatar
lh1887 committed
287
    typename std::array<int,dim> s;
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
    std::fill(s.begin(), s.end(), 1);

    GridType grid(L, s);

    for (int i = 0; i < refine; ++i)
        grid.globalRefine(1);

    std::cout << "Testing with YaspGrid<" << dim << ">" << std::endl;
    std::cout << "Number of levels: " << (grid.maxLevel()+1) << std::endl;
    std::cout << "Number of leaf nodes: " << grid.leafGridView().size(dim) << std::endl;

    passed = passed and checkWithGrid(grid, fileName);


    return passed;
}




int main(int argc, char** argv) try
{
    Dune::MPIHelper::instance(argc, argv);
    bool passed(true);


//    int refine1d = 16;
    int refine1d = 1;
316
    int refine2d = 6;
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
    int refine3d = 1;

    if (argc>1)
        std::istringstream(argv[1]) >> refine1d;
    if (argc>2)
        std::istringstream(argv[2]) >> refine2d;
    if (argc>3)
        std::istringstream(argv[3]) >> refine3d;

    passed = passed and checkWithYaspGrid<1>(refine1d, "obstacletnnmgtest_yasp_1d_solution");
    passed = passed and checkWithYaspGrid<2>(refine2d, "obstacletnnmgtest_yasp_2d_solution");
//    passed = passed and checkWithYaspGrid<3>(refine3d, "obstacletnnmgtest_yasp_3d_solution");

    return passed ? 0 : 1;
}
332
catch (const Dune::Exception& e) {
333
334

    std::cout << e << std::endl;
335
    return 1;
336
337
338

}