[petsc-users] 2 Dirichlet conditions for one Element in PetscFE

Yann Jobic yann.jobic at univ-amu.fr
Thu Nov 23 07:09:52 CST 2017


Hello,


Le 23/11/2017 à 13:45, Matthew Knepley a écrit :
> On Thu, Nov 23, 2017 at 3:39 AM, Yann Jobic <yann.jobic at univ-amu.fr 
> <mailto:yann.jobic at univ-amu.fr>> wrote:
>
>     Hello,
>
>     I checked out  the branch knepley/fix-plex-bc-multiple, but i now
>     have a strange problem.
>     I splited my labels as in ex69.c of snes. It may be the right way
>     to do it.
>     In petsc 3.8.2, i have the same behavior as before, the element
>     containing the face is called by PetscDSAddBoundary.
>     In the git version, PetscDSAddBoundary does not call my boundary
>     function at all.
>     The call :
>       ierr = PetscDSAddBoundary(prob, DM_BC_ESSENTIAL, "wallL",
>     "markerLeft",   0, Ncomp, components, (void (*)(void)) uIn,  1,
>     &id, ctx);CHKERRQ(ierr);
>
>     I checked that my labels are correct :
>       markerLeft: 1 strata with value/size (1 (11))
>       markerTop: 1 strata with value/size (1 (11))
>       markerRight: 1 strata with value/size (1 (11))
>       markerBottom: 1 strata with value/size (1 (11))
>
>     What am i doing wrong ?
>
>
> So you call AddBoundary() and then the boundary values are never 
> inserted?
Yes exactly.
> The call looks correct. Can you send me an example to check?

My code is simple. It's ex46.c from ts, but i removed the temporal 
contributions in order to debug, and add the boundary.

Thanks a lot for the help !


Yann

> Obviously this works for my simple examples in the repo. I can't see 
> by looking what might be wrong for you.
>
>   Thanks,
>
>     Matt
>
>     Thanks,
>
>     Regards,
>
>     Yann
>
>     Le 22/11/2017 à 18:51, Matthew Knepley a écrit :
>>     On Wed, Nov 22, 2017 at 12:39 PM, Yann Jobic
>>     <yann.jobic at univ-amu.fr <mailto:yann.jobic at univ-amu.fr>> wrote:
>>
>>         Hello,
>>
>>         I've found a strange behavior when looking into a bug for the
>>         pressure convergence of a simple Navier-Stokes problem using
>>         PetscFE.
>>
>>         I followed many examples for labeling boundary faces. I first
>>         use DMPlexMarkBoundaryFaces, (label=1 to the faces).
>>         I find those faces using DMGetStratumIS and searching 1 as it
>>         is the value of the marked boundary faces.
>>         Finally i use DMPlexLabelComplete over the new label.
>>         I then use :
>>           ierr = PetscDSAddBoundary(prob, DM_BC_ESSENTIAL, "in",
>>         "Faces", 0, Ncomp, components, (void (*)(void)) uIn, NWest,
>>         west, NULL);CHKERRQ(ierr);
>>         in order to impose a dirichlet condition for the faces
>>         labeled by the correct value (west=1, south=3,...).
>>
>>         However, the function "uIn()" is called in all the Elements
>>         containing the boundary faces, and thus impose the values at
>>         nodes that are not in the labeled faces.
>>         Is it a normal behavior ? I then have to test the position of
>>         the node calling uIn, in order to impose the good value.
>>         I have this problem for a Poiseuille flow, where at 2 corner
>>         Elements i have a zero velocity dirichlet condition (wall)
>>         and a In flow velocity one.
>>
>>
>>     I believe I have fixed this in knepley/fix-plex-bc-multiple which
>>     should be merged soon. Do you know how to merge that branch and try?
>>
>>       Thanks,
>>
>>          Matt
>>
>>         The pressure is then very high at the corner nodes of those 2
>>         Elements.
>>         Do you think my pressure problem comes from there ? (The
>>         velocity field is correct)
>>
>>         Many thanks,
>>
>>         Regards,
>>
>>         Yann
>>
>>         PS : i'm using those runtime options :
>>         -vel_petscspace_order 2 -pres_petscspace_order 1 \
>>         -ksp_type fgmres -pc_type fieldsplit -pc_fieldsplit_type
>>         schur -pc_fieldsplit_schur_fact_type full  \
>>         -fieldsplit_velocity_pc_type lu -fieldsplit_pressure_ksp_rtol
>>         1.0e-10 -fieldsplit_pressure_pc_type jacobi
>>
>>
>>         ---
>>         L'absence de virus dans ce courrier électronique a été
>>         vérifiée par le logiciel antivirus Avast.
>>         https://www.avast.com/antivirus <https://www.avast.com/antivirus>
>>
>>
>>
>>
>>     -- 
>>     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.caam.rice.edu/%7Emk51/>
>
>
>
>
>
> -- 
> 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.caam.rice.edu/%7Emk51/>

