#pragma once #include #include #include #include #include #include #include #include namespace AMDiS { template class PfcPrecon : public Dune::Preconditioner { using SubMatrix = std::decay_t()[_0][_0] )>; using SubVector = std::decay_t()[_0] )>; public: PfcPrecon(Matrix const& matrix, double* tauPtr, double M0) : matrix(matrix) , tauPtr(tauPtr) , M0(M0) , matL0(matrix[_1][_0]) , matL (matrix[_2][_1]) , matM (matrix[_2][_2]) {} virtual void pre(Vector& x, Vector& b) override { double delta = std::sqrt(M0 * (*tauPtr)); matMpL = matM; matMpL.axpy(1.0/delta, matL0); // => MpL = M + 1/delta * L0 matMpL2 = matM; matMpL2.axpy(std::sqrt(delta), matL); y0.resize(matM.N()); y1.resize(matM.N()); tmp.resize(matM.N()); } virtual void apply(Vector& x, const Vector& b) override { double delta = std::sqrt(M0 * (*tauPtr)); solve(matM, y0, b[_0]); // M*y0 = b0 matL0.mv(y0, y1); // y1 = K*y0 tmp = b[_1]; tmp-= y1; // tmp := b1 - tau*y1 solve(matMpL, y1, tmp); // (M + delta*K) * y1 = tmp x[_0] = y0; x[_0].axpy(1.0/delta, y1); // x0 = y0 + (1/delta)*y1 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 x[_0].axpy(-1.0/delta, x[_1]); // x0 = x0 - (1/delta)*x1 = y0 + (1/delta)*(y1 - x1) matL.mv(x[_1], y1); tmp = b[_2]; tmp-= y1; // tmp := b2 - K*x1 solve(matM, x[_2], tmp); } virtual void post(Vector& x) override { } private: template void solve(Mat const& A, Vec& x, Vec b) { Dune::MatrixAdapter op(A); Dune::SeqJac precon(A, 1, 1.0); Dune::CGSolver solver(op, precon, 1.e-3, 20, 0); Dune::InverseOperatorResult statistics; solver.apply(x, b, statistics); } private: Matrix const& matrix; double* tauPtr; double M0; SubMatrix const& matL0; SubMatrix const& matL; SubMatrix const& matM; SubMatrix matMpL; SubMatrix matMpL2; SubVector y0, y1, tmp; }; } // end namespace AMDiS