[petsc-users] How to set boundary conditions for a structured finite volume (FV) discretization when using TS/SNES?
Barry Smith
bsmith at petsc.dev
Fri Nov 26 09:31:08 CST 2021
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>
More information about the petsc-users
mailing list