Re: [deal.II] Re: imposing boundary conditions for a three-fields formulation

[Wolfgang Bangerth]

1 - Can I manage the exception in a way that in case CG does not converge I
can try GMRes rather than aborting? Is the usual try-catch instruction working?
2 - I see that the class PETScWrappers implements both BiCGStab and GMRes. Can
you pinpoint to an example?


Yes, this works. Here is an example where we switch from one preconditioner to 
another if the first is not good enough to lead to convergence:


https://github.com/geodynamics/aspect/blob/master/source/simulator/solver.cc#L745

Best
 W.

--

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

--
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: [deal.II] Re: imposing boundary conditions for a three-fields formulation

[Alberto Salvadori]

Thank you Daniel and Jean-Paul.

Indeed I am using PETSc wrappers and the columns manipulation is not
active. And indeed if I run a little less academic example, i.e. more than
a single element, CG converges. Just to know, in case of a single element
CG does not, likely this is not unexpected. The error message is quite
self-explanatory (note that the solution is actually zero...):

**

*Exception on processing: *


**

*An error occurred in line <151> of file <../source/lac/petsc_solver.cc> in
function*

*void dealii::PETScWrappers::SolverBase::solve(const
dealii::PETScWrappers::MatrixBase &, dealii::PETScWrappers::VectorBase &,
const dealii::PETScWrappers::VectorBase &, const
dealii::PETScWrappers::PreconditionerBase &)*

*The violated condition was: *

*false*

*The name and call sequence of the exception was:*

*SolverControl::NoConvergence (solver_control.last_step(),
solver_control.last_value())*

*Additional Information: *

*Iterative method reported convergence failure in step 0. The residual in
the last step was nan.*


*This error message can indicate that you have simply not allowed a
sufficiently large number of iterations for your iterative solver to
converge. This often happens when you increase the size of your problem. In
such cases, the last residual will likely still be very small, and you can
make the error go away by increasing the allowed number of iterations when
setting up the SolverControl object that determines the maximal number of
iterations you allow.*


*The other situation where this error may occur is when your matrix is not
invertible (e.g., your matrix has a null-space), or if you try to apply the
wrong solver to a matrix (e.g., using CG for a matrix that is not symmetric
or not positive definite). In these cases, the residual in the last
iteration is likely going to be large.*

**


*Aborting!*

**

*Program ended with exit code: 1*

I am again very grateful with deal.II community, this is a great plus and I
hope one day to contribute, once my expertise will allow me to. Final
questions to close this call:

1 - Can I manage the exception in a way that in case CG does not converge I
can try GMRes rather than aborting? Is the usual try-catch instruction
working?
2 - I see that the class PETScWrappers implements both BiCGStab and GMRes.
Can you pinpoint to an example?

Thanks
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.salvad...@unibs.it
web-pages:
 http://m4lab.unibs.it/faculty.html
 http://dicata.ing.unibs.it/salvadori

On Thu, Apr 20, 2017 at 5:20 AM, Daniel Arndt <
daniel.ar...@iwr.uni-heidelberg.de> wrote:

> Alberto,
>
> Now I understand your concern. The boolean flag that is passed as the
>> final parameter to MatrixTools::apply_boundary_values
>> 
>> determines whether column elimination is performed. You've set it to false,
>> so this is not done and the RHS is unaltered. So switching it to true (the
>> default) should preserve the symmetry of the system and make the necessary
>> adjustments to the RHS vector.
>>
> While this is true in general, for PETSc and Trilinos matrices only the
> case "eliminate_columns = false" is implemented [1].
> This means that your system matrix will not be symmetric if you are using
> MatrixTools::apply_boundary_values, but a CG will likely still work [2,
> Global elimination].
>
> Best,
> Daniel
>
> [1] https://www.dealii.org/developer/doxygen/deal.II/
> namespaceMatrixTools.html#a96d8a143c787e0dca41ce6168152dc04
> [2] https://www.dealii.org/developer/doxygen/deal.II/
> namespaceMatrixTools.html
>
> --
> 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.
>

-- 

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

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

Re: [deal.II] Re: imposing boundary conditions for a three-fields formulation

