Nan for stress values

[Raghunandan Pratoori]

Hello,

I am running a code to study phase transformations in metals. The code runs without any problem in release mode, but gives nan when outputting stress values using the same code snippet as in Step-18.

When I run the same code in debug mode, it gives the following error:

An error occurred in line <2183> of file </home/rnp/Software/dealii-build/deal.II-v9.3.3/include/deal.II/lac/affine_constraints.h> in function
    dealii::types::global_dof_index dealii::AffineConstraints<number>::calculate_line_index(dealii::AffineConstraints<number>::size_type) const [with number = double; dealii::types::global_dof_index = unsigned int; dealii::AffineConstraints<number>::size_type = unsigned int]
The violated condition was:
    local_lines.is_element(line_n)
Additional information:
    The index set given to this constraints object indicates constraints
    for degree of freedom 31328 should not be stored by this object, but a
    constraint is being added.

Stacktrace:
-----------
#0  microPF_polycrystal: dealii::AffineConstraints<double>::calculate_line_index(unsigned int) const
#1  microPF_polycrystal: dealii::AffineConstraints<double>::is_constrained(unsigned int) const
#2  /home/rnp/Software/dealii-build/deal.II-v9.3.3/lib/libdeal_II.g.so.9.3.3: dealii::AffineConstraints<double>::make_sorted_row_list(std::vector<unsigned int, std::allocator<unsigned int> > const&, std::vector<unsigned int, std::allocator<unsigned int> >&) const
#3  /home/rnp/Software/dealii-build/deal.II-v9.3.3/lib/libdeal_II.g.so.9.3.3: void dealii::AffineConstraints<double>::add_entries_local_to_global<dealii::DynamicSparsityPattern>(std::vector<unsigned int, std::allocator<unsigned int> > const&, dealii::DynamicSparsityPattern&, bool, dealii::Table<2, bool> const&, std::integral_constant<bool, false>) const
#4  /home/rnp/Software/dealii-build/deal.II-v9.3.3/lib/libdeal_II.g.so.9.3.3: void dealii::AffineConstraints<double>::add_entries_local_to_global<dealii::DynamicSparsityPattern>(std::vector<unsigned int, std::allocator<unsigned int> > const&, dealii::DynamicSparsityPattern&, bool, dealii::Table<2, bool> const&) const
#5  /home/rnp/Software/dealii-build/deal.II-v9.3.3/lib/libdeal_II.g.so.9.3.3: void dealii::DoFTools::make_sparsity_pattern<3, 3, dealii::DynamicSparsityPattern, double>(dealii::DoFHandler<3, 3> const&, dealii::DynamicSparsityPattern&, dealii::AffineConstraints<double> const&, bool, unsigned int)
#6  microPF_polycrystal: PhaseField::Solid<3>::setup_system()
#7  microPF_polycrystal: PhaseField::Solid<3>::run()
#8  microPF_polycrystal: main
--------------------------------------------------------

It points to a problem in defining constraints in setup_system. But I believe I defined the constraints correctly. I have a similar code that works perfectly and the only thing I have changed is adding grain orientations which do not effect the constraint definition.

Is my understanding of the error message correct? Is there something that I am missing in debugging the code?

Below is the code snippet of setup_system -

