[petsc-users] Understanding PETSC: boundary layer flow with SNES

[Barry Smith]

  Sorry no one answered sooner.

> On Nov 17, 2015, at 9:06 AM, Pavel Schor <schor.pavel at gmail.com> wrote:
> 
> Dear PETSC users,
> I am a newbie to PETSC and I try to understand it. My first attempt is to solve 2D laminar boundary layer equations:
> 
>    u*du/dx +v*du/dy  = Ue*(dUe/dx). -nu*d2u/dy2
>    dudx + dvdy = 0
> 
> B.C.:
>    y=0:    u=0,    v=0
>    y=H:    u(x)=Ue(x)
> Where H is height of the domain, Ue(x) is edge velocity given by a potential flow solver. I assume Ue(x)=10.0
> 
> 
> I have several questions regarding the workflow in PETSC, I took SNES example29 (2D cavity flow) and I tried to modify it:
> 
> 1) There are vectors x and f. I suppose vector x to represent the values of velocity and vector f to represent velocity residuals. Please is this correct?
> 
> Based on example29 and my assumption, I calculated the velocity:
>  dhx = (PetscReal)(info->mx-1);  dhy = (PetscReal)(info->my-1);
>  dx = 1.0/dhx;                   dy = 1./dhy;
>    u          = x[j][i].u;
>    v          = x[j][i].v;
>    dudx        = (x[j][i+1].u - x[j][i-1].u)/2.0/dx;
>    dudy        = (x[j+1][i].u - x[j-1][i].u)/2.0/dy;
>    dvdy        = (x[j+1][i].v - x[j-1][i].v)/2.0/dy;
>    d2udy2        = (-2.0*u + x[j-1][i].u + x[j+1][i].u)/dy/dy;
> 
>    /* U velocity */
>    f[j][i].u  = u*dudx +v*dudy -Ue*dUedx -nu*d2udy2;
> 
>    /* V velocity */
>    f[j][i].v  = dudx + dvdy;
> 
> The code does not work. With initial conditions x[j][i].u=Ue; and x[0][i].u=0.0; the result is the same as the initial conditions.

  I haven't look in detailed but your general approach seems ok. You can run for a very small grid and print out the input and out vectors to see why they are not changing.

> 
> 2) In SNES example29, there are boundary conditions specified on vector f? For example:
> 
>  /* Test whether we are on the top edge of the global array */
>  if (yinte == info->my) {
>    j = info->my - 1;
>    yinte = yinte - 1;
>    /* top edge */
>    for (i=info->xs; i<info->xs+info->xm; i++) {
>        f[j][i].u     = x[j][i].u - lid;
>        f[j][i].v     = x[j][i].v;
> 
> I don't understand the last two lines. I just deleted the conditions for left and right edges and replaced f[j][i].u     = x[j][i].u - lid; with f[j][i].u     = x[j][i].u - Ue;

  Basically on the boundaries instead of using finite differences to satisfy the PDE we use an equation to satisfy the boundary conditions so what you do seem reasonable.

> 
> 3) Please could someone explain the normalization on following lines? (Taken from example 29)
> 
>      /*
>         Define mesh intervals ratios for uniform grid.
> 
>         Note: FD formulae below are normalized by multiplying through by
>         local volume element (i.e. hx*hy) to obtain coefficients O(1) in two dimensions.
> 
>      */
>      dhx = (PetscReal)(info->mx-1);  dhy = (PetscReal) (info->my-1);
>      hx = 1.0/dhx;                      hy = 1.0/dhy;
>      hxdhy = hx*dhy;                 hydhx = hy*dhx;

  When using geometric multigrid one must be careful that the interpolation function used between levels is correct. The easiest way to insure this is to use the same scaling for equations as the one gets with the finite element method. If you are not using multigrid the scaling is not important. We always like to use finite element scaling with our examples so one can use multigrid to solve the linear systems.

  Barry

> 
> 
> Thanks in advance & Kind regards
> Pavel Schor
> PhD. student, Institute of aerospace engineering, Brno University of technology