[petsc-users] How to set boundary conditions for a structured finite volume (FV) discretization when using TS/SNES?

[zhfreewill]

Hi,Barry,
Thanks for reply. I meant I cannot get the right result as the FD version
did, although the code did not report any error when running.

Indeed, there were mistakes and I fixed it in the code (I also attached
additional derivation and code this time, please see the attachment).
Still, I cannot get the right output so far.

Regards,
Qw

Barry Smith <bsmith at petsc.dev> 于2021年11月26日周五 下午11:31写道：

>
>   Could you be more specific in saying how it fails (can you share the
> actual prototype code?)
>
>   Note here you seem to be missing a /dx term
>
> *G[i] = (u[i+1] - u[i]) - q_A;           **//Neumann BC:   du/dx(x=0.0) =
> q_A;*
>
>
>    Barry
>
> On Nov 26, 2021, at 4:02 AM, zhfreewill <zhfreewill at gmail.com> wrote:
>
> Hi all,
> I have a basic question on how to apply B.C. for a cell-centered
> structured FV discretization when using the SNES and Ts object.
> For instance, for a 1-D transient diffusion equation
>
> *u_t = u_xx, 0 < x < 1*
>
> with IC:
>
> *u(x,t=0) = 0;*
>
> and Dirichlet BC:
>
> *u(x=0) = u_A = 0.0;*
> *u(x=1) = u_B = 1.0;*
>
> or Neumann BC:
>
> *du/dx(x=0.0) = q_A = 0.0;*
> *du/dx(x=1.0) = q_B = 1.0;*
>
>
> When using TS of PETSc, the equation can be organized in form of
>     *U_t = G(t,u)*
> where *G* is the right-hand-side (RHS) function regardless of the
> specific method of discretization.
> *(1) For FD discretization*
>
> *A               |---dx--|               B*
> *|-------|-------|-------|-------|-------| grid size = nx*
> *0              i-1      i      i+1     nx-1*
>
> the vertex-centered FD discretization of the domain is
>
> *x(i) = i*dx (i = 0,..., nx)), where dx = 1/nx*
>
> and the RHS function G of the equation is
>
> *G[i] = (u[i+1]-2*u[i]+u[i-1])/dx^2*
>
> There have been dozens of examples demonstrating way of BC setting with a
> finite difference medthod, e.g.,
>
> https://github.com/petsc/petsc/blob/386ebab930a846575609b49074a40c46e7f8ed75/src/ts/tutorials/ex26.c#L293-L343
> which I think I have understood how it works and the following lines of
> code evaluate the RHS function *G(t,u)* in PETSc's *TSSetRHSFunction()*
> routine
> /*------------------FD discretization-------------------------*/
> *for (i = xs; i < xs + xm; i++) {*
>
> *// left boundary*
> *if (i == 0) {*
>
> *G[i] = u[i] - u_A;                **// Dirichlet BC: u(x=0) = u_A;*
> *G[i] = (u[i] - u[i+1])/dx - q_A;  **// Neumann BC:   du/dx(x=0.0) = q_A;*
>
> *} *
> *// right boundary*
> *else if (i == mx - 1) {*
>
> *G[i] = u[i] - u_B;                **// Dirichlet BC: u(x=1.0) = u_B;*
> *G[i] = (u[i]-u[i-1])/dx - q_B;    **// Neumann BC:   du/dx(x=0.0) = q_A;*
>
> *}*
> *// internal*
> *else {*
>
> *G[i] = (u[i-1] - 2 * u[i] + u[i + 1])/dx^2; **//u_xx*
>
> *}*
>
> *}**/*------------------------------------------------------------*/*
>
> The BCs, e.g., at the left boudnary, can be straightforwardly set as
> follows:
> Dirichlet BC:
>     *G[0] = u[0] - u_A*
>     which indicates *G[0] = u[0] - u_A = 0* and Dirichlet BC *u[0] = u_A* is
> thus applied
> Neumann BC:
>     *G[0] = (u[0] - u[1])/dx - q_A*
>     which indicates *G[0] = (u[0] - u[1])/dx - q_A = 0* and Neumann BC *du/dx
> = (u[0] - u[1])/dx = q_A* is applied
> This is the way the above examples apply BCs.
>
> For a structured, cell-centered FVM discretization, however, I tried and
> found the BCs can not be applied in a similar way. I will explain how I
> tried this as follows.
> *(2) FV discretization*
>
> *A               |---dx--|               B*
> *|---o---|---o---|---o---|---o---|---o---| cell size = nx*
> *W   w   P   e   E           *
> *0      i-1      i      i+1     nx-1*
>
> (NOTE: the center cell is denoted as P and its neighbouring cells as W, E; cell
> faces are denoted as w, e; the left and right boundary faces are A and B,
> respectively.)
>
> The cell-centered FV discretization of the domain is
>
> *x(i) = dx/2 + i*dx (i = 0,..., nx-1)), where dx = 1/nx*
>
> the diffusion term u_xx can be reorganized by using divergence theorem
>
> *(u_E-u_P)   (u_P-u_W)*
> *--------- - ---------*
> *    dx         dx*
>
> and the RHS function G of the equation is
>
> *G[i] = (u_E-u_P)/dx - (u_P-u_W)/dx = (u[i+1]-2*u[i]+u[i-1])/dx*
>
> the key lines of code in RHS function G(t,u) in PETSc's
> *TSSetRHSFunction()* routine can be
>
> */*------------------FV discretization-------------------------*/*
> *for (i = xs; i < xs + xm; i++) {*
>
> *// left boundary*
> *if (i == 0) {*
>
> *G[i] = 3.0 * u[i] - u[i+1] - 2.0 * u_A; **//Dirichlet BC: u(x=0.0) =
> u_A, does not work*
> *G[i] = (u[i+1] - u[i]) - q_A;           **//Neumann BC:   du/dx(x=0.0) =
> q_A;*
>
> *}*
> *// right boundary*
> *else if (i == mx - 1) {*
>
> *G[i] = 3.0 * u[i] - u[i-1] - 2.0 * u_B; **//Dirichlet BC: u(x=1.0) =
> u_B, does not work*
> *G[i] = (u[i]-u[i-1])/dx - q_B;          **//Neumann BC:   du/dx(x=1.0) =
> q_B, does not work*
>
> *}*
> *// internal*
> *else {*
>
> *G[i] = (u[i-1] - 2*u[i] + u[i+1])/dx;   **//u_x*
>
> *}*
>
> *}**/*------------------------------------------------------------*/*
>
> I tried to apply BCs in a way similar to that of the FD discretization.
> For the Dirichlet BC on the left boundary, because of the cell-centered
> arrangement of u[i], the left BC is given as *u_A* at cell face A instead
> of cell center.
> So, I first applied a linear interpolation to get the cell-centered value
> of *u[0]* from the left boundary value *u_A* and neighbouring *u[1]*
>
>     *u[0] = (u[1]+2*u_A)/3.0 (equivalent to 3*u[0] - u[1] - 2*u_A = 0.0)*
>
> *    |-dx/2-|------dx-----|*
> *    |------o------|------o------| linear interpolation: u[0] =
> (u[1]+2*u_A)/3.0*
> *    A      P      e      E*
> *   u_A   u[0]          u[1]*
>
> Dirichlet BC:
>     *G[0] = 3.0*u[0] - u[1] - 2.0*u_A*
>     in which *u[0] = (u[1]+2*u_A)/3.0* is implicitly calculated and I
> expect the Dirichlet BC *u(x=0) = u_A* is applied
>
> For the Neumann BC, left BC is given as a gradient of *u*, *q_A*, at cell
> face A
> Neumann BC:
>     *G[i] = (u[1] - u[0])/dx - q_A*
>     in which Neumann BC *du/dx = (u[1] - u[0])/dx = q_A* is expected to
> be applied.
>
> However, this FV discretization failed to work. Is there anything wrong or
> could you give me a hint?
>
> Best,
> Qw
>
>
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20211127/4133bc33/attachment-0001.html>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: UnsteadyDiffusion1D.c
Type: application/octet-stream
Size: 6555 bytes
Desc: not available
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20211127/4133bc33/attachment-0001.obj>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: Derivation.zip
Type: application/x-zip-compressed
Size: 2607803 bytes
Desc: not available
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20211127/4133bc33/attachment-0001.bin>