Re: [deal.II] Solving a system defined up to a constant (stabilized Navier-Stokes pressure definition)
On Thursday, 7 March 2019 23:35:15 UTC-5, Wolfgang Bangerth wrote: > > On 3/7/19 7:45 AM, Bruno Blais wrote: > > > > However, my GMRES stops very quickly after a certain number of newton > > iteration with the following : AztecOO::Iterate error code -4: GMRES > > Hessenberg ill-conditioned > > This is surprising and suggests that something is wrong with the matrix, > not > the right hand side. At least that's what I think this probably means. Did > you > have the same error when you did not modify the right hand side? What > happens > if you just artificially create a right hand side that you know is in the > range, for example by picking a random vector and multiplying it by the > matrix? > I did not have the same error when I did not modify the RHS. Which I also find very surprising. I also found that I did not get the same error if I modify the RHS, but that I also modify the Newton correction (if I project my newton correction such that corr = corr - P*corr) I will look at what you suggested, that's an extremely good idea, and get back at you. Thanks! > > 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+unsubscr...@googlegroups.com. For more options, visit https://groups.google.com/d/optout.
Re: [deal.II] Solving a system defined up to a constant (stabilized Navier-Stokes pressure definition)
On 3/7/19 7:45 AM, Bruno Blais wrote: > > However, my GMRES stops very quickly after a certain number of newton > iteration with the following : AztecOO::Iterate error code -4: GMRES > Hessenberg ill-conditioned This is surprising and suggests that something is wrong with the matrix, not the right hand side. At least that's what I think this probably means. Did you have the same error when you did not modify the right hand side? What happens if you just artificially create a right hand side that you know is in the range, for example by picking a random vector and multiplying it by the matrix? 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] Solving a system defined up to a constant (stabilized Navier-Stokes pressure definition)
Dear Wolfgang,
Sorry for an additional question.
I have tried to implement by myself something similar to what is done in
ASPECT (I guess?, it is a bit harder for me to grasp ASPECT concepts
because of the introspection and etc.).
However, I seem to have failed miserably at one point.
My understanding is that I can get
A x = f- Pf
By:
- First I calculate my projection operator. It is the integral of the
pressure shape function over ALL elements associated with a DOF. I
calculate it by looping over the elements and integrating with a quadrature
(see code below)
- Then I assemble my system and my right-hand side as usual
- Before I solve my linear system, I modify my Right-hand-side so that it
does not include the constant mode. I litteraly calculate f= f- Pf
- I solve my linear system for this new RHS.
However, my GMRES stops very quickly after a certain number of newton
iteration with the following : AztecOO::Iterate error code -4: GMRES
Hessenberg ill-conditioned
Clearly I am doing something wrong / I understood something wrong regarding
this, because I know the rest of my code works.
In my system there is no block matrices and everything is lumped in the
same matrix. I havent re-ordered by components either.
Calculation of the projection:
template
void
GLSNavierStokesSolver::initializePressureRHSCorrection()
{
TimerOutput::Scope t(computing_timer, "pressure_RHS_projector");
pressure_shape_function_integrals=0;
QGauss quadrature_formula(degreeQuadrature_);
FEValues fe_values (fe,
quadrature_formula,
update_values |
update_quadrature_points |
update_JxW_values );
const unsigned intdofs_per_cell = fe.dofs_per_cell;
const unsigned intn_q_points= quadrature_formula.
size();
const FEValuesExtractors::Scalar pressure (dim);
Vector
local_pressure_shape_function_integrals(dofs_per_cell);
std::vector local_dof_indices (dofs_per_cell);
typename DoFHandler::active_cell_iterator
cell = dof_handler.begin_active(),
endc = dof_handler.end();
for (; cell!=endc; ++cell)
{
if (cell->is_locally_owned())
{
fe_values.reinit(cell);
local_pressure_shape_function_integrals=0;
for (unsigned int i=0; iget_dof_indices (local_dof_indices);
zero_constraints.distribute_local_to_global(
local_pressure_shape_function_integrals,
local_dof_indices,
pressure_shape_function_integrals);
}
}
pressure_shape_function_integrals.compress (VectorOperation::add);
}
Projection of the System RHS onto it's space without the constant:
template
void GLSNavierStokesSolver::correctRHSClosedSystem()
{
// calculate projection of system_rhs on pressure_shape_function_integrals
for (unsigned int dof_i=0 ; dof_i < pressure_shape_function_integrals.size
() ; ++dof_i)
pressure_projection[dof_i] = pressure_shape_function_integrals[dof_i]*
system_rhs[dof_i];
// Add it to RHS
system_rhs.add(-1.,pressure_projection);
}
Thank you again for everything,
This is greatly appreciated.
Bruno
On Tuesday, 5 March 2019 16:10:09 UTC-5, Wolfgang Bangerth wrote:
>
>
> > A quick question. I think I understand what is done in
> >
> https://github.com/geodynamics/aspect/blob/master/source/simulator/helper_functions.cc#L1041
>
> > Weirdfully, I found the "pickaxe" version easier to understand (Love the
> > comments by the way, this is awesome)
>
> :-) I suspect that Timo gets credit for that one!
>
>
> > A question I have is related to the _pressure_shape_function_integrals_
> member.
> > My understanding is that this is the integration over the element cellI
> of the
> > shape function associated with pressure.
>
> I believe it's actually the integral over all cells the shape function is
> associated with.
>
>
> > However, I have looked over the ASPECT documentation and I have not been
> able
> > to find where this is implemented. Is it a vector that is manually
> filled or
> > is there a helper function that can be used to automatically generate
> it?
>
> Yes, that required a bit of searching. This is filled here:
>
> https://github.com/geodynamics/aspect/blob/master/source/simulator/assemblers/stokes.cc#L601
>
>
> 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
Re: [deal.II] Solving a system defined up to a constant (stabilized Navier-Stokes pressure definition)
> A quick question. I think I understand what is done in > https://github.com/geodynamics/aspect/blob/master/source/simulator/helper_functions.cc#L1041 > Weirdfully, I found the "pickaxe" version easier to understand (Love the > comments by the way, this is awesome) :-) I suspect that Timo gets credit for that one! > A question I have is related to the _pressure_shape_function_integrals_ > member. > My understanding is that this is the integration over the element cellI of > the > shape function associated with pressure. I believe it's actually the integral over all cells the shape function is associated with. > However, I have looked over the ASPECT documentation and I have not been able > to find where this is implemented. Is it a vector that is manually filled or > is there a helper function that can be used to automatically generate it? Yes, that required a bit of searching. This is filled here: https://github.com/geodynamics/aspect/blob/master/source/simulator/assemblers/stokes.cc#L601 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] Solving a system defined up to a constant (stabilized Navier-Stokes pressure definition)
Dear Wolfgang, A quick question. I think I understand what is done in https://github.com/geodynamics/aspect/blob/master/source/simulator/helper_functions.cc#L1041 Weirdfully, I found the "pickaxe" version easier to understand (Love the comments by the way, this is awesome) A question I have is related to the *pressure_shape_function_integrals* member. My understanding is that this is the integration over the element cellI of the shape function associated with pressure. However, I have looked over the ASPECT documentation and I have not been able to find where this is implemented. Is it a vector that is manually filled or is there a helper function that can be used to automatically generate it? Thanks for everything Bruno -- 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] Solving a system defined up to a constant (stabilized Navier-Stokes pressure definition)
Dear Bruno and Wolfgang, Thank you for your answers. I believe Wolfgang's answer is exactly what i had in mind (but said in clear words...). I will look at the Aspect code and try to port that to mine. Thank you for the very detailed answer. Best Bruno On Tuesday, 5 March 2019 10:08:29 UTC-5, Wolfgang Bangerth wrote: > > > > One of the issue of my system is that the pressure is defined up to a > > constant. On coarse mesh this does not affect the GMRES solver. However, > on > > finer mesh, it seems that the GMRES Solver is greatly affected by this > > near-singularity of the matrix system. > > I have often read in the literature that for stabilized method, the best > way > > was to remove the "mode" associated to a constant pressure. I believe > this > > implies some sort of projection of the residual in a space without a > pressure > > constant? > There are two parts of this problem: > > 1/ A deep theorem in linear algebra states that because the constant > pressures > are in the kernel of the matrix (i.e., Ay=0 for all vectors y that > correspond > to a constant pressure and zero velocity), that the *range* of the matrix > A > has dimension at most n-1. As a consequence, if you have a linear system >A x = f > then it is only possible to find a vector x with >|| A x - f || = 0 > if f is in the range of the matrix A. That will not generally be the case, > due > to discretization and integration errors. If f is not in the range of A, > there > is no way for GMRES to reduce the residual below a certain threshold, no > matter how many iterations one runs. > > The way to solve this is to solve >A x = f - Pf > instead where Pf is the projection onto the subspace not reachable by A. > This > is easy enough if you have a uniform mesh, but it's a bit complicated if > that's not the case. It's also complicated if one uses an enriched > pressure > space. Here is ASPECT's implementation of this step: > > https://github.com/geodynamics/aspect/blob/master/source/simulator/helper_functions.cc#L1041 > > You will be able to copy this and simplify it for your case. > > > 2/ There is now a null space and GMRES will find one of the infinitely > many > solutions of >A x = f > This may not be the solution you are looking for, so you probably want to > update the additive constant in the pressure after solution. How you want > to > do this depends on the application. Here again is the implementation in > ASPECT: > > https://github.com/geodynamics/aspect/blob/master/source/simulator/helper_functions.cc#L753 > > > I hope this helps! 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+unsubscr...@googlegroups.com. For more options, visit https://groups.google.com/d/optout.
Re: [deal.II] Solving a system defined up to a constant (stabilized Navier-Stokes pressure definition)
> One of the issue of my system is that the pressure is defined up to a > constant. On coarse mesh this does not affect the GMRES solver. However, on > finer mesh, it seems that the GMRES Solver is greatly affected by this > near-singularity of the matrix system. > I have often read in the literature that for stabilized method, the best way > was to remove the "mode" associated to a constant pressure. I believe this > implies some sort of projection of the residual in a space without a pressure > constant? There are two parts of this problem: 1/ A deep theorem in linear algebra states that because the constant pressures are in the kernel of the matrix (i.e., Ay=0 for all vectors y that correspond to a constant pressure and zero velocity), that the *range* of the matrix A has dimension at most n-1. As a consequence, if you have a linear system A x = f then it is only possible to find a vector x with || A x - f || = 0 if f is in the range of the matrix A. That will not generally be the case, due to discretization and integration errors. If f is not in the range of A, there is no way for GMRES to reduce the residual below a certain threshold, no matter how many iterations one runs. The way to solve this is to solve A x = f - Pf instead where Pf is the projection onto the subspace not reachable by A. This is easy enough if you have a uniform mesh, but it's a bit complicated if that's not the case. It's also complicated if one uses an enriched pressure space. Here is ASPECT's implementation of this step: https://github.com/geodynamics/aspect/blob/master/source/simulator/helper_functions.cc#L1041 You will be able to copy this and simplify it for your case. 2/ There is now a null space and GMRES will find one of the infinitely many solutions of A x = f This may not be the solution you are looking for, so you probably want to update the additive constant in the pressure after solution. How you want to do this depends on the application. Here again is the implementation in ASPECT: https://github.com/geodynamics/aspect/blob/master/source/simulator/helper_functions.cc#L753 I hope this helps! 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.
