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
