Dear deal.ii,
I am having a problem with the CG Solver. Here is a code snippet:
SolverControl solver_control
(solution[ts].block(1).size(),1e-8*schur_rhs.l2_norm());
std::cout << "iterations = " << solution[ts].block(1).size() <<
std::endl;
std::cout << "tolerance = " << 1e-8*schur_rhs.l2_norm()<<std::endl;
SolverCG<> cg (solver_control);
cg.solve (schur_complement, solution[ts].block(1), schur_rhs,
preconditioner);
This is the output:
iterations = 48
tolerance = 4.92837e-10
----------------------------------------------------
Exception on processing:
Iterative method reported convergence failure in step 48 with residual
3.18991e-17
Aborting!
----------------------------------------------------
This seems strange as the algorithm is clearly achieving the tolerance.
If I reduce the tolerance to 1e-7*schur_rhs.l2_norm() then the algorithm
reaches the tolerance in 11 iterations and continues with the code.
Further interesting behaviour is that if I adjust the iterations by
hand, e.g., SolverControl solver_control (100,1e-8*schur_rhs.l2_norm())
the code still terminates after 48 iterations (or
solution[ts].block(1).size() if the problem size is altered). This is
also the case if I just use a default solver_control.
I will post the full solver code below which is based on that of step-21.
Many thanks,
John Chapman
template <int dim>
void TwoPhaseFlowProblem<dim>::solve_velocity_pressure (const unsigned
int ts)
{
const InverseMatrix<SparseMatrix<double> >
m_inverse(system_matrix.block(0,0));
// vector size of velocity components
Vector<double> tmp (solution[ts].block(0).size());
// vector size of pressure components
Vector<double> schur_rhs (solution[ts].block(1).size());
// Solve for pressure using Schur Complement
{
m_inverse.vmult (tmp, system_rhs.block(0));
system_matrix.block(1,0).vmult (schur_rhs, tmp);
schur_rhs -= system_rhs.block(1);
SchurComplement schur_complement (system_matrix, m_inverse);
ApproximateSchurComplement approximate_schur_complement
(system_matrix);
InverseMatrix<ApproximateSchurComplement>
preconditioner (approximate_schur_complement);
SolverControl solver_control (solution[ts].block(1).size(),
1e-8*schur_rhs.l2_norm());
std::cout << "iterations = " << solution[ts].block(1).size() <<
std::endl;
std::cout << "tolerance = " << 1e-8*schur_rhs.l2_norm()<<std::endl;
SolverCG<> cg (solver_control);
cg.solve (schur_complement, solution[ts].block(1), schur_rhs,
preconditioner);
std::cout << " "
<< solver_control.last_step()
<< " CG Schur complement iterations for pressure." << std::endl;
}
// Now the velocity using an inverse
{
system_matrix.block(0,1).vmult (tmp, solution[ts].block(1));
tmp *= -1;
tmp += system_rhs.block(0);
m_inverse.vmult (solution[ts].block(0), tmp);
}
constraints.distribute (solution[ts]);
}
_______________________________________________
dealii mailing list http://poisson.dealii.org/mailman/listinfo/dealii
