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

Sun Nov 28 08:17:41 CST 2021

* 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/20211128/d1a11154/attachment-0001.html>