Re: [deal.II] Choosing the appropriate parameters for TrilinosWrappers::PreconditionILU for GLS stabilized Navier-Stokes

2019-02-12 Thread Wolfgang Bangerth
On 2/12/19 6:01 AM, Bruno Blais wrote:
> 
> The documentation suggest an atol between 1e-2 and 1e-5 and a rtol of 
> the order of 1.01. However, what should my ilu_fill and ilut_drop be? I 
> understand that increasing the fill will increase the memory consumption 
> and that I should try to reach a compromise, but right now I find that 
> it's mostly the speed of my iterative solve that is the problem. Are 
> there any guidelines to orient my choice of parameters? Am I going with 
> the right type of pre-conditioner or am I in the wrong here?

Others may have more experience with concrete values, but yours is a 
problem that is amendable to just trying things out: run the program for 
ten values for each of these parameters and plot the run-time for each 
run. That will give you an idea for how the behavior of the program 
depends on these values.

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.


[deal.II] Choosing the appropriate parameters for TrilinosWrappers::PreconditionILU for GLS stabilized Navier-Stokes

2019-02-12 Thread Bruno Blais
Hello everyone,
I am currently solving a stabilized version of the Navier-Stokes equation 
using dealii and everything is working well when I use a direct solver.
I would like (obviously?) to extend it to use a FGMRES iterative solver to 
be able to tackle very large systems. I have ported those aspects to use 
Trilinos.

The weak form  is (for now we can neglect the temperature related term) :


Concretely, the additional term due to the stabilization affect the 
velocity block in a SUPG manner and the pressure block by adding an 
element-dependent stiffness matrix.
I have read that ILU and ILUT type of preconditioners are the best-suited  
for this type of systems, but I find myself unable to orient my choice of 
parameters.
When we look at the constructor for the AdditionalData of ILU and ILUT 
preconditioner there are the following parameters:
ILU :
  AdditionalData(const unsigned int ilu_fill = 0,
 const double   ilu_atol = 0.,
 const double   ilu_rtol = 1.,
 const unsigned int overlap  = 0);
ILUT:
  AdditionalData(const double   ilut_drop = 0.,
 const unsigned int ilut_fill = 0,
 const double   ilut_atol = 0.,
 const double   ilut_rtol = 1.,
 const unsigned int overlap   = 0);

The documentation suggest an atol between 1e-2 and 1e-5 and a rtol of the 
order of 1.01. However, what should my ilu_fill and ilut_drop be? I 
understand that increasing the fill will increase the memory consumption 
and that I should try to reach a compromise, but right now I find that it's 
mostly the speed of my iterative solve that is the problem. Are there any 
guidelines to orient my choice of parameters? Am I going with the right 
type of pre-conditioner or am I in the wrong here?
Thank you very much for your help.

-- 
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] different results when compute integration with adaptive mesh

2019-02-12 Thread chucui1016
Dear Jean-Paul,

Thank you very much for you answer! It helps me a lot. The attached is the 
implemention of your idea, and the result is:
Cycle 0:
   Number of active cells:   20
   Number of degrees of freedom: 89
u_square:  0.0305621  grad_u_square:  0.256542
Cycle 1:
   Number of active cells:   44
   Number of degrees of freedom: 209
u_square:  0.0425496  grad_u_square:  0.315907
Cycle 2:
   Number of active cells:   92
   Number of degrees of freedom: 449
u_square:  0.0489012  grad_u_square:  0.341164
Cycle 3:
   Number of active cells:   200
   Number of degrees of freedom: 921
u_square:  0.0521656  grad_u_square:  0.353416
Cycle 4:
   Number of active cells:   440
   Number of degrees of freedom: 2017
u_square:  0.0540783  grad_u_square:  0.36169
Cycle 5:
   Number of active cells:   956
   Number of degrees of freedom: 4425
u_square:  0.0544531  grad_u_square:  0.363144
Cycle 6:
   Number of active cells:   1916
   Number of degrees of freedom: 8993
u_square:  0.0552556  grad_u_square:  0.366492
Cycle 7:
   Number of active cells:   3860
   Number of degrees of freedom: 18353
