# [petsc-users] Updating TS solution outside PETSc

Smith, Barry F. bsmith at mcs.anl.gov
Thu Dec 5 13:38:39 CST 2019

```   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
>>

```