template <int dim>
    void Solid<dim>::setup_system()
    {
        TimerOutput::Scope t(computing_timer, "setup");

        dof_handler.distribute_dofs(fe);
        dof_handler_c.distribute_dofs(fe_c);
        history_dof_handler.distribute_dofs(history_fe);

        // DoFRenumbering::Cuthill_McKee(dof_handler);

        const unsigned int n_dofs   = dof_handler.n_dofs();
        const unsigned int n_dofs_c = dof_handler_c.n_dofs();

        // Number of active cells and DoFs
        pcout   << "   Number of active cells: "
                << triangulation.n_active_cells()
                << std::endl
                << "   Number of degrees of freedom: "
                << n_dofs + n_dofs_c
                << " (" << n_dofs << '+' << n_dofs_c << ')'
                << std::endl;

        locally_owned_dofs = dof_handler.locally_owned_dofs();
        DoFTools::extract_locally_relevant_dofs(dof_handler,  locally_relevant_dofs);

        constraints.clear();
        constraints.reinit(locally_relevant_dofs);
        {
            DoFTools::make_periodicity_constraints(dof_handler, 0, 1, 0, constraints);
            DoFTools::make_periodicity_constraints(dof_handler, 2, 3, 1, constraints);
            DoFTools::make_periodicity_constraints(dof_handler, 4, 5, 2, constraints);
        }
        constraints.close();
        // constraints.print(std::cout);

        DynamicSparsityPattern dsp(locally_relevant_dofs);
        DoFTools::make_sparsity_pattern(dof_handler, dsp, constraints, false);
        Utilities::MPI::all_gather(mpi_communicator, dof_handler.locally_owned_dofs());
        SparsityTools::distribute_sparsity_pattern (dsp, dof_handler.locally_owned_dofs(),
                                                    mpi_communicator, locally_relevant_dofs);
       
        tangent_matrix.reinit(locally_owned_dofs, locally_owned_dofs, dsp, mpi_communicator);
        solution.reinit(locally_owned_dofs, locally_relevant_dofs, mpi_communicator);
        solution_update.reinit(locally_owned_dofs, locally_relevant_dofs, mpi_communicator);
        system_rhs.reinit(locally_owned_dofs, mpi_communicator);

        locally_owned_dofs_c = dof_handler_c.locally_owned_dofs();
        DoFTools::extract_locally_relevant_dofs(dof_handler_c,  locally_relevant_dofs_c);

        constraints_c.clear();
        constraints_c.reinit(locally_relevant_dofs);
        {
            // DoFTools::make_periodicity_constraints(dof_handler_c, 0, 1, 0, constraints_c);
            // DoFTools::make_periodicity_constraints(dof_handler_c, 2, 3, 1, constraints_c);
            // DoFTools::make_periodicity_constraints(dof_handler_c, 4, 5, 2, constraints_c);
        }
        constraints_c.close();

        DynamicSparsityPattern dsp_c(locally_relevant_dofs_c);
        DoFTools::make_sparsity_pattern(dof_handler_c, dsp_c, constraints_c, true);
        Utilities::MPI::all_gather(mpi_communicator, dof_handler_c.locally_owned_dofs());
        SparsityTools::distribute_sparsity_pattern (dsp_c, dof_handler_c.locally_owned_dofs(),
                                                    mpi_communicator, locally_relevant_dofs_c);
       
        mass_matrix.reinit(locally_owned_dofs_c, locally_owned_dofs_c, dsp_c, mpi_communicator);
        solution_c0.reinit(locally_owned_dofs_c, locally_relevant_dofs_c, mpi_communicator);
        solution_c1.reinit(locally_owned_dofs_c, locally_relevant_dofs_c, mpi_communicator);
        solution_c2.reinit(locally_owned_dofs_c, locally_relevant_dofs_c, mpi_communicator);
        solution_c3.reinit(locally_owned_dofs_c, locally_relevant_dofs_c, mpi_communicator);
        system_rhs_c1.reinit(locally_owned_dofs_c, mpi_communicator);
        system_rhs_c2.reinit(locally_owned_dofs_c, mpi_communicator);
        system_rhs_c3.reinit(locally_owned_dofs_c, mpi_communicator);  }

Thanks in advance,

Raghunandan.

[Wolfgang Bangerth]

Pratoori:

> I am running a code to study phase transformations in metals. The code runs
> without any problem in release mode, but gives nan when outputting stress
> values using the same code snippet as in Step-18.
>
> When I run the same code in debug mode, it gives the following error:

If your program generates an error in debug mode, nothing you compute in
release mode can be trusted. In other words, it is not worth thinking about
the NaN issue until you have figured this out.

> It points to a problem in defining constraints in setup_system.

Actually, it states that the error happens in building the sparsity pattern.
which appears to reference constrained degrees of freedom that aren't stored
in the AffineConstraints object. Why this is so, I can't say -- we'd need to
see the full code, or even better a reduced test case that has everything not
necessary to recreate the code removed.


As for the code snippet, I don't see anything specific that looks wrong except...

...this line: Utilities::MPI::all_gather() returns the combined results of the
operation as a vector. Calling the function without capturing the return does
nothing useful and only serves to burn CPU time.

>         locally_owned_dofs_c = dof_handler_c.locally_owned_dofs();
>         DoFTools::extract_locally_relevant_dofs(dof_handler_c,
>  locally_relevant_dofs_c);
>
>         constraints_c.clear();
>         constraints_c.reinit(locally_relevant_dofs);

And this looks wrong too: You want to initialize constraints_c with
locally_relevant_dofs_c, not locally_relevant_dofs.

Best
W.

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