[petsc-users] How to set boundary conditions for a structured finite volume (FV) discretization when using TS/SNES?
Matthew Knepley
knepley at gmail.com
Fri Nov 26 14:04:31 CST 2021
What we do in the FVM examples is to use ghost cells around the boundary.
Then the state in these cells is set so that the flux on the
boundary surface is exactly what you want. You could also prescribe the
fluxes directly, but that is not as natural in the setup I have used
for FVM.
Thanks,
Matt
On Fri, Nov 26, 2021 at 10:31 AM Barry Smith <bsmith at petsc.dev> wrote:
>
> 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
> 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
>
>
>
--
What most experimenters take for granted before they begin their
experiments is infinitely more interesting than any results to which their
experiments lead.
-- Norbert Wiener
https://www.cse.buffalo.edu/~knepley/ <http://www.cse.buffalo.edu/~knepley/>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20211126/94ba43de/attachment-0001.html>
More information about the petsc-users
mailing list