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.
[deal.II] Re: Solving a system defined up to a constant (stabilized Navier-Stokes pressure definition)
Bruno, On Tuesday, March 5, 2019 at 9:00:13 AM UTC-5, Bruno Blais wrote: > Are there any examples in DEALII where a Poisson problem is solved but > with strictly Neumann boundaries? That would be a very similar problem that > could guide me in the right direction. > Yes, take a look at step-11. I think one of the tutorial solving the Stokes equation also shows how to deal with that problem. If you can use it, method 1, i.e. just fix one dof, is the easiest. Best, 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.
[deal.II] Solving a system defined up to a constant (stabilized Navier-Stokes pressure definition)
Hello everyone, I am currently solving a GLS stabilized form of the Navier-Stokes equation using DEALII. The residual of the system looks similar to the regular Incompressible Navier-Stokes, except that a stiffness matrix that is dependent on the element size is added to the P-P block. I have a great success / scaling when solving Manufactured Solutions and I have been able to reach relatively decent P1-P! meshes (say 2056x2056 or 256x256x256) on a workstation computer. I am currently using the GMRES Trilinos Wrapper with ILU preconditioning with a relatively high fill-in level (4). 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? I know this is a very broad question, but how would I go and implement such a thing in my solver? Are there any examples in DEALII where a Poisson problem is solved but with strictly Neumann boundaries? That would be a very similar problem that could guide me in the right direction. Thank you very much for your help. The support of this forum is greatly appreciated. -- 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] Re: Conversion of the LDG-example from the code gallery for using MeshWorker
As initial debug approach, I used the function local_assemble_system instead of assemble_cell_terms, which worked (left the other functions alone). Then I tried to assemble just the boundary terms, by adding local_assemble_boundaries, uncommenting assemble_Dirichlet_boundary_terms() and assemble_Neumann_boundary_terms() and adding MeshWorker::AssembleFlags::assemble_boundary_faces to the flags in the mesh_loop. Here it failed, but I was not able to find out why yet. -- 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] Shared memory parallelism of Trilinos
Dear Bruno, > Thank you for the rapid and very detailed answer. This makes this community > so great. Bruno T. has given you far more insight into the inner workings of the Trilinos linear algebra packages than I’d have been able to (thanks :-) ), so I’m really glad that we were able to help you and we appreciate the positive feedback! Best, Jean-Paul -- 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.
