# [petsc-users] How to apply boundary conditions when using TS?

Matthew Knepley knepley at gmail.com
Tue Nov 16 05:11:28 CST 2021

```On Mon, Nov 15, 2021 at 10:24 PM zhfreewill <zhfreewill at gmail.com> wrote:

> Hi,
> I'm new to PETSc. I have been confused about how to apply boundary
> conditions (BCs) (e.g., the Dirichlet and the Neumann type) when using the
> time-stepping (TS) object of PETSc.
>
> I have read the example codes (
> https://petsc.org/release/src/ts/tutorials/index.html) and vaguely found
> that BCs seems to be set  as follows:
>
> 1. when evaluating a matrix of the right hand side function g(x) (e.g.,
> using a function like RHSMatrixHeat)
> 2. when evaluating a matrix of the Jacobian matrix of g'(x) (e.g., using a
> function like FormJacobian)
>
> For instance, in the ex1.4 (
> the beginning states
>
>  17: /*
> ------------------------------------------------------------------------
>  19:    This program solves the one-dimensional heat equation (also called
> the
>  20:    diffusion equation),
>  21:        u_t = u_xx,
>  22:    on the domain 0 <= x <= 1, with the boundary conditions
>  23:        u(t,0) = 0, u(t,1) = 0,
>  24:    and the initial condition
>  25:        u(0,x) = sin(6*pi*x) + 3*sin(2*pi*x).
>  26:    This is a linear, second-order, parabolic equation.
>
>  28:    We discretize the right-hand side using finite differences with
>  29:    uniform grid spacing h:
>  30:        u_xx = (u_{i+1} - 2u_{i} + u_{i-1})/(h^2)
> ...
>  47:
> ------------------------------------------------------------------------- */
>
>  and Dirichlet BCs, u(t,0) = 0 and u(t,1) = 0, are to set on the both side
> of the discretized domain, the corresponding codes can be found on the
> following lines in a function RHSMatrixHeat:
>
> /*
> ------------------------------------------------------------------------- */
> 463: RHSMatrixHeat - User-provided routine to compute the right-hand-side
> matrix for the heat equation. */
> 491: PetscErrorCode RHSMatrixHeat(TS ts,PetscReal t,Vec X,Mat AA,Mat
> BB,void *ctx)
> 492: {
> ...
> 505:   /* Set matrix rows corresponding to boundary data */
> 509:   if (mstart == 0) {  /* first processor only */
> 510:     v[0] = 1.0;
> 511:     MatSetValues(A,1,&mstart,1,&mstart,v,INSERT_VALUES);
> 512:     mstart++;
> 513:   }
> 515:   if (mend == appctx->m) { /* last processor only */
> 516:     mend--;
> 517:     v[0] = 1.0;
> 518:     MatSetValues(A,1,&mend,1,&mend,v,INSERT_VALUES);
> 519:   }
>
> 521:   /* Set matrix rows corresponding to interior data.  We construct
> the matrix one row at a time. */
> 525:   v[0] = sone; v[1] = stwo; v[2] = sone;
> 526:   for (i=mstart; i<mend; i++) {
> 527:     idx[0] = i-1; idx[1] = i; idx[2] = i+1;
> 528:     MatSetValues(A,1,&i,3,idx,v,INSERT_VALUES);
> 529:   }
> ...
> 550: }
> /*
> ------------------------------------------------------------------------- */
>
>
> My questions are:
>
> 1. How do these lines of code embody the Dirichlet BCs u(t,0) = 0 and
> u(t,1) = 0.
>

This is the finite difference version of BC. It is just moving the term for
a known value to the forcing. However, since
the BC is 0, then we just do not compute the term. If the BC were nonzero,
you would also see the term in the
residual. For example,

https://gitlab.com/petsc/petsc/-/blob/main/src/snes/tutorials/ex19.c#L271

> 2. Similarly, if I want to apply a Neumann type BCs in these lines, How
> can I do?
>

For finite differences, you usually just approximate the derivative as the
one-sided difference of the last two values,
and set this equal to the known value. For finite elements, Neumann
conditions are just an addition weak form
integrated over the boundary.

Thanks,

Matt

> Any suggestions or documents?
>
> Best,
> Qw
>

--
What most experimenters take for granted before they begin their
experiments is infinitely more interesting than any results to which their