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

[zhfreewill]

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