Little correction: I wrote "In the vmult() function remove the mean value (i.e., project the rhs on the orthogonal complement of the kernel of the kernel)", I meant "In the vmult() function remove the mean value after you multiply (i.e., project the rhs on the orthogonal complement of the kernel of the system matrix)"
This paper <https://epubs.siam.org/doi/abs/10.1137/S0036144503426074> by Bochev and Lehoucq describes the technique with a toy model (Poisson with purely natural BCs). It was not new at that time but to the best of my knowledge they were the first ones to write this down in a proper way. The technique I mentioned above works for other kernels in system matrices as well. For example in linearized elasticity one maybe needs to remove not the mean value but rigid motions (well not really rigid but infinitesimal rotations - invariance under rotations cannot be linear) from the kernel. Such motions do not affect stress but constitute of nontrivial displacements. Also in this case one can project the rhs onto the orthogonal complement of the kernel of the system matrix. In this case the kernel has dimension 3 in 2D (one for scaling+90 degree rotation and two independent shifts) and dimension 6 in 3D. Once you project the rhs a Krylov solver can deal with your singular problem. Cheers, Konrad On Saturday, December 26, 2020 at 11:06:36 AM UTC+1 Konrad Simon wrote: > Hi, > > On Saturday, December 26, 2020 at 6:43:21 AM UTC+1 smetca...@gmail.com > wrote: > >> Hi all, >> >> Does anyone here have any experience applying mean value constraints >> (specifically with periodic boundary conditions)? I'm having some trouble. >> As far as I can tell, there are two approaches (both of which I have been >> able to get working). The first approach is to just add the mean value >> constraint straight into the constraint matrix. I got this working just >> fine and was able to get an output but the problem is that this condition >> seems to make my matrices dense(?) or, at the very least, require a huge >> amount of memory hindering this approach for the types of fine spatial >> meshes I require for this application. >> > > If you think of the mean value constraint as an additional row it looks > like (for Poisson equation) u_0+u_1+...+u_n=0 and hence it couples all > DoFs. This is mathematically ok but not efficient since this constraint is > non-local. What one can also do is just constrain one DoF to a specific > value (this would also remove rigid motion in elasticity). But think about > your solution variable: If it is in the Sobolev space H^1 then point > evaluations may not be defined for dimension larger than 2. Similarly if, > for example, the pressure in a mechanical or fluid problem is often just in > L^2. Point evaluations do not make sense there at all. > > >> The other approach is to just apply periodic boundary conditions to the >> matrices, try to solve the system anyway and, because UMFPACK is awesome, >> somehow manage to get a solution out of all of this (which will of course >> only be defined up to a constant). The problem here is that the solution I >> am getting is of order 10^12! When I postprocess and subtract the mean >> value, what I get looks correct but I am losing a lot of accuracy in the >> process but it does work for the fine meshes that I require. Anyone have >> any ideas or do I just have to live with the reduced accuracy here? Thanks! >> > > This is of course like solving an ill-posed problem but some direct > solvers are smart enough to deal with zero pivots (redundant rows). > > Both approaches above do work but only for small and non-parallel > problems. This is what you can do: > > Many people do not think of the fact that even iterative solvers can deal > with singular problems if the right hand side is in the orthogonal > complement of the kernel. For the mean value constraint do the following: > > 1. Setup boundary conditions in constraints if necessary. > 2. Assemble the system matrix and rhs and apply your constraints > 3. Create a class that wraps your system matrix. This class must provide a > vmult() function > 4. Before starting the iterative solver remove the mean value, see step-57 > <https://www.dealii.org/current/doxygen/deal.II/step_56.html#StokesProblemprocess_solution> > > 4. In the vmult() function remove the mean value (i.e., project the rhs on > the orthogonal complement of the kernel of the kernel) > 5. use the wrapped matrix in the solver. > > Hope that helps. > > Best, > Konrad > -- 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. To view this discussion on the web visit https://groups.google.com/d/msgid/dealii/d013d489-b432-4551-bb71-b600fe8dd21bn%40googlegroups.com.