[Daniel Arndt]

Alberto,

Now I understand your concern. The boolean flag that is passed as the final 
> parameter to MatrixTools::apply_boundary_values 
> 
>  
> determines whether column elimination is performed. You've set it to false, 
> so this is not done and the RHS is unaltered. So switching it to true (the 
> default) should preserve the symmetry of the system and make the necessary 
> adjustments to the RHS vector.
>
While this is true in general, for PETSc and Trilinos matrices only the 
case "eliminate_columns = false" is implemented [1]. 
This means that your system matrix will not be symmetric if you are using 
MatrixTools::apply_boundary_values, but a CG will likely still work [2, 
Global elimination].

Best,
Daniel

[1] 
https://www.dealii.org/developer/doxygen/deal.II/namespaceMatrixTools.html#a96d8a143c787e0dca41ce6168152dc04
[2] 
https://www.dealii.org/developer/doxygen/deal.II/namespaceMatrixTools.html

-- 
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: [deal.II] Re: imposing boundary conditions for a three-fields formulation

[Jean-Paul Pelteret]

Dear Alberto,

Now I understand your concern. The boolean flag that is passed as the final 
parameter to MatrixTools::apply_boundary_values 

 
determines whether column elimination is performed. You've set it to false, 
so this is not done and the RHS is unaltered. So switching it to true (the 
default) should preserve the symmetry of the system and make the necessary 
adjustments to the RHS vector.

J-P

