"Ellen M. Price" <[email protected]> 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 <[email protected]> >>> 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 >>