u_square:  0.0555461  grad_u_square:  0.36768


Thank you very much! And also thank Prof.Bangerth and Prof.Arndt! And I 
will think about the constraints much more.

Best,
Chucui 


在 2019年2月12日星期二 UTC+8下午4:01:10,chucu...@gmail.com写道:
>
> Dear Jean-Paul,
>
> Thank you very much for your detailed answer! I will write it as you say 
> right now and report the results as soon as possible !
>
> Best,
> Chucui
>
> 在 2019年2月12日星期二 UTC+8下午3:26:34,Jean-Paul Pelteret写道:
>>
>> Dear Chucui, 
>>
>> Too add to all of the other comments, you can "sanity check” the result 
>> that you’re getting using this approach by performing the same calculations 
>> another way and comparing results. Specifically, you could easily calculate 
>> these norms by looping over all cells and integrating the quantities 
>> yourself. As pseudocode it would look like this: 
>> 1. Loop over all cells in DoFHandler 
>> 2. Reinit FEValues with the update_values and update_gradients flags 
>> enabled 
>> 3. Compute local solution and its gradient at the cell quadrature points 
>> using fe_values.get_solution_[values/gradients]() 
>> 4. Loop over all quadrature points 
>> 5. u_square +=  local_solution[q_point]*local_solution[q_point] * 
>> fe_values.JxW() ; grad_u_square +=  local_gradient[q_point]* 
>> local_gradient[q_point] * fe_values.JxW() 
>>
>> Maybe this approach could help you debug your issue. 
>> Best, 
>> Jean-Paul 
>>
>> > On 12 Feb 2019, at 06:49, Wolfgang Bangerth  
>> wrote: 
>> > 
>> > On 2/11/19 12:20 AM, chucu...@gmail.com wrote: 
>> >> Dear Prof. Bangerth: 
>> >> 
>> >> Thank you very much for your quick reply! 
>> >> 
>> >> When I apply hanging node constraints to my matrix, as the step-6 
>> >> says: https://www.dealii.org/developer/doxygen/deal.II/step_6.html 
>> >> I make 4 steps: 
>> >>   1.Create a Constraints Class: 
>> >> | 
>> >> ConstraintMatrixconstraints; 
>> >> | 
>> >> 
>> >>   2.Fill this object using the function 
>> >> DoFTools::make_hanging_node_constraints() to ensure continuity of the 
>> elements 
>> >> of the finite element space. 
>> >> | 
>> >> DoFTools::make_hanging_node_constraints (dof_handler, 
>> >>  constraints); 
>> >> | 
>> >> 
>> >>   3.Copy the local contributions to the matrix and right hand side 
>> into the 
>> >> global objects: 
>> >> | 
>> >> constraints.distribute_local_to_global (cell_matrix, 
>> >> cell_rhs, 
>> >> local_dof_indices, 
>> >> system_matrix, 
>> >> system_rhs); 
>> >> constraints.distribute_local_to_global (cell_volume_matrix, 
>> >> local_dof_indices, 
>> >> volume_matrix); 
>> >> constraints.distribute_local_to_global (cell_gradient_matrix, 
>> >> local_dof_indices, 
>> >> gradient_matrix); 
>> >> | 
>> >> 
>> >> 
>> >>   4.Make sure that the degrees of "freedom" located on hanging nodes 
>> get 
>> >> their correct (constrained) value: 
>> >> | 
>> >> constraints.distribute (solution); 
>> >> | 
>> > 
>> > Constraints are funny and sometimes require deep thought about what 
>> exactly 
>> > they mean. What happens if you don't apply constraints to the 
>> > cell_volume_matrix and cell_gradient_matrix -- i.e., you copy the 
>> elements 1:1 
>> > into the matrix, without the 'constraints' object? (Like in step-4.) 
>> > 
>> > Alternatively, you could do as you show above and after calling 
>> > 'constraint.distribute(solution)' you manually set all constrained 
>> degrees of 
>> > freedom in the solution vector to zero. That's not what you want to use 
>> f

