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

[Barry Smith]

  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 <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/20211126/37069bae/attachment-0001.html>