-- 
___________________________

Yann JOBIC
HPC engineer
IUSTI-CNRS UMR 7343 - Polytech Marseille
Technopôle de Château Gombert
5 rue Enrico Fermi
13453 Marseille cedex 13
Tel : (33) 4 91 10 69 43
Fax : (33) 4 91 10 69 69



---
L'absence de virus dans ce courrier électronique a été vérifiée par le logiciel antivirus Avast.
https://www.avast.com/antivirus
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20171123/99fc3adb/attachment-0001.html>
-------------- next part --------------
static char help[] = "Navier-Stokes problem in 2d and 3d with finite elements.\n\
We solve the Navier-Stokes in a rectangular\n\
domain, using a parallel unstructured mesh (DMPLEX) to discretize it.\n\
This example supports discretized auxiliary fields (Re) as well as\n\
multilevel nonlinear solvers.\n\n\n";

#include <petscdmplex.h>
#include <petscdmforest.h>
#include <petscsnes.h>
#include <petscts.h>
#include <petscds.h>

/*
  Navier-Stokes equation:

  u . grad u - \Delta u - grad p = f
  div u  = 0
*/

typedef struct {
  PetscInt          dim;
  PetscInt          cells[2];
  char              filename[2048];    /* The optional mesh file */
  PetscErrorCode (**init_zero)(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar *u, void *ctx);
} AppCtx;

#define REYN 1.0

PetscErrorCode init_u(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar *u, void *ctx)
{
  u[0] = 0.;
  u[1] = 0.;
  return 0;
}

PetscErrorCode init_p(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar *p, void *ctx)
{
  p[0] = 1.;
  return 0;
}

PetscErrorCode zero_vector(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar *u, void *ctx)
{
  PetscInt d;
  for (d = 0; d < dim; ++d) u[d] = 0.0;
  return 0;
}

PetscErrorCode constant_p(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar *p, void *ctx)
{
  *p = 1.0;
  return 0;
}

