[petsc-users] ts behavior question

[Smith, Barry F.]

> On Nov 12, 2019, at 2:48 PM, Gideon Simpson <gideon.simpson at gmail.com> wrote:
> 
> I think I'm almost with you.  In my code, I make a local copy of the vector (with DMGetLocalVector) and after calling GlobaltoLocal, I call DMDAVecGetArray on the local vector.  I use the array I obtain of this local copy in populating by right hand side function.  Is that consistent with your the approach that you guys recommend?

  When you need ghost points yes. 

>  If I were to do as you say and have a separate set of calls for populating the ghost points, where would this fit in the ts framework?

   You could do it after you after the DMDAVecGetArray() call. So it is in the same routine. It would just "separate" more clearly the "setting the ghost values" from the "competing the RHS function".


  Barry


>  Are are you saying this would be done at the beginning of the RHS function?
> 
> On Tue, Nov 12, 2019 at 3:41 PM Smith, Barry F. <bsmith at mcs.anl.gov> wrote:
> 
> 
> > On Nov 12, 2019, at 2:09 PM, Gideon Simpson <gideon.simpson at gmail.com> wrote:
> > 
> > So this might be a resolution/another question.  Part of the reason to use the da is that it provides you with ghost points.  If you're only accessing the dependent variables entries with DMDAVecGetArrayRead, then you can't modify the ghost points.  If you can't modify the ghost points here, where would you do so in the context of a problem with, for instance, time dependent boundary conditions?
> 
>    In that case, as I say below, you have a ghosted local copy and you can put whatever values you wish into those ghosted locations. That is, when using ghosted local vectors you don't need to use the Read() version.
> 
>   Barry
> 
>   Note: if I were writing the code I would open the ghosted local input vector as writeable to put in the ghost values. Close it and then separately open it again as Read() to use in compute the needed TS functions. This is certainly not necessary but it helps with code maintainability and to decreases the likelihood of bugs. You have one set of access where you are legitimately changing values thus should not use Read() and another where you should not be changing values and thus should use read().
> 
> 
> 
> 
> > 
> > On Tue, Nov 12, 2019 at 10:43 AM Smith, Barry F. <bsmith at mcs.anl.gov> wrote:
> > 
> >   For any vector you only read you should use the read version.
> > 
> >   Sometimes the vector may not be locked and hence the other routine can be used but that may change as we add more locks and improve the code. So best to do it right
> > 
> > > On Nov 12, 2019, at 9:26 AM, Gideon Simpson <gideon.simpson at gmail.com> wrote:
> > > 
> > > So, in principle, should we actually be using DMDAVecGetArrayRead in this context?  I seem to be able to get away with DMDAVecGetArray with all time steppers.
> > 
> >   I am not sure why DMDAVecGetArray would  work if VecGetArray did not work. Internally it calls VecGetArray() that will do the check. If you call it on local ghosted vectors it doesn't check if the vector is locked since the ghosted version is a copy of the true locked vector.
> > 
> >    Barry
> > 
> > > 
> > > On Tue, Nov 12, 2019 at 12:33 AM Smith, Barry F. <bsmith at mcs.anl.gov> wrote:
> > > 
> > > 
> > > > On Nov 11, 2019, at 7:00 PM, Gideon Simpson via petsc-users <petsc-users at mcs.anl.gov> wrote:
> > > > 
> > > > I noticed that when I am solving a problem with the ts and I am *not* using a da, if I want to use an implicit time stepping routine:
> > > > 1. I have to explicitly provide the Jacobian
> > > 
> > >    Yes
> > > 
> > > > 2. When I do provide the Jacobian, if I want to access the elements of x(t) to construct f(t,x), I need to use a const PetscScalar and a VecGetArrayRead to get it to work.
> > > 
> > >   Presumably you call VecGetArray() instead? 
> > > > 
> > > >  
> > > > 3.  My code works without declaring const when I'm using an explicit scheme.
> > > > 
> > > > In contrast, if I solve a problem using a da, my code works, I can use implicit schemes without having to provide the Jacobian, and I don't have to use const anywhere.
> > > 
> > >   The use with DMDA provides automatic routines for computing the needed Jacobians using finite differencing of your provided function and coloring of the Jacobian. This results in reasonably efficient computation of Jacobians that work in most  (almost all) cases.
> > > > 
> > > > Can someone clarify what is expected/preferred?
> > > 
> > >   You should always use VecGetArrayRead() for vectors you are accessing but NOT changing the values in. There is no reason not and it provides the potential for higher performance.
> > > 
> > >   The algebraic solvers have additional checks to prevent peopled from inadvertently changing the entries in x (which would produce bugs). Presumably this results in generating an error when you call VecGetArray(). At least some of the TS explicit calls do not have such checks. They could be added and should be added.  https://gitlab.com/petsc/petsc/issues/493
> > > 
> > >   Thanks for pointing out the inconsistency
> > > 
> > >   Barry
> > > 
> > > > 
> > > > -- 
> > > > gideon
> > > 
> > > 
> > > 
> > > -- 
> > > gideon
> > 
> > 
> > 
> > -- 
> > gideon
> 
> 
> 
> -- 
> gideon