step-42 clarification

[Alberto Salvadori]

Dear community,

I am studying with pleasure step-42 and I got a bit confused in the method PlasticityContactProblem::solve_newton. At the very end of it, we find these instructions:

old_active_set = active_set;

previous_residual_norm = residual_norm;

 

if (Utilities::MPI::sum((active_set == old_active_set) ? 0 : 1,

mpi_communicator) == 0)

{

pcout << " Active set did not change!" << std::endl;

if (residual_norm < 1e-10)

break;

}

Now, I am a bit confused by the if statement and I am clearly missing something. The first instruction makes old_active_set equal to active_set. Accordingly, wouldn't active_set == old_active_set  always be true? What am I missing?

Thanks for your help

Alberto

[Wolfgang Bangerth]

On 01/16/2018 12:59 AM, Alberto Salvadori wrote:
> Dear community,
>
> I am studying with pleasure step-42 and I got a bit confused in the
> method PlasticityContactProblem::solve_newton. At the very end of it, we
> find these instructions:
>
> old_active_set = active_set;
> previous_residual_norm = residual_norm;
> if (Utilities::MPI::sum

> <http://www.dealii.org/8.5.0/doxygen/deal.II/namespaceUtilities_1_1MPI.html#ab544a3bf3301a6dd3e705ee352c5551b>((active_set

> == old_active_set) ? 0 : 1,
> mpi_communicator) == 0)
> {
> pcout << " Active set did not change!" << std::endl;
> if (residual_norm < 1e-10)
> break;
> }
>
>
> Now, I am a bit confused by the if statement and I am clearly missing
> something. The first instruction makes old_active_set equal to
> active_set. Accordingly, wouldn't active_set == old_active_set  always
> be true? What am I missing?

You're right -- the code would suggest that the if-condition is always
true and that what matter is only the inner if.

I'm not sure I have the current email address of Joerg Frohne who wrote
step-42. He would be the one who could possibly provide a better answer.

Best
W.

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

[Frohne, Joerg]

Dear Wolfgang, dear community,

obviously this is a mistake. I checked the source files which we have used for results in corresponding paper.
There I have found the following coding:

resid_old = resid;

resid_vector = system_rhs_newton;
resid_vector.compress(VectorOperation::insert);

int is_my_set_changed = (active_set == active_set_old) ? 0 : 1;
int num_changed = Utilities::MPI::sum(is_my_set_changed,
MPI_COMM_WORLD);
if (num_changed == 0 && resid < 1e-8)
break;
active_set_old = active_set;

Unfortunately I have not been working with deal.ii for over two years now.
For this reason I would be glad if someone else could fix this bug.
Otherwise I would need a short manual how to do it by myself.

With the best wishes for the future for deal.ii and their developers!
Joerg

________________________________________
Von: Wolfgang Bangerth [bang...@colostate.edu]
Gesendet: Dienstag, 16. Januar 2018 21:48
An: dea...@googlegroups.com; Frohne, Joerg
Cc: Suttmeier, Franz-Theo
Betreff: Re: [deal.II] step-42 clarification

[Wolfgang Bangerth]

On 01/17/2018 02:45 PM, Frohne, Joerg wrote:
> obviously this is a mistake. I checked the source files which we have used for results in corresponding paper.
> There I have found the following coding:
>
> resid_old = resid;
>
> resid_vector = system_rhs_newton;
> resid_vector.compress(VectorOperation::insert);
>
> int is_my_set_changed = (active_set == active_set_old) ? 0 : 1;
> int num_changed = Utilities::MPI::sum(is_my_set_changed,
> MPI_COMM_WORLD);
> if (num_changed == 0 && resid < 1e-8)
> break;
> active_set_old = active_set;
>
> Unfortunately I have not been working with deal.ii for over two years now.
> For this reason I would be glad if someone else could fix this bug.
> Otherwise I would need a short manual how to do it by myself.

I can do that for you. Do I understand correctly that the only change
that needs to happen is to move
active_set_old = active_set;
to *after* the if-statement?

Hope you're doing alright! Best
Wolfgang

[Frohne, Joerg]

