[petsc-users] Updating TS solution outside PETSc

[Smith, Barry F.]

   Let's look at the 1d case with vertex centered differencing 

                              |                  |     
      u-1*                u0             u1    ...


  u-1* is the ghost value, u0 the edge value etc. 

       Say you are solving        u_n on boundary is g(t),   u inside satisfies    u_t =  U_xx   

   Then if you provide a function of TSSetRHSFunction() it would be something like 

              DMGetLocalVector(DM,&ulocal);
              DMGlobalToLocal(dm, U, ulocal) ;

              DMDAGetArray(ulocal,&u)
              DMDAGetArray(fglobal,&f)

              for (i=0,....
                 f[i]  =  (u[i+1] - 2*u[i] + u[i-1])/h 


              f[0] += g(t)

               Note that f[0] =  (u[i+1] - 2*u[i] + u[i-1])/h^2 = (u[1] - u[0])/h + g(t)

               So you are satisfying the PDE inside the domain and on the boundary you are using one-sided differencing to approximate the Neumann boundary conditions.

               I may have some details wrong.

     Barry
 


         



> On Dec 5, 2019, at 11:30 AM, Ellen M. Price <ellen.price at cfa.harvard.edu> wrote:
> 
> I think I'm still unclear on exactly how this should work. Suppose, in
> my RHS function for TS, I'm processing the grid to compute its time
> derivative and get to an edge. What do I do?
> 
> If I set the derivative there to zero, the value will never change, but
> it *should* change so that the spatial derivative there is zero.
> 
> If I set it to the value it would get if it wasn't an edge, then the
> derivative isn't preserved anymore.
> 
> This is where I get stuck.
> 
> Ellen
> 
> 
> On 12/5/19 10:16 AM, Smith, Barry F. wrote:
>> 
>>   Are you using cell-centered or vertex centered discretization ( makes a slight difference)?
>> 
>>   Our model is to use DM_BOUNDARY_MIRROR  DMBoundaryType. This means that u_first_real_grid_point - u_its_ghost_point = 0 (since DMGlobalToLocal will automatically put into the physical ghost location the appropriate mirror values) thus u_n is zero along the edge; zero Neumann  conditions, for non-zero Neuman you need to put something in the "local rhs" to represent that, I'm not sure it is clear, but think about it in one dimension for the non-zero Neumann case.
>> 
>>   Bad news  not yet implemented for 3d.
>> 
>>   If you are using 3d we should fix this for you (or you fix it and make a MR because we should have this support). 
>> 
>>   If you have a 2d version of your code I would test with 3d your model etc and then let us know and request how we can get the 3d mirror written.
>> 
>>   Others may have alternative advice for Neumann with DMDA,
>> 
>>   Barry
>> 
>> 
>>> On Dec 5, 2019, at 10:00 AM, Ellen M. Price <ellen.price at cfa.harvard.edu> wrote:
>>> 
>>> Hi PETSc users,
>>> 
>>> I am working with a code that solves a set of PDEs on a rectangular
>>> domain with Neumann boundary conditions. My understanding of
>>> implementing the boundary condition is that I should set the boundary
>>> value to be the value that makes the finite difference derivative go to
>>> zero (or some other prescribed value) on that boundary.
>>> 
>>> I was attempting to update the solution from TS using a pre- or
>>> post-step/stage function and TSGetSolution, but this does not appear to
>>> work as expected. What would be the correct way to prescribe the
>>> boundary condition, given that I'm using TSRK for timestepping and a
>>> DMDA for discretization?
>>> 
>>> Looking forward to any help!
>>> 
>>> Ellen Price
>>