Re: [deal.II] different results when compute integration with adaptive mesh

2019-02-12 Thread chucui1016
Dear Jean-Paul,

Thank you very much for your detailed answer! I will write it as you say 
right now and report the results as soon as possible !

Best,
Chucui

在 2019年2月12日星期二 UTC+8下午3:26:34,Jean-Paul Pelteret写道:
>
> Dear Chucui, 
>
> Too add to all of the other comments, you can "sanity check” the result 
> that you’re getting using this approach by performing the same calculations 
> another way and comparing results. Specifically, you could easily calculate 
> these norms by looping over all cells and integrating the quantities 
> yourself. As pseudocode it would look like this: 
> 1. Loop over all cells in DoFHandler 
> 2. Reinit FEValues with the update_values and update_gradients flags 
> enabled 
> 3. Compute local solution and its gradient at the cell quadrature points 
> using fe_values.get_solution_[values/gradients]() 
> 4. Loop over all quadrature points 
> 5. u_square +=  local_solution[q_point]*local_solution[q_point] * 
> fe_values.JxW() ; grad_u_square +=  local_gradient[q_point]* 
> local_gradient[q_point] * fe_values.JxW() 
>
> Maybe this approach could help you debug your issue. 
> Best, 
> Jean-Paul 
>
> > On 12 Feb 2019, at 06:49, Wolfgang Bangerth  > wrote: 
> > 
> > On 2/11/19 12:20 AM, chucu...@gmail.com  wrote: 
> >> Dear Prof. Bangerth: 
> >> 
> >> Thank you very much for your quick reply! 
> >> 
> >> When I apply hanging node constraints to my matrix, as the step-6 
> >> says: https://www.dealii.org/developer/doxygen/deal.II/step_6.html 
> >> I make 4 steps: 
> >>   1.Create a Constraints Class: 
> >> | 
> >> ConstraintMatrixconstraints; 
> >> | 
> >> 
> >>   2.Fill this object using the function 
> >> DoFTools::make_hanging_node_constraints() to ensure continuity of the 
> elements 
> >> of the finite element space. 
> >> | 
> >> DoFTools::make_hanging_node_constraints (dof_handler, 
> >>  constraints); 
> >> | 
> >> 
> >>   3.Copy the local contributions to the matrix and right hand side into 
> the 
> >> global objects: 
> >> | 
> >> constraints.distribute_local_to_global (cell_matrix, 
> >> cell_rhs, 
> >> local_dof_indices, 
> >> system_matrix, 
> >> system_rhs); 
> >> constraints.distribute_local_to_global (cell_volume_matrix, 
> >> local_dof_indices, 
> >> volume_matrix); 
> >> constraints.distribute_local_to_global (cell_gradient_matrix, 
> >> local_dof_indices, 
> >> gradient_matrix); 
> >> | 
> >> 
> >> 
> >>   4.Make sure that the degrees of "freedom" located on hanging nodes 
> get 
> >> their correct (constrained) value: 
> >> | 
> >> constraints.distribute (solution); 
> >> | 
> > 
> > Constraints are funny and sometimes require deep thought about what 
> exactly 
> > they mean. What happens if you don't apply constraints to the 
> > cell_volume_matrix and cell_gradient_matrix -- i.e., you copy the 
> elements 1:1 
> > into the matrix, without the 'constraints' object? (Like in step-4.) 
> > 
> > Alternatively, you could do as you show above and after calling 
> > 'constraint.distribute(solution)' you manually set all constrained 
> degrees of 
> > freedom in the solution vector to zero. That's not what you want to use 
> for 
> > visualization, error estimation, ... but it might work for the operation 
> > you're trying here. 
> > 
> > I'm pretty sure that both of these solutions should work individually 
> (but not 
> > in combination). But explaining why this is so would take a page or two, 
> I'm 
> > afraid. 
> > 
> > Best 
> > W. 
> > 
> > 
> > -- 
> >  
> > Wolfgang Bangerth  email: bang...@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+un...@googlegroups.com . 
> > For more options, visit https://groups.google.com/d/optout. 
>
>

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