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

[zhfreewill]

Thank you all! I think I did it. I used the ghost cells and the code
behaves properly. Attached is my document with detailed derivation and the
code.

Best,
Qw

Mark Adams <mfadams at lbl.gov> 于2021年11月28日周日 下午10:17写道：

> * I have used Matt's approach and it works but you do have to remember to
> update the ghost cells after you update your solution and before you
> use the operator again.
> * I would use a convergence test to determine if the result is correct.
> * Check that a constant function is a (trivial) solution to your
> discretization.
>
> On Sat, Nov 27, 2021 at 2:54 AM zhfreewill <zhfreewill at gmail.com> wrote:
>
>> 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/20211129/ad6abee7/attachment-0001.html>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: UnsteadyDiffusion1D_ver2.c
Type: application/octet-stream
Size: 6424 bytes
Desc: not available
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20211129/ad6abee7/attachment-0001.obj>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: Diffusion_FV_derivation.pdf
Type: application/pdf
Size: 2962554 bytes
Desc: not available
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20211129/ad6abee7/attachment-0001.pdf>