[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