Dear John,
At a glance I think that the problem could be that you're trying to work
to a tolerance of less than the machine tolerance. The smallest
difference between two values that can be represented on a standard
machine is ~1.1e-16. Although one can store a value smaller than this,
it may be that the operations required to solve with the CG method are
not able to produce a solution with the end precision that you require
and the final residual is tiny but cannot decrease to the margin you
want. This would explain why if you raise the tolerance by an order of
magnitude then convergence can be achieved - the norm of the schur RHS
vector must be very small already.
You can check what the actual value is with the numerical limits call:
#include< limits>
std::numerical_limits<double>::epsilon()
See http://www.cplusplus.com/reference/std/limits/numeric_limits/ ,
http://www.cplusplus.com/reference/clibrary/cfloat/ for more info.
It is unclear to me as to why the solver wouldn't run to full number of
allowed iterations though -- between the CG solver and solver control
classes, it doesn't appear as if there is a check/throw for divergence
or slow convergence that would stop the solver iterating.
Regards,
J-P
On 20/12/2011 13:15, John Chapman wrote:
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
_______________________________________________
dealii mailing list http://poisson.dealii.org/mailman/listinfo/dealii
