[petsc-users] Using DM_BOUNDARY_GHOSTED in DMSTAG

[Barry Smith]

   The linear solvers (unless you write your own matrix-free beast, which I don't think you want to do) has no concept of the "ghost" locations and cannot access or manipulate them because the linear solver (and most standard preconditioners) work only with the global vector (not the local ghosted vector).

   This means you have one equation for each "point" on the DMStag or DMDA (for simplicity, assuming a scalar problem) and similarly one variable for each of those DM "points".  Resulting in a square matrix. 

   It sounds like you would like to have one equation for each point on the DMStag but variables on both the points on the DMStag and the extra ghost points, resulting in a rectangular matrix. The rectangular matrix is nice because it has the same regular stencil on all of its rows; no rows connected to the boundary are missing part of the stencil.  We don't support linear solvers that can work with the PETSc parallel matrices that can directly work with this form.

   Since you know the values on the ghost points, you can eliminate them (in theory or practice) by updating the right-hand side for the points on the DM. This elimination results in a square matrix (which now has the annoying boundary rows), which can then be solved.  This is the model that we work with for linear problems, having unknowns and equations only for variables in the global vector.

   Now, for nonlinear problems, it is a different story; here, using ghost points is very useful in evaluating f(x) (what becomes the right-hand side in Newton's method) and the Jacobian J(x). So both the function evaluation and Jacobian evaluation start by scattering the global x into local x (the ghost points between processes), the other ghost points (for boundary conditions) are filled by us as appropriate, and then a local function evaluation is done filling in the local points of the global vector using the values in the local vectors. Many of our SNES examples use this style.   

   If you do want to use your scheme for a linear problem directly you can do it. Do not have DM ghost boundary locations; instead, increase the size of the domain by one stencil width on each side and put the ghost boundary locations in the global solution vector, and make identity equations for each ghost boundary locations in the linear system. Conceptually, you have your nice rectangular matrix embedded in a square matrix by just having the other rows be rows of the identity matrix. And put in the right-hand side for those equations the values of your ghost locations, 

  There are a bunch of ways of thinking about these issues if you get into it,

   Barry


> On Aug 14, 2024, at 6:36 PM, Colton Bryant <coltonbryant2021 at u.northwestern.edu> wrote:
> 
> Hello,
> 
> I'm trying to understand the use of DM_BOUNDARY_GHOSTED and am a little confused. Is there any way for the linear solver to access and manipulate the ghost point value during the solve? I currently have a code using DM_BOUNDARY_PERIODIC and at the periodic boundary I simply apply the same discretization as I do everywhere else and as I understand it the value at e.g. i=-1 is set automatically by the periodic boundary condition. I would like to use DM_BOUNDARY_GHOSTED to set my own condition by which the point at i=-1 is set (a Neumann type condition). I have seen some matrix free examples but is there an easy way to "add" such a condition to the linear system in this case?
> 
> Thanks for any help you can provide.
> 
> Best,
> Colton Bryant