AztecOO::Iterate error code -3: loss of precision for offset values in step-15

[Maxi Miller]

I rewrote example 15 for a MPI-environment using Trilinos, and solve it with

IndexSet solution_relevant_partitioning(dof_handler.n_dofs());
        DoFTools::extract_locally_relevant_dofs(dof_handler, solution_relevant_partitioning);
        LinearAlgebraTrilinos::MPI::Vector completely_distributed_solution(dof_handler.locally_owned_dofs(), mpi_communicator);
        LinearAlgebraTrilinos::MPI::Vector completely_distributed_update(dof_handler.locally_owned_dofs(), mpi_communicator);

        completely_distributed_solution = present_solution;

        SolverControl solver_control (dof_handler.n_dofs(),
                                      (system_rhs.l2_norm() > 0) ? 1e-6 * system_rhs.l2_norm() : 1e-6);
        LinearAlgebraTrilinos::SolverGMRES  solver (solver_control);

        LinearAlgebraTrilinos::MPI::PreconditionAMG preconditioner;
        LinearAlgebraTrilinos::MPI::PreconditionAMG::AdditionalData data;

        print_status_update("Initializing system matrix with preconditioner\n");
        preconditioner.initialize(system_matrix, data);
        print_status_update("Solving\n");
        solver.solve (system_matrix, completely_distributed_update, system_rhs,
                      preconditioner);
        print_status_update("Solving done\n");

        hanging_node_constraints.distribute (completely_distributed_update);


        print_status_update("Adding to newton\n");
        newton_update = completely_distributed_update;

        newton_update.compress(VectorOperation::insert);

        const double alpha = determine_step_length();
        completely_distributed_solution.add (alpha, completely_distributed_update);

        present_solution = completely_distributed_solution;

which works fine. Now, as next step I wanted to test if I could introduce an offset, i.e. the resulting surface is not located at z=0, but rather at z=OFFSET, with OFFSET a double value set at the beginning. Thus I also initialize the solution accordingly with

    template <int dim>
    class InitialValues : public Function<dim>
    {
        public:
            InitialValues () : Function<dim>() {}
            virtual double value(const Point<dim> &p, const unsigned int component) const;
    };

    template <int dim>
    double InitialValues<dim>::value(const Point<dim> &p, const unsigned int component) const
    {
        return OFFSET;
    }                

//In run-function
                LinearAlgebraTrilinos::MPI::Vector local_solution;
                local_solution.reinit(dof_handler.locally_owned_dofs(), mpi_communicator);
                VectorTools::project (dof_handler,
                                      boundary_constraints,
                                      QGauss<dim>(fe.degree+2),
                                      InitialValues<dim>(),
                                      local_solution);
                //boundary_constraints.distribute(local_solution);
                present_solution = local_solution;

Now I notice that if my offset goes above a certain level, I have to reduce the step length, else I get the error

An error occurred in line <522> of file </home/roland/Downloads/dealii/source/lac/trilinos_solver.cc> in function
    void dealii::TrilinosWrappers::SolverBase::do_solve(const Preconditioner&) [with Preconditioner = dealii::TrilinosWrappers::PreconditionBase]
The violated condition was:
    false
Additional information:
    AztecOO::Iterate error code -3: loss of precision

I do not understand the reason why that happens. Is it a mathematical problem? Or rather a problem in my code? Afaik the gradients should not depend on the offset, thus I do not know where to look for the problem here.
Thanks!

[Wolfgang Bangerth]

On 09/25/2017 01:25 PM, 'Maxi Miller' via deal.II User Group wrote:
>
> I do not understand the reason why that happens. Is it a mathematical
> problem? Or rather a problem in my code? Afaik the gradients should not
> depend on the offset, thus I do not know where to look for the problem here.

Did you adjust the boundary values in the same way as the initial
values? Because if you don't, your initial guess will have a sharp
boundary layer where it drops from OFFSET to zero, and that leads to a
very large non-linearity.

Best
W.

--
------------------------------------------------------------------------
Wolfgang Bangerth email: bang...@colostate.edu
www: http://www.math.colostate.edu/~bangerth/

[Maxi Miller]

My BoundaryValues-function looks like

template <int dim>
    class BoundaryValues : public Function<dim>
    {
        public:
            BoundaryValues () : Function<dim>() {}

            virtual double value (const Point<dim>   &p,
                                  const unsigned int  component = 0) const;
    };


    template <int dim>
    double BoundaryValues<dim>::value (const Point<dim> &p,
                                       const unsigned int /*component*/) const
    {
        return std::sin(2 * numbers::PI * (p[0]+p[1])) + OFFSET;
    }

Thus I do not think that that is the problem. Assembly of the matrix and the residual is done using a second constraint matrix called "newton_constraints", defined as

newton_constraints.clear();
        newton_constraints.reinit(solution_relevant_partitioning);
        DoFTools::make_hanging_node_constraints(dof_handler, newton_constraints);
        VectorTools::interpolate_boundary_values(dof_handler, 0, ZeroFunction<dim>(), newton_constraints);
        newton_constraints.close();

in the following way:

newton_constraints.distribute_local_to_global(cell_residual, local_dof_indices, local_residual);

or

newton_constraints.distribute_local_to_global(cell_matrix, cell_rhs, local_dof_indices, system_matrix, system_rhs);

in order to keep the zero boundary conditions for the newton update.

[Maxi Miller]

Your assumption is correct, I checked the shape of the surface after initialization. It looks (for OFFSET=2 and refine_global(2)) like

Setting the OFFSET-value to 10 while leaving refine_global at 2 results in

Here the drop is clearly visible (which presumably results in the calculation problems), but can not be removed with higher values for refine_global. As example using refine_global(4) and OFFSET=10:

I initialize my solution using the code shown above, but still the sharp drop remains. Why?
Thank you!

Am Dienstag, 26. September 2017 03:15:47 UTC+2 schrieb Wolfgang Bangerth:

[Wolfgang Bangerth]

On 09/26/2017 02:24 AM, 'Maxi Miller' via deal.II User Group wrote:
>
>
> I initialize my solution using the code shown above, but still the sharp drop
> remains. Why?

I don't know. But at least now you know what is happening, and that gives you
an angle to debugging the problem. Are you outputting the initial guess for
the solution, or the solution after the first Newton step?

Btw, the reason your solver fails in the presence of such sharp drops is that
the coefficient in the equation has the form
1 / sqrt( 1 + |grad u|^2 )
In the pictures you show, you have areas where grad u = 0 and other areas
where grad u is very large. So the coefficient varies greatly. The condition
number of your matrix is proportional to the ratio of the largest to the
smallest value of the coefficient, so a sharp drop in your solution means that
you end up with a problem with a large condition number. That's why your
solver fails.