Re: AW: [deal.II] step-42 clarification

[Timo Heister]

> in the code and re-implemented it. In serial version, all works fine so far.
> However, when running in parallel, I am seeing an issue in the method
> PlasticityContactProblem::update_solution_and_constraints.
>
> In particular, it turns out that the value of
>
> const unsigned int index_z = dof_indices[q_point];
>
> might be out of the range of

If you do a loop over all locally owned and locally relevant cells
than all dof values of a ghosted vector should exist. If you see an
error, something else must be incorrect (like the IndexSets).

>   PETScWrappers::MPI::Vector lambda( this->locally_relevant_dofs, 
> this->mpi_communicator);

This looks suspicious. Does this really create a ghosted vector in
PETSc? I thought this would fail (at least in debug mode).

Finally, it looks like you modified it to only look at locally owned
cells to build constraints. The problem with this is that processors
also need to know about constraints on ghost cells, not only locally
owned cells. You no longer compute them, which means the solution
might become incorrect around processor boundaries. It probably
(hopefully?) works without adaptivity because each locally owned DoF
is within at least one locally owned cell, but imagine a case where a
dof on a ghost cells is constrained and interacts with a hanging node
the current processor owns. You will not handle this case correctly.

I don't quite remember if there is an easy way to do this, but I
remember writing a debug function that checks if a ConstraintMatrix is
consistent in parallel. This was a while back, but I can try to find
it.

-- 
Timo Heister
http://www.math.clemson.edu/~heister/

-- 
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see 
https://groups.google.com/d/forum/dealii?hl=en
--- 
You received this message because you are subscribed to the Google Groups 
"deal.II User Group" group.
To unsubscribe from this group and stop receiving emails from it, send an email 
to dealii+unsubscr...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

Re: AW: [deal.II] step-42 clarification

[Alberto Salvadori]

Hi Wolfgang and Jörg

I am taking advantage for a second question. I have made a few minor
changes in the code and re-implemented it. In serial version, all works
fine so far. However, when running in parallel, I am seeing an issue in the
method PlasticityContactProblem::update_solution_and_constraints.

In particular, it turns out that the value of

const unsigned int index_z = dof_indices[q_point];

might be out of the range of

TrilinosWrappers::MPI::Vector

lambda(locally_relevant_dofs, mpi_communicator);
lambda = newton_rhs_uncondensed;

I have solved run-time issues by checking if

if ( lambda.in_local_range(index_z))

but I wonder if the issue is possible (in principle, I mean) for the
published code or is it my implementation that is not correct? Note that I
am basically implementing step-42 upon the code built in step-20 using
PETScWrappers::MPI::Vector, rather than step-40. If of use, I am writing my
method below.

Thanks for your help.

Alberto


::update_solution_and_constraints ()


// The following function is the first function we call in each outer loop
of the Newton iteration.

// What it does is to project the solution onto the feasible set and update
the active set for the degrees of freedom that touch or penetrate the
obstacle.


{



  std::pair  vrange;



  this->IO.log() << " In update_solution_and_constraints. " << std::flush;



  // We need to write into the solution vector (which we can only do with
fully distributed vectors without ghost elements) and we need to read the
Lagrange multiplier and the elements of the diagonal mass matrix from their
respective vectors (which we can only do with vectors that do have ghost
elements), so we create the respective vectors. We then also initialize the
constraints object that will contain constraints from contact, as well as
an object that contains an index set of all locally owned degrees of
freedom that are part of the contact:



  std::vector dof_touched( this->dof_handler.n_dofs(), false);



  PETScWrappers::MPI::Vector distributed_solution( this->locally_owned_dofs,
this->mpi_communicator );

  distributed_solution = this->accumulated_displacement;

  vrange = distributed_solution.local_range();

  this->IO.log() << " distributed_solution local range =  " << vrange.first
<< ", " << vrange.second << "\n"  << std::flush;




  PETScWrappers::MPI::Vector lambda( this->locally_relevant_dofs, this
->mpi_communicator);

  lambda = this->system_rhs_uncondensed;

  vrange = lambda.local_range();

  this->IO.log() << " lambda local range =  " << vrange.first << ", " <<
vrange.second << "\n"  << std::flush;




  PETScWrappers::MPI::Vector diag_mass_matrix_vector_relevant(
this->locally_relevant_dofs,
this->mpi_communicator);

  diag_mass_matrix_vector_relevant = this->diag_mass_matrix_vector;

  vrange = diag_mass_matrix_vector_relevant.local_range();

  this->IO.log() << " diag_mass_matrix_vector_relevant local range =  " <<
vrange.first << ", " << vrange.second << "\n"  << std::flush;




  this->all_constraints.reinit( this->locally_relevant_dofs );

  this->contact_constraints.reinit( this->locally_relevant_dofs );

  this->active_set.clear();



  // The second part is a loop over all cells in which we look at each
point where a degree of freedom is defined whether the active set condition
is true and we need to add this degree of freedom to the active set of
contact nodes. As we always do, if we want to evaluate functions at
individual points, we do this with an FEValues object (or, here, an
FEFaceValues object since we need to check contact at the surface) with an
appropriately chosen quadrature object. We create this face quadrature
object by choosing the "support points" of the shape functions defined on
the faces of cells (for more on support points, see this glossary entry).
As a consequence, we have as many quadrature points as there are shape
functions per face and looping over quadrature points is equivalent to
looping over shape functions defined on a face. With this, the code looks
as follows:



  Quadrature face_quadrature( this
->fe.get_unit_face_support_points());

  FEFaceValues fe_values_face( this->fe, face_quadrature,
update_quadrature_points );

  const unsigned int dofs_per_face = this->fe.dofs_per_face;

  const unsigned int n_face_q_points = face_quadrature.size();

  std::vector dof_indices(dofs_per_face);



  typename DoFHandler::active_cell_iterator

  cell = this->dof_handler.begin_active(),

  endc = this->dof_handler.end();

  for (; cell != endc; ++cell)

if (!cell->is_artificial())

{

  if (cell->is_locally_owned())  // MY OWN CODE

  {

for (unsigned int face=0; faceface(face)->at_boundary()

  &&