>>On 01/17/2018 02:45 PM, Frohne, Joerg wrote:
>> obviously this is a mistake. I checked the source files which we have used for results in corresponding paper.
>> There I have found the following coding:
>>
>> resid_old = resid;
>>
>> resid_vector = system_rhs_newton;
>> resid_vector.compress(VectorOperation::insert);
>>
>> int is_my_set_changed = (active_set == active_set_old) ? 0 : 1;
>> int num_changed = Utilities::MPI::sum(is_my_set_changed,
>> MPI_COMM_WORLD);
>> if (num_changed == 0 && resid < 1e-8)
>> break;
>> active_set_old = active_set;
>>
>> Unfortunately I have not been working with deal.ii for over two years now.
>> For this reason I would be glad if someone else could fix this bug.
>> Otherwise I would need a short manual how to do it by myself.

>I can do that for you. Do I understand correctly that the only change
>that needs to happen is to move
> active_set_old = active_set;
>to *after* the if-statement?

Exactly.

Best,
Joerg

[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 <unsigned, unsigned> 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<bool> 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<dim-1> face_quadrature( this->fe.get_unit_face_support_points());

  FEFaceValues<dim> 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<types::global_dof_index> dof_indices(dofs_per_face);

  

  typename DoFHandler<dim>::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; face<GeometryInfo<dim>::faces_per_cell; ++face)

          

          // the Contact boundary is set to value 3, i.e. set_boundary_id (3);

          

          if (cell->face(face)->at_boundary()

              &&

              cell->face(face)->boundary_id() == 3)

          {

            fe_values_face.reinit(cell, face);

            cell->face(face)->get_dof_indices(dof_indices);

            for (unsigned int q_point=0; q_point<n_face_q_points; ++q_point)

            {

              // At each quadrature point (i.e., at each support point of a degree of freedom located on the contact boundary), we then ask whether it is part of the X2-displacement degrees of freedom and if we haven't encountered this degree of freedom yet (which can happen for those on the edges between faces), we need to evaluate the gap between the deformed object and the obstacle. If the active set condition is true, then we add a constraint to the ConstraintMatrix object that the next Newton update needs to satisfy, set the solution vector's corresponding element to the correct value, and add the index to the IndexSet object that stores which degree of freedom is part of the contact:

              

              const unsigned int component = this->fe.face_system_to_component_index(q_point).first;

              const unsigned int index_X2 = dof_indices[q_point];

              

              // MY OWN CODE: component has been changed to 1 since the direction of contact is X2 for the cell.

              

              if ((component == 1) && (dof_touched[index_X2] == false))

              {

                dof_touched[index_X2] = true;

                const Point<dim> this_support_point = fe_values_face.quadrature_point(q_point);

                const double obstacle_value = this->obstacle->value(this_support_point, 1);

                

                const double solution_here = this->accumulated_displacement(index_X2); // solution(index_X2);

                

                // I believe it should be undeformed_gap >=0, i.e. :

                const double undeformed_gap = this_support_point(1) - obstacle_value;  // rather than obstacle_value - this_support_point(1);

                

                double e_modulus;

                if (  cell->material_id() == 0 ) e_modulus = this->cytoplasm_model.youngmodulus();

                else if (  cell->material_id() == 1 ) e_modulus = this->nucleus_model.youngmodulus();

                else e_modulus = 1;

                if ( lambda.in_local_range(index_X2) && diag_mass_matrix_vector_relevant.in_local_range(index_X2) )

                {

                  

                  lambda(index_X2);

                  diag_mass_matrix_vector_relevant(index_X2);

                  this->hanging_node_constraints.is_constrained(index_X2);

                  

                  const double c = 100.0 * e_modulus;

                  if ((lambda(index_X2) / diag_mass_matrix_vector_relevant(index_X2)

                       +

                       c * ( - solution_here - undeformed_gap )   // I believe it should be this, rather than (solution_here - undeformed_gap)

                       > 0)

                      &&

                      !(this->hanging_node_constraints.is_constrained(index_X2)) )

                  {

                    this->contact_constraints.add_line( index_X2 );   // it was this->all_constraints.add_line(index_X2);

                    this->contact_constraints.set_inhomogeneity(index_X2, - undeformed_gap); // it was this->all_constraints.set_inhomogeneity // I believe it should be this, rather than undeformed_gap);

                    distributed_solution( index_X2 ) = - undeformed_gap;  // I believe it should be this, rather than undeformed_gap;

                    this->active_set.add_index( index_X2 );

                  }

                };

                

              }

            }

          }

      }

    }

  

  // At the end of this function, we exchange data between processors updating those ghost elements in the solution variable that have been written by other processors.

  

  distributed_solution.compress(VectorOperation::insert);

  this->accumulated_displacement = distributed_solution;

  

  // We then merge the constraints from hanging nodes into the ConstraintMatrix object that already contains the active set.

  

  this->contact_constraints.close();    // it was this->all_constraints.close();

  this->all_constraints.merge(this->hanging_node_constraints);

  this->all_constraints.merge(this->contact_constraints);

  this->pcout << "   Size of active set: " << Utilities::MPI::sum((this->active_set & this->locally_owned_dofs).n_elements(), this->mpi_communicator) << std::endl;

  

}

