[petsc-users] Updating TS solution outside PETSc

[Jed Brown]

"Ellen M. Price" <ellen.price at cfa.harvard.edu> writes:

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

This can be implemented by creating algebraic constraints (which yields
a DAE) or by formulating the problem in terms of ghost cells.

Note that it is not required that a given one-sided difference equation
be satisfied at a boundary.  Indeed, if you enforce

  (u[1] - u[0])/h = 0

you'll see only first order accuracy.  The DAE formulation (even if you
use a high-order difference formula) can also lead to order reduction in
time, depending on the time integrator.  It's often preferable to use
ghost points

  u[-1] = u[1]

and then discretize your internal equations, e.g.,

  udot[0] + (-u[-1] + 2*u[0] - u[1])/h2 = 0

which often reduces to incorporating a transient term with the one-sided
algebraic formula up top

  udot[0] * h/2 - (u[1] - u[0])/h = 0

but this yields an ODE and retains second order accuracy.

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