static void f0_bd_u_2d(PetscInt dim, PetscInt Nf, PetscInt NfAux,
                       const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
                       const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
                       PetscReal t, const PetscReal x[], const PetscReal n[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[])
{
  f0[0] = 1;
  f0[1] = 0;
}

static void f1_bd_u(PetscInt dim, PetscInt Nf, PetscInt NfAux,
                    const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
                    const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
                    PetscReal t, const PetscReal x[], const PetscReal n[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f1[])
{
  const PetscInt Ncomp = dim;
  PetscInt       comp, d;
  for (comp = 0; comp < Ncomp; ++comp) {
    for (d = 0; d < dim; ++d) {
      f1[comp*dim+d] = 0.0;
    }
  }
}


static void f0_u(PetscInt dim, PetscInt Nf, PetscInt NfAux,
                      const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
                      const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
                      PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[])
{
  const PetscInt  Ncomp = dim;
  PetscInt        c, d;

  for (c = 0; c < Ncomp; ++c) {
    for (d = 0; d < dim; ++d) {
      f0[c] += u[d] * u_x[c*dim+d];
    }
  }
}

static void f1_u(PetscInt dim, PetscInt Nf, PetscInt NfAux,
                 const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
                 const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
                 PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f1[])
{
  const PetscReal Re    = REYN;
  const PetscInt  Ncomp = dim;
  PetscInt        comp, d;

  for (comp = 0; comp < Ncomp; ++comp) {
    for (d = 0; d < dim; ++d) {
      f1[comp*dim+d] = 1.0/Re * u_x[comp*dim+d];
    }
    f1[comp*dim+comp] -= u[Ncomp];
  }
}

static void f0_p(PetscInt dim, PetscInt Nf, PetscInt NfAux,
                 const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
                 const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
                 PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[])
{
  PetscInt d;
  for (d = 0, f0[0] = 0.0; d < dim; ++d) f0[0] += u_x[d*dim+d];
}

static void f1_p(PetscInt dim, PetscInt Nf, PetscInt NfAux,
                 const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
                 const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
                 PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f1[])
{
  PetscInt d;
  for (d = 0; d < dim; ++d) f1[d] = 0.0;
}

/*
  (psi_i, u_j grad_j u_i) ==> (\psi_i, \phi_j grad_j u_i)
*/
static void g0_uu(PetscInt dim, PetscInt Nf, PetscInt NfAux,
                  const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
                  const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
                  PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g0[])
{
  PetscInt NcI = dim, NcJ = dim;
  PetscInt fc, gc;
  PetscInt d;

  for (d = 0; d < dim; ++d) {
    g0[d*dim+d] = 0;
  }

  for (fc = 0; fc < NcI; ++fc) {
    for (gc = 0; gc < NcJ; ++gc) {
      g0[fc*NcJ+gc] += u_x[fc*NcJ+gc];
    }
  }
}

/*
  (psi_i, u_j grad_j u_i) ==> (\psi_i, \u_j grad_j \phi_i)
*/
static void g1_uu(PetscInt dim, PetscInt Nf, PetscInt NfAux,
                  const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
                  const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
                  PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g1[])
{
  PetscInt NcI = dim;
  PetscInt NcJ = dim;
  PetscInt fc, gc, dg;
  for (fc = 0; fc < NcI; ++fc) {
    for (gc = 0; gc < NcJ; ++gc) {
      for (dg = 0; dg < dim; ++dg) {
        /* kronecker delta */
        if (fc == gc) {
          g1[(fc*NcJ+gc)*dim+dg] += u[dg];
        }
      }
    }
  }
}

/* < q, \nabla\cdot u >
   NcompI = 1, NcompJ = dim */
static void g1_pu(PetscInt dim, PetscInt Nf, PetscInt NfAux,
                  const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
                  const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
                  PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g1[])
{
  PetscInt d;
  for (d = 0; d < dim; ++d) g1[d*dim+d] = 1.0; /* \frac{\partial\phi^{u_d}}{\partial x_d} */
}

/* -< \nabla\cdot v, p >
    NcompI = dim, NcompJ = 1 */
static void g2_up(PetscInt dim, PetscInt Nf, PetscInt NfAux,
                  const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
                  const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
                  PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g2[])
{
  PetscInt d;
  for (d = 0; d < dim; ++d) g2[d*dim+d] = -1.0; /* \frac{\partial\psi^{u_d}}{\partial x_d} */
}

/* < \nabla v, \nabla u + {\nabla u}^T >
   This just gives \nabla u, give the perdiagonal for the transpose */
static void g3_uu(PetscInt dim, PetscInt Nf, PetscInt NfAux,
                  const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
                  const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
                  PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g3[])
{
  const PetscReal Re    = REYN;
  const PetscInt  Ncomp = dim;
  PetscInt        compI, d;

  for (compI = 0; compI < Ncomp; ++compI) {
    for (d = 0; d < dim; ++d) {
      g3[((compI*Ncomp+compI)*dim+d)*dim+d] = 1.0/Re;
    }
  }
}

static PetscErrorCode ProcessOptions(MPI_Comm comm, AppCtx *options)
{
  PetscInt       n;
  PetscBool      flg;
  PetscErrorCode ierr;

  PetscFunctionBeginUser;
  options->dim           = 2;
  options->cells[0]      = 2;
  options->cells[1]      = 2;
  options->filename[0]   = '\0';

  ierr = PetscOptionsBegin(comm, "", "Navier-Stokes Equation Options", "DMPLEX");CHKERRQ(ierr);
  n = 2;
  ierr = PetscOptionsIntArray("-cells", "The initial mesh division", "ex.c", options->cells, &n, NULL);CHKERRQ(ierr);
  ierr = PetscOptionsString("-f", "Mesh filename to read", "ex.c", options->filename, options->filename, sizeof(options->filename), &flg);CHKERRQ(ierr);
  ierr = PetscOptionsInt("-dim", "The topological mesh dimension", "ex46.c", options->dim, &options->dim, NULL);CHKERRQ(ierr);
  ierr = PetscOptionsEnd();
  PetscFunctionReturn(0);
}

/* pris sur ex56.c de snes */
/* "boundary" == 1 => WEST */
/* "boundary" == 2 => EST */
/* "boundary" == 3 => SUD */
/* "boundary" == 4 => NORD */
static PetscErrorCode MarkBoundaryFaces(DM dm)
{
  DMLabel         label,label2;
  PetscInt        dim;
  IS              is;
  PetscErrorCode ierr;
  
  PetscFunctionBeginUser;
  ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr);
  ierr = DMCreateLabel(dm, "boundary");CHKERRQ(ierr);
  ierr = DMGetLabel(dm, "boundary", &label);CHKERRQ(ierr);
  ierr = DMPlexMarkBoundaryFaces(dm, label);CHKERRQ(ierr);

  ierr = DMGetStratumIS(dm, "boundary", 1,  &is);CHKERRQ(ierr);
  ierr = DMCreateLabel(dm,"Faces");CHKERRQ(ierr);
  if (is) {
    PetscInt        d, f, Nf;
    const PetscInt *faces;
    PetscInt        csize;
    PetscSection    cs;
    Vec             coordinates ;
    DM              cdm;
    ISGetLocalSize(is, &Nf);
    ISGetIndices(is, &faces);
    DMGetCoordinatesLocal(dm, &coordinates);
    DMGetCoordinateDM(dm, &cdm);
    DMGetDefaultSection(cdm, &cs);
    for (f = 0; f < Nf; ++f) {
      PetscReal   faceCoord;
      PetscInt    b,v;
      PetscScalar *coords = NULL;
      PetscInt    Nv;
      DMPlexVecGetClosure(cdm, cs, coordinates, faces[f], &csize, &coords);
      Nv   = csize/dim; 
      for (d = 0; d < dim; ++d) {
        faceCoord = 0.0;
        for (v = 0; v < Nv; ++v) faceCoord += PetscRealPart(coords[v*dim+d]);
        faceCoord /= Nv;
        for (b = 0; b < 2; ++b) {
          if (PetscAbs(faceCoord - b*1.0) < PETSC_SMALL) {
            DMSetLabelValue(dm, "Faces", faces[f], d*2+b+1);
            PetscPrintf(PETSC_COMM_WORLD, "face %d facecoords : %g label : %d\n", f,(double)faceCoord,d*2+b+1);
          }
        }
      }
      DMPlexVecRestoreClosure(cdm, cs, coordinates, faces[f], &csize, &coords);
    }
    ISRestoreIndices(is, &faces);
  }
  ISDestroy(&is);
  DMGetLabel(dm, "Faces", &label2);
  DMPlexLabelComplete(dm, label2);

  /* Make split labels so that we can have corners in multiple labels */
  {
    const char *names[4] = {"markerBottom", "markerRight", "markerTop", "markerLeft"};
    PetscInt    ids[4]   = {3, 2, 4, 1}, Nf;
    DMLabel     label;
    IS          is;
    PetscInt    f;

    for (f = 0; f < 4; ++f) {
      ierr = DMGetStratumIS(dm, "Faces", ids[f],  &is);CHKERRQ(ierr);
      if (!is) continue;
      ierr = ISGetLocalSize(is, &Nf);CHKERRQ(ierr);
      ierr = DMCreateLabel(dm, names[f]);CHKERRQ(ierr);
      ierr = DMGetLabel(dm, names[f], &label);CHKERRQ(ierr);
      if (is) { 
        ierr = PetscPrintf(PETSC_COMM_WORLD, "Nombre de face %s : %d ou l'on insere 1\n",names[f],Nf);
        ierr = DMLabelInsertIS(label, is, 1);CHKERRQ(ierr);
      }
      ierr = ISDestroy(&is);CHKERRQ(ierr);
    }
  }
  PetscFunctionReturn(0);
}

static PetscErrorCode CreateMesh(MPI_Comm comm, DM *dm, AppCtx *ctx)
{
  DM             pdm = NULL;
  const PetscInt dim = ctx->dim;
  const char    *filename  = ctx->filename;
  PetscBool      hasLabel;
  size_t         len;
  PetscErrorCode ierr;

  PetscFunctionBeginUser;
  ierr = PetscStrlen(filename, &len);CHKERRQ(ierr);
  if (!len) {
/*  ierr = DMPlexCreateHexBoxMesh(comm, dim, ctx->cells, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, dm);CHKERRQ(ierr); */
    ierr = DMPlexCreateBoxMesh   (comm, dim, 0, ctx->cells, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, 1, dm);CHKERRQ(ierr);
    ierr = PetscObjectSetName((PetscObject) *dm, "Mesh");CHKERRQ(ierr);
  } else {
    ierr = DMPlexCreateFromFile(comm, filename, 1, dm);CHKERRQ(ierr);
    ierr = DMPlexSetRefinementUniform(*dm, PETSC_FALSE);CHKERRQ(ierr);
  }
 
  ierr = PetscObjectSetName((PetscObject) *dm, "Mesh");CHKERRQ(ierr);
  /* Distribute mesh over processes */
  ierr = DMPlexDistribute(*dm, 0, NULL, &pdm);CHKERRQ(ierr);
  if (pdm) {
    ierr = DMDestroy(dm);CHKERRQ(ierr);
    *dm  = pdm;
  }
  ierr = DMSetFromOptions(*dm);CHKERRQ(ierr);
  /* If no boundary marker exists, mark the whole boundary, after refine of dm_refine from line options */
  ierr = DMHasLabel(*dm, "boundary", &hasLabel);CHKERRQ(ierr);
  if (!hasLabel) {ierr = MarkBoundaryFaces(*dm);CHKERRQ(ierr);}

  ierr = DMViewFromOptions(*dm, NULL, "-dm_view");CHKERRQ(ierr);
  PetscFunctionReturn(0);
}

PetscErrorCode zero(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar *u, void *ctx)
{
  PetscPrintf(PETSC_COMM_WORLD, "UIN  : coord : [%g,%g]\n",x[0],x[1]);
  u[0] = 0.;
  u[1] = 0.;
  return 0;
}

PetscErrorCode uIn(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar *u, void *ctx)
{
  PetscPrintf(PETSC_COMM_WORLD, "UIN  : coord : [%g,%g]\n",x[0],x[1]);
  u[0] = 0.005;
  u[1] = 0.;
  return 0;
}

static PetscErrorCode SetupProblem(PetscDS prob, AppCtx *ctx)
{
  const PetscInt id  = 1;
  const PetscInt Ncomp = 2; /*  dim = 2 */
  const PetscInt components[] = {0,1,2};
  PetscErrorCode ierr;
  
  PetscFunctionBeginUser;
  ierr = PetscDSSetResidual(prob, 0, f0_u, f1_u);CHKERRQ(ierr);
  ierr = PetscDSSetResidual(prob, 1, f0_p, f1_p);CHKERRQ(ierr);
  ierr = PetscDSSetJacobian(prob, 0, 0, g0_uu, g1_uu, NULL,  g3_uu);CHKERRQ(ierr);
  ierr = PetscDSSetJacobian(prob, 0, 1, NULL, NULL, g2_up, NULL);CHKERRQ(ierr);
  ierr = PetscDSSetJacobian(prob, 1, 0, NULL, g1_pu, NULL,  NULL);CHKERRQ(ierr);
  /* Neumann */
  /*ierr = PetscDSSetBdResidual(prob, 0, f0_bd_u_2d, f1_bd_u);CHKERRQ(ierr);*/

  ctx->init_zero[0] = init_u;
  ctx->init_zero[1] = init_p;
 
  ierr = PetscDSAddBoundary(prob, DM_BC_ESSENTIAL, "wallB", "markerBottom", 0, Ncomp, components, (void (*)(void)) zero, 1, &id, ctx);CHKERRQ(ierr);
  ierr = PetscDSAddBoundary(prob, DM_BC_ESSENTIAL, "wallT", "markerTop",    0, Ncomp, components, (void (*)(void)) zero, 1, &id, ctx);CHKERRQ(ierr);
  ierr = PetscDSAddBoundary(prob, DM_BC_ESSENTIAL, "wallL", "markerLeft",   0, Ncomp, components, (void (*)(void)) uIn,  1, &id, ctx);CHKERRQ(ierr);

  /*ierr = PetscDSAddBoundary(prob, DM_BC_NATURAL, "out", "Faces", 0, Ncomp, components, NULL, NEst, est, NULL);CHKERRQ(ierr);*/
  PetscFunctionReturn(0);
}

static PetscErrorCode SetupDiscretization(DM dm, MatNullSpace *nullSpace, AppCtx *ctx)
{
  DM              cdm = dm;
  const PetscInt  dim = ctx->dim;
  PetscDS         prob;
  PetscFE         fe[2];
  PetscQuadrature q;
  PetscObject     pressure;
  PetscErrorCode  ierr;

  PetscFunctionBeginUser;
  /* Create finite element */
  ierr = PetscFECreateDefault(dm, dim, dim, 0, "vel_", PETSC_DEFAULT, &fe[0]);CHKERRQ(ierr);
  ierr = PetscObjectSetName((PetscObject) fe[0], "velocity");CHKERRQ(ierr);
  ierr = PetscFEGetQuadrature(fe[0], &q);CHKERRQ(ierr);
  ierr = PetscFECreateDefault(dm, dim, 1, 0, "pres_", PETSC_DEFAULT, &fe[1]);CHKERRQ(ierr);
  ierr = PetscFESetQuadrature(fe[1], q);CHKERRQ(ierr);
  ierr = PetscObjectSetName((PetscObject) fe[1], "pressure");CHKERRQ(ierr);
  /* Set discretization and boundary conditions for each mesh */
  ierr = DMGetDS(dm, &prob);CHKERRQ(ierr);
  ierr = PetscDSSetDiscretization(prob, 0, (PetscObject) fe[0]);CHKERRQ(ierr);
  ierr = PetscDSSetDiscretization(prob, 1, (PetscObject) fe[1]);CHKERRQ(ierr);
  ierr = SetupProblem(prob, ctx);CHKERRQ(ierr);


  ierr = DMGetField(cdm, 1, &pressure);CHKERRQ(ierr);
  ierr = MatNullSpaceCreate(PetscObjectComm(pressure), PETSC_TRUE, 0, NULL, nullSpace);CHKERRQ(ierr);
  ierr = PetscObjectCompose(pressure, "nullspace", (PetscObject) *nullSpace);CHKERRQ(ierr);


  ierr = PetscFEDestroy(&fe[0]);CHKERRQ(ierr);
  ierr = PetscFEDestroy(&fe[1]);CHKERRQ(ierr);
  PetscFunctionReturn(0);
}

PetscErrorCode CreatePressureNullSpace(DM dm, AppCtx *user, Vec *v, MatNullSpace *nullSpace)
{
  Vec              vec;
  PetscErrorCode (*funcs[2])(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar *u, void* ctx) = {zero_vector, constant_p};
  PetscErrorCode   ierr;

  PetscFunctionBeginUser;
  ierr = DMGetGlobalVector(dm, &vec);CHKERRQ(ierr);
  ierr = DMProjectFunction(dm, 0.0, funcs, NULL, INSERT_ALL_VALUES, vec);CHKERRQ(ierr);
  ierr = VecNormalize(vec, NULL);CHKERRQ(ierr);

  ierr = MatNullSpaceCreate(PetscObjectComm((PetscObject)dm), PETSC_FALSE, 1, &vec, nullSpace);CHKERRQ(ierr);
  if (v) {
    ierr = DMCreateGlobalVector(dm, v);CHKERRQ(ierr);
    ierr = VecCopy(vec, *v);CHKERRQ(ierr);
  }
  ierr = DMRestoreGlobalVector(dm, &vec);CHKERRQ(ierr);
  /* New style for field null spaces */
  {
    PetscObject  pressure;
    MatNullSpace nullSpacePres;

    ierr = DMGetField(dm, 1, &pressure);CHKERRQ(ierr);
    ierr = MatNullSpaceCreate(PetscObjectComm(pressure), PETSC_TRUE, 0, NULL, &nullSpacePres);CHKERRQ(ierr);
    ierr = PetscObjectCompose(pressure, "nullspace", (PetscObject) nullSpacePres);CHKERRQ(ierr);
    ierr = MatNullSpaceDestroy(&nullSpacePres);CHKERRQ(ierr);
  }
  PetscFunctionReturn(0);
}


int main(int argc, char **argv)
{
  AppCtx         ctx;
  DM             dm;
  SNES           snes;
  Mat            A,J;
  PetscInt       its;
  Vec            u, r;
  MatNullSpace   nullSpace;            /* necessary for pressure */  
  PetscErrorCode ierr;

  ierr = PetscInitialize(&argc, &argv, NULL, help);CHKERRQ(ierr);
  ierr = ProcessOptions(PETSC_COMM_WORLD, &ctx);CHKERRQ(ierr);
  ierr = CreateMesh(PETSC_COMM_WORLD, &dm, &ctx);CHKERRQ(ierr);
  ierr = DMSetApplicationContext(dm, &ctx);CHKERRQ(ierr);
  ierr = PetscMalloc1(2, &ctx.init_zero);CHKERRQ(ierr);
  ierr = SetupDiscretization(dm, &nullSpace, &ctx);CHKERRQ(ierr);
  ierr = DMPlexCreateClosureIndex(dm, NULL);CHKERRQ(ierr);


  ierr = DMCreateGlobalVector(dm, &u);CHKERRQ(ierr);
  ierr = VecDuplicate(u, &r);CHKERRQ(ierr);

  ierr = SNESCreate(PETSC_COMM_WORLD, &snes);CHKERRQ(ierr);
  ierr = SNESSetDM(snes, dm);CHKERRQ(ierr);
  ierr = DMSetApplicationContext(dm, &ctx);CHKERRQ(ierr);

  ierr = DMCreateMatrix(dm, &J);CHKERRQ(ierr);
  A = J;
/*  ierr = CreatePressureNullSpace(dm, &ctx, NULL, &nullSpace);CHKERRQ(ierr);*/
/*  ierr = MatSetNullSpace(A, nullSpace);CHKERRQ(ierr);*/

  ierr = DMPlexSetSNESLocalFEM(dm,&ctx,&ctx,&ctx);CHKERRQ(ierr);
  ierr = SNESSetJacobian(snes, A, J, NULL, NULL);CHKERRQ(ierr);
  ierr = SNESSetFromOptions(snes);CHKERRQ(ierr);

  ierr = DMProjectFunction(dm, 0.0, ctx.init_zero, NULL, INSERT_ALL_VALUES, u);CHKERRQ(ierr);
  ierr = VecViewFromOptions(u, NULL, "-sol_init");CHKERRQ(ierr);

  ierr = SNESSolve(snes, NULL, u);CHKERRQ(ierr);
  ierr = SNESGetIterationNumber(snes, &its);CHKERRQ(ierr);
  ierr = PetscPrintf(PETSC_COMM_WORLD, "Number of SNES iterations = %D\n", its);CHKERRQ(ierr);

  ierr = VecViewFromOptions(u, NULL, "-sol_vec_view");CHKERRQ(ierr);

  ierr = VecDestroy(&u);CHKERRQ(ierr);
  ierr = VecDestroy(&r);CHKERRQ(ierr);

  if (A != J) {ierr = MatDestroy(&A);CHKERRQ(ierr);}
  ierr = MatDestroy(&J);CHKERRQ(ierr);
  ierr = VecDestroy(&u);CHKERRQ(ierr);
  ierr = VecDestroy(&r);CHKERRQ(ierr);
  ierr = SNESDestroy(&snes);CHKERRQ(ierr);
  ierr = DMDestroy(&dm);CHKERRQ(ierr);
  ierr = PetscFree(ctx.init_zero);CHKERRQ(ierr);
  ierr = PetscFinalize();
  return ierr;
}

/*TEST

mpirun -np 1 ./ns_stat -cells 10,10 -vel_petscspace_order 2 -pres_petscspace_order 1 -ksp_type fgmres -pc_type fieldsplit -pc_fieldsplit_type schur -pc_fieldsplit_schur_fact_type full -fieldsplit_velocity_pc_type lu -fieldsplit_pressure_ksp_rtol 1.0e-10 -fieldsplit_pressure_pc_type jacobi -ksp_monitor_short -ksp_converged_reason -snes_monitor_short -snes_converged_reason -sol_vec_view vtk:sol_ns.vtu:VTK_VTU 
mpirun -np 1 ./ns_stat -cells 10,10 -vel_petscspace_order 2 -pres_petscspace_order 1 -ksp_type fgmres -ksp_gmres_restart 10 -ksp_rtol 1.0e-9 -pc_type fieldsplit -pc_fieldsplit_type schur -pc_fieldsplit_schur_factorization_type full -fieldsplit_pressure_ksp_rtol 1e-10 -fieldsplit_velocity_ksp_type gmres -fieldsplit_velocity_pc_type lu -fieldsplit_pressure_pc_type jacobi -snes_error_if_not_converged -ksp_error_if_not_converged -snes_view -sol_vec_view vtk:sol_ns.vtu:VTK_VTU
bidon :
mpirun -np 1 ./ns_stat -cells 10,10 -vel_petscspace_order 2 -pres_petscspace_order 1 -pc_type fieldsplit -pc_fieldsplit_type schur -pc_fieldsplit_schur_fact_type full -fieldsplit_velocity_pc_type lu -fieldsplit_pressure_ksp_rtol 1.0e-10 -fieldsplit_pressure_pc_type jacobi -ksp_monitor_short  -ksp_type fgmres 
mpirun -np 1 ./ns_stat -cells 10,10 -vel_petscspace_order 2 -pres_petscspace_order 1 -ksp_type fgmres -pc_type fieldsplit -pc_fieldsplit_type schur -pc_fieldsplit_schur_fact_type lower -fieldsplit_0_ksp_type gmres -fieldsplit_0_pc_type bjacobi -fieldsplit_1_pc_type jacobi -fieldsplit_1_inner_ksp_type preonly -fieldsplit_1_inner_pc_type jacobi -fieldsplit_1_upper_ksp_type preonly -fieldsplit_1_upper_pc_type jacobi -ksp_converged_reason -snes_monitor_short -snes_converged_reason -sol_vec_view vtk:sol_ns.vtu:VTK_VTU

TEST*/


More information about the petsc-users mailing list