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 <zhfreew...@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