Alberto Salvadori
 Dipartimento di Ingegneria Civile, Architettura, Territorio, Ambiente e di Matematica (DICATAM)
 Universita` di Brescia, via Branze 43, 25123 Brescia
 Italy
 tel 030 3711239
 fax 030 3711312

e-mail:
 alberto....@unibs.it
web-pages:
 http://m4lab.unibs.it/faculty.html
 http://dicata.ing.unibs.it/salvadori

--
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+unsubscribe@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

Informativa sulla Privacy: http://www.unibs.it/node/8155

[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/

[Wolfgang Bangerth]

>> I can do that for you. Do I understand correctly that the only change
>> that needs to happen is to move
>> active_set_old = active_set;
>> to *after* the if-statement?
>
> Exactly.

Great, thanks for confirming, Joerg!

The patch is now here:
https://www.dealii.org/8.5.1/doxygen/deal.II/Tutorial.html

Best
W.

[Alberto Salvadori]

Thank you, Timo. Your remarks have been very useful.

It turned out that I made a mistake in the way the mesh was prepared,  specifically some hanging nodes were not properly dealt with.

This caused also a related issue, that I shared here some time ago (Nov. 25).

This leads to another question, that I take the opportunity to ask. Suppose that a run is too long for a time slot allocated on a large scale computer. In such a case, one wants to restart the computations

from a given time. To this aim, the history of the computations up to a certain time shall be stored. All data stored in cell->user_pointer(), data of the loads and mesh.

How can one store information on the latter, - as for the partition and the hanging nodes - correctly? I understand that saving the mesh in a "ucd" or alternative form may not be the right strategy.

Thank you very much.

Alberto

Alberto Salvadori
 Dipartimento di Ingegneria Civile, Architettura, Territorio, Ambiente e di Matematica (DICATAM)
 Universita` di Brescia, via Branze 43, 25123 Brescia
 Italy
 tel 030 3711239
 fax 030 3711312

e-mail:
 alberto....@unibs.it
web-pages:
 http://m4lab.unibs.it/faculty.html
 http://dicata.ing.unibs.it/salvadori

--
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+unsubscribe@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

[Timo Heister]

> In such a case, one wants to restart the computations

We have support for snapshot&resume in parallel with adaptive meshes
by using the SolutionTransfer class. See "Use for serialization" in
https://urldefense.proofpoint.com/v2/url?u=http-3A__www.dealii.org_developer_doxygen_deal.II_classparallel-5F1-5F1distributed-5F1-5F1SolutionTransfer.html&d=DwIBaQ&c=Ngd-ta5yRYsqeUsEDgxhcqsYYY1Xs5ogLxWPA_2Wlc4&r=4k7iKXbjGC8LfYxVJJXiaYVu6FRWmEjX38S7JmlS9Vw&m=88KR9zh1f_gaVgWHpZwwXXlaOCduKgZ1yzsuERztdmU&s=Q549G3WmW1DtS06Yq9GyH98gIg8Rl44Ezac7I1Fc1Tc&e=

You can also check out how we do things in ASPECT:
https://urldefense.proofpoint.com/v2/url?u=https-3A__github.com_geodynamics_aspect_blob_master_source_simulator_checkpoint-5Frestart.cc&d=DwIBaQ&c=Ngd-ta5yRYsqeUsEDgxhcqsYYY1Xs5ogLxWPA_2Wlc4&r=4k7iKXbjGC8LfYxVJJXiaYVu6FRWmEjX38S7JmlS9Vw&m=88KR9zh1f_gaVgWHpZwwXXlaOCduKgZ1yzsuERztdmU&s=cUF_2LwSAUqoivqPWcoBBBpepBLrCfq8fDSkWlrkBTY&e=