PfcPrecon.hpp 2.43 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
#pragma once

#include <tuple>
#include <cmath>

#include <dune/istl/bcrsmatrix.hh>
#include <dune/istl/bvector.hh>
#include <dune/istl/multitypeblockmatrix.hh>
#include <dune/istl/multitypeblockvector.hh>
#include <dune/istl/preconditioner.hh>

12
13
#include <dune/amdis/common/Mpl.hpp>

14
15
16
17
namespace AMDiS {

template <class Matrix, class Vector>
class PfcPrecon : public Dune::Preconditioner<Vector, Vector>
18
{
19
20
    using SubMatrix = std::decay_t<decltype( std::declval<Matrix>()[_0][_0] )>;
    using SubVector = std::decay_t<decltype( std::declval<Vector>()[_0] )>;
21
22

public:
23
    PfcPrecon(Matrix const& matrix, double* tauPtr, double M0)
24
25
26
27
28
29
	  : matrix(matrix)
	  , tauPtr(tauPtr)
	  , M0(M0)
	  , matL0(matrix[_1][_0])
	  , matL (matrix[_2][_1])
	  , matM (matrix[_2][_2])
30
    {}
31

32
33
    virtual void pre(Vector& x, Vector& b) override
    {
34
		double delta = std::sqrt(M0 * (*tauPtr));
35

36
37
		matMpL = matM;
		matMpL.axpy(1.0/delta, matL0); // => MpL = M + 1/delta * L0
38

39
40
		matMpL2 = matM;
		matMpL2.axpy(std::sqrt(delta), matL);
41

42
43
44
		y0.resize(matM.N());
		y1.resize(matM.N());
		tmp.resize(matM.N());
45
46
47
48
    }

    virtual void apply(Vector& x, const Vector& b) override
    {
49
		double delta = std::sqrt(M0 * (*tauPtr));
50

51
52
53
54
		solve(matM, y0, b[_0]);		// M*y0 = b0
		matL0.mv(y0, y1);		// y1 = K*y0
		tmp = b[_1];
		tmp-= y1;			// tmp := b1 - tau*y1
55

56
57
58
		solve(matMpL, y1, tmp);		// (M + delta*K) * y1 = tmp
		x[_0] = y0;
		x[_0].axpy(1.0/delta, y1);		// x0 = y0 + (1/delta)*y1
59

60
61
62
63
		matM.mv(y1, tmp);		// tmp := M*y1
		solve(matMpL2, y1, tmp);	// (M+eps*sqrt(delta)K) * y1 = tmp
		matM.mv(y1, tmp);		// tmp := M*y1
		solve(matMpL2, x[_1], tmp);	// (M+eps*sqrt(delta)K) * x1 = tmp
64

65
		x[_0].axpy(-1.0/delta, x[_1]);	// x0 = x0 - (1/delta)*x1 = y0 + (1/delta)*(y1 - x1)
66

67
68
69
70
		matL.mv(x[_1], y1);
		tmp = b[_2];
		tmp-= y1;			// tmp := b2 - K*x1
		solve(matM, x[_2], tmp);
71
72
73
74
    }

    virtual void post(Vector& x) override
    {
75

76
    }
77

78

79
private:
80

81
82
83
    template <class Mat, class Vec>
    void solve(Mat const& A, Vec& x, Vec b)
    {
84
85
86
		Dune::MatrixAdapter<Mat, Vec, Vec> op(A);
		Dune::SeqJac<Mat, Vec, Vec> precon(A, 1, 1.0);
		Dune::CGSolver<Vec> solver(op, precon, 1.e-3, 20, 0);
87

88
89
		Dune::InverseOperatorResult statistics;
		solver.apply(x, b, statistics);
90
    }
91

92

93
private:
94

95
    Matrix const& matrix;
96
97
    double* tauPtr;
    double M0;
98

99
100
101
    SubMatrix const& matL0;
    SubMatrix const& matL;
    SubMatrix const& matM;
102

103
104
    SubMatrix matMpL;
    SubMatrix matMpL2;
105

106
107
108
109
    SubVector y0, y1, tmp;
};

} // end namespace AMDiS