On Thursday, April 20, 2017 at 2:19:50 AM UTC+2, Alberto Salvadori wrote:
>
> Dear Jean-Paul,
>
> thanks for your reply. 
>
> I see your point. It was my understanding (wrong, I guess) that in 
> imposing Dirichlet BC deal.II would modify the rows but also the columns, 
> adjusting the rhs accordingly. In this way, the symmetry of the system 
> matrix would have been preserved. I see now that this is not happening, in 
> all cases, is it correct? Then my question is about the solver, can we 
> still use CG as a solver if the system matrix looses symmetry? Does deal.II 
> make some further elaboration?
>
> Thank you again
>
> 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.salvad...@unibs.it
> web-pages:
>  http://m4lab.unibs.it/faculty.html
>  http://dicata.ing.unibs.it/salvadori
>
> On Tue, Apr 18, 2017 at 11:10 AM, Jean-Paul Pelteret <> wrote:
>
>> Dear Alberto,
>>
>> The correctness of my answer depends on the numbering assigned to the 
>> degrees-of-freedom, which is something that you've not mentioned above. But 
>> from what I can discern the 9th row relates to the pressure test function, 
>> so K_{8,[0-7]} are the components from the coupling block K_{pu}, K_{88} is 
>> K_{pp} and K_{89} is K_{pJ}. Since these degree's of freedom are not 
>> constrained, this row should remain untouched after boundary values are 
>> applied. So, to me, nothing necessarily appears amiss on row 9 (as you've 
>> indicated).
>>
>> Do you agree with this?
>>
>> Kind regards,
>> Jean-Paul
>>
>> On Tuesday, April 18, 2017 at 12:32:44 PM UTC+2, Alberto Salvadori wrote:
>>>
>>> Hi, 
>>>
>>> sorry for asking your help again, that I really appreciated. 
>>> Today's issue is about imposing boundary conditions for a three-fields 
>>> formulation.
>>>
>>> In brief, I defined a 4 nodes square element, of type:
>>>
>>> fe(FE_Q(poly_degree), dim, // displacement
>>>
>>>FE_DGPMonomial(poly_degree - 1), 1, // pressure
>>>
>>>FE_DGPMonomial(poly_degree - 1), 1), // dilatation
>>>
>>> with  poly_degree = 1.
>>>
>>> While debugging on a single element, which has 10 dofs in 2D ( 4x2 
>>> displacements + 1 pressure + 1 dilatation), I aim at imposing Dirichlet 
>>> boundary conditions on displacements, whereas body forces are zero. 
>>>
>>> The cell matrix and the cell_rhs are the following.
>>>
>>> *Cell matrix = *
>>>
>>> *{  0.7037,   0.0278,  -0.2037,  -0.4722,  -0.1481,   0.4722,  -0.3519,  
>>> -0.0278,  -1.,   0.}, *
>>>
>>> *{  0.0278,   0.7037,   0.4722,  -0.1481,  -0.4722,  -0.2037,  -0.0278,  
>>> -0.3519,  -1.,   0.}, *
>>>
>>> *{ -0.2037,   0.4722,   0.7037,  -0.0278,  -0.3519,   0.0278,  -0.1481,  
>>> -0.4722,   1.,   0.}, *
>>>
>>> *{ -0.4722,  -0.1481,  -0.0278,   0.7037,   0.0278,  -0.3519,   0.4722,  
>>> -0.2037,  -1.,   0.}, *
>>>
>>> *{ -0.1481,  -0.4722,  -0.3519,   0.0278,   0.7037,  -0.0278,  -0.2037, 
>>>   0.4722,  -1.,   0.}, *
>>>
>>> *{  0.4722,  -0.2037,   0.0278,  -0.3519,  -0.0278,   0.7037,  -0.4722,  
>>> -0.1481,   1.,   0.}, *
>>>
>>> *{ -0.3519,  -0.0278,  -0.1481,   0.4722,  -0.2037,  -0.4722,   0.7037, 
>>>   0.0278,   1.,   0.}, *
>>>
>>> *{ -0.0278,  -0.3519,  -0.4722,  -0.2037,   0.4722,  -0.1481,   0.0278, 
>>>   0.7037,   1.,   0.}, *
>>>
>>> *{ -1.,  -1.,   1.,  -1.,  -1.,   1.,   1., 
>>>   1.,   0.,  -4.}, *
>>>
>>> *{  0.,   0.,   0.,   0.,   0.,   0.,   0., 
>>>   0.,  -4.,   6.6667}, *
>>>
>>> *Cell rhs =0.000e+00 0.000e+00 0.000e+00 0.000e+00 0.000e+00 0.000e+00 
>>> 0.000e+00 0.000e+00 0.000e+00 0.000e+00 *
>>>
>>> They seem to be correct. After calling
>>>
>>>   system_matrix.compress(VectorOperation::add);
>>>
>>>   system_rhs.compress(VectorOperation::add);
>>>
>>> the system matrix coincides with the cell matrix, the same for the rhs. 
>>> To impose BCs I followed the scheme of step-18, as follows:
>>>
>>>   std::map boundary_values;
>>>
>>>
>>>   ComponentMask displacements_mask = fe.component_mask (displacements);
>>>
>>>   VectorTools::interpolate_boundary_values (dof_handler,
>>>
>>> 0,
>>>
>>> 

Re: [deal.II] Re: imposing boundary conditions for a three-fields formulation

[Alberto Salvadori]

Dear Jean-Paul,

thanks for your reply.

I see your point. It was my understanding (wrong, I guess) that in imposing
Dirichlet BC deal.II would modify the rows but also the columns, adjusting
the rhs accordingly. In this way, the symmetry of the system matrix would
have been preserved. I see now that this is not happening, in all cases, is
it correct? Then my question is about the solver, can we still use CG as a
solver if the system matrix looses symmetry? Does deal.II make some further
elaboration?

Thank you again

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.salvad...@unibs.it
web-pages:
 http://m4lab.unibs.it/faculty.html
 http://dicata.ing.unibs.it/salvadori

On Tue, Apr 18, 2017 at 11:10 AM, Jean-Paul Pelteret 
wrote:

> Dear Alberto,
>
> The correctness of my answer depends on the numbering assigned to the
> degrees-of-freedom, which is something that you've not mentioned above. But
> from what I can discern the 9th row relates to the pressure test function,
> so K_{8,[0-7]} are the components from the coupling block K_{pu}, K_{88} is
> K_{pp} and K_{89} is K_{pJ}. Since these degree's of freedom are not
> constrained, this row should remain untouched after boundary values are
> applied. So, to me, nothing necessarily appears amiss on row 9 (as you've
> indicated).
>
> Do you agree with this?
>
> Kind regards,
> Jean-Paul
>
> On Tuesday, April 18, 2017 at 12:32:44 PM UTC+2, Alberto Salvadori wrote:
>>
>> Hi,
>>
>> sorry for asking your help again, that I really appreciated.
>> Today's issue is about imposing boundary conditions for a three-fields
>> formulation.
>>
>> In brief, I defined a 4 nodes square element, of type:
>>
>> fe(FE_Q(poly_degree), dim, // displacement
>>
>>FE_DGPMonomial(poly_degree - 1), 1, // pressure
>>
>>FE_DGPMonomial(poly_degree - 1), 1), // dilatation
>>
>> with  poly_degree = 1.
>>
>> While debugging on a single element, which has 10 dofs in 2D ( 4x2
>> displacements + 1 pressure + 1 dilatation), I aim at imposing Dirichlet
>> boundary conditions on displacements, whereas body forces are zero.
>>
>> The cell matrix and the cell_rhs are the following.
>>
>> *Cell matrix = *
>>
>> *{  0.7037,   0.0278,  -0.2037,  -0.4722,  -0.1481,   0.4722,  -0.3519,
>> -0.0278,  -1.,   0.}, *
>>
>> *{  0.0278,   0.7037,   0.4722,  -0.1481,  -0.4722,  -0.2037,  -0.0278,
>> -0.3519,  -1.,   0.}, *
>>
>> *{ -0.2037,   0.4722,   0.7037,  -0.0278,  -0.3519,   0.0278,  -0.1481,
>> -0.4722,   1.,   0.}, *
>>
>> *{ -0.4722,  -0.1481,  -0.0278,   0.7037,   0.0278,  -0.3519,   0.4722,
>> -0.2037,  -1.,   0.}, *
>>
>> *{ -0.1481,  -0.4722,  -0.3519,   0.0278,   0.7037,  -0.0278,  -0.2037,
>> 0.4722,  -1.,   0.}, *
>>
>> *{  0.4722,  -0.2037,   0.0278,  -0.3519,  -0.0278,   0.7037,  -0.4722,
>> -0.1481,   1.,   0.}, *
>>
>> *{ -0.3519,  -0.0278,  -0.1481,   0.4722,  -0.2037,  -0.4722,   0.7037,
>> 0.0278,   1.,   0.}, *
>>
>> *{ -0.0278,  -0.3519,  -0.4722,  -0.2037,   0.4722,  -0.1481,   0.0278,
>> 0.7037,   1.,   0.}, *
>>
>> *{ -1.,  -1.,   1.,  -1.,  -1.,   1.,   1.,
>> 1.,   0.,  -4.}, *
>>
>> *{  0.,   0.,   0.,   0.,   0.,   0.,   0.,
>> 0.,  -4.,   6.6667}, *
>>
>> *Cell rhs =0.000e+00 0.000e+00 0.000e+00 0.000e+00 0.000e+00 0.000e+00
>> 0.000e+00 0.000e+00 0.000e+00 0.000e+00 *
>>
>> They seem to be correct. After calling
>>
>>   system_matrix.compress(VectorOperation::add);
>>
>>   system_rhs.compress(VectorOperation::add);
>>
>> the system matrix coincides with the cell matrix, the same for the rhs.
>> To impose BCs I followed the scheme of step-18, as follows:
>>
>>   std::map boundary_values;
>>
>>
>>   ComponentMask displacements_mask = fe.component_mask (displacements);
>>
>>   VectorTools::interpolate_boundary_values (dof_handler,
>>
>> 0,
>>
>> IncrementalDirichletBoundary
>> Values( TIDM, NRit, GPDofs ),
>>
>> boundary_values,
>>
>> displacements_mask);
>>
>>
>>
>>   PETScWrappers::MPI::Vector tmp (locally_owned_dofs,mpi_communicator);
>>
>>   MatrixTools::apply_boundary_values (boundary_values,
>>
>>   system_matrix,
>>
>>   tmp,
>>
>>   system_rhs,
>>
>>   false);
>>
>>   incremental_displacement = tmp;
>>
>> where the class IncrementalDirichletBoundaryValues implements the
>> Dirichlet BC and GPDofs = 4.
>> The boundary values after VectorTools::interpolate_boundary_values