[petsc-users] DMPlex with AMR
Yann Jobic
yann.jobic at univ-amu.fr
Sun Nov 5 12:45:27 CST 2017
Dear PETSc expert,
I first correctly solve the advection/diffusion equation, with the
advection/diffusion of a gaussian, in a rotate field.
I'm using finit element, with PetscFE and the velocity field is in an
auxillary one, in order to correctly set the residual and the jacobian.
Then, i would like to use p4est, with AMR. To do that, I re-use the
example : "petsc/examples/src/dm/impls/forest/examples/tests/ex2.c"
I create a "base" DM, which has everything : PetscFE, PetscDS, the
boundaries, ... Then, i correctly adapt the mesh, by using :
ierr = DMForestTemplate(base,comm,&postForest);CHKERRQ(ierr);
ierr = DMForestSetAdaptivityLabel(postForest,adaptLabel);CHKERRQ(ierr);
My problem is that nothing happens when i'm solving the system, i.e. my
initial solution does not move, or diffuse.
That would mean that i didn't correctly transfer the problem definition.
I tried to attach everything to the adapted mesh, without succes.
What is the correct way to transfer the problem definition (PetscFE and
PetscDS) from one DM to another ?
Many thanks in advance,
Regards,
Yann
-------------- next part --------------
static char help[] = "Advection/diffusion of a gaussian using finite elements.\n\
We solve the problem in a rectangular\n\
domain, using a parallel unstructured mesh (DMPLEX) to discretize it.\n\n\n\n";
#include <petscdmplex.h>
#include <petscdmforest.h>
#include <petscds.h>
#include <petscts.h>
#include "petscpc.h" /* pour le precondditionnement LU */
#include <petscsys.h>
typedef struct {
PetscLogEvent createMeshEvent;
/* Domain and mesh definition */
PetscInt dim; /* The topological mesh dimension */
PetscInt cells[2]; /* The initial domain division */
char filename[2048]; /* The optional mesh file */
PetscBool interpolate; /* Generate intermediate mesh elements */
PetscInt monitorStepOffset; /* Pour que le TS monitor ajoute bien step */
/* Gaussian parameters */
PetscReal pos_init[2]; /* Position initiale de la gaussienne */
PetscReal sigma; /* ecart type de la gaussienne */
PetscReal dx;
PetscScalar Diffu; /* Diffusivite */
/* Exact solution */
PetscErrorCode (**exactConc) (PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx);
PetscErrorCode (**bord_d) (PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx);
} AppCtx;
/* refine_level :
1: coarse
2: normal
3: finner
*/
static PetscErrorCode AddIdentityLabel(DM dm)
{
PetscInt cStart, cEnd, c, medium;
PetscErrorCode ierr;
PetscFunctionBegin;
medium = 2;
ierr = DMCreateLabel(dm, "refine_level");CHKERRQ(ierr);
ierr = DMForestGetCellChart(dm,&cStart,&cEnd);CHKERRQ(ierr);
for (c = cStart; c < cEnd; c++) {ierr = DMSetLabelValue(dm, "refine_level", c, medium);CHKERRQ(ierr);} /* 2 par defaut */
PetscFunctionReturn(0);
}
static PetscErrorCode CreateAdaptivityLabel(DM forest, DMLabel *adaptLabel, Vec u, PetscReal ValMax)
{
const PetscInt debug = 0;
DM plex;
DMLabel RefineLevelLabel;
PetscInt cStart, cEnd, c, Nb, niv1, niv2, niv3; /* niv1 : coarse niv3 : finner */
PetscErrorCode ierr;
Vec localX;
PetscSection section;
PetscReal ValRefRefine,ValRefCoarse;
PetscFunctionBegin;
ValRefRefine = ValMax/10.; niv1 = 1; niv2 = 2; niv3 = 3;
ValRefCoarse = ValMax/100.;
ierr = DMLabelCreate("adapt",adaptLabel);CHKERRQ(ierr);
ierr = DMLabelSetDefaultValue(*adaptLabel,DM_ADAPT_DETERMINE);CHKERRQ(ierr);
ierr = DMGetLabel(forest,"refine_level",&RefineLevelLabel);CHKERRQ(ierr);
ierr = DMConvert(forest, DMPLEX, &plex);CHKERRQ(ierr);
DMGetDefaultSection(plex, §ion);
DMGetLocalVector(plex, &localX);
DMGlobalToLocalBegin(plex, u, INSERT_VALUES, localX);
DMGlobalToLocalEnd(plex, u, INSERT_VALUES, localX);
DMPlexGetHeightStratum(plex, 0, &cStart, &cEnd);
ierr = PetscPrintf(PETSC_COMM_WORLD, "COUCOU : \n\n");CHKERRQ(ierr);
for (c = cStart; c < cEnd; ++c) {
PetscFE fe;
PetscScalar *x = NULL;
PetscReal valEle = 0.0;
DMGetField(forest, 0, (PetscObject *) &fe);
PetscFEGetDimension(fe, &Nb);
DMPlexVecGetClosure(plex, 0, localX, c, NULL, &x);
if (debug) {
char title[1024];
PetscSNPrintf(title, 1023, "Solution for Field %d", 0);
DMPrintCellVector(c, title, Nb, &x[0]);
}
{
PetscInt i;
valEle = 0;
for(i=0; i<Nb; i++) valEle += (double) x[i];
valEle /= Nb;
if (valEle > ValRefRefine) {
ierr = DMLabelSetValue(*adaptLabel,c,DM_ADAPT_REFINE);CHKERRQ(ierr);
ierr = DMSetLabelValue(plex, "refine_level", c, niv3);CHKERRQ(ierr);
}
if (valEle < ValRefCoarse) {
ierr = DMLabelSetValue(*adaptLabel,c,DM_ADAPT_COARSEN);CHKERRQ(ierr);
ierr = DMSetLabelValue(plex, "refine_level", c, niv1);CHKERRQ(ierr);
}
if (debug) {PetscPrintf(PETSC_COMM_SELF, " elem %d val %g\n", c, valEle);}
}
DMPlexVecRestoreClosure(plex, NULL, localX, c, NULL, &x);
}
DMRestoreLocalVector(plex, &localX);
ierr = DMDestroy(&plex);CHKERRQ(ierr);
return(0);
}
static PetscErrorCode dirichlet (PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar *u, void *ctx)
{
u[0] = 0.;
return 0;
}
static PetscErrorCode rotate_field(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *v, void *ctx)
{
v[0] = -4.*(x[1]-0.5);
v[1] = 4.*(x[0]-0.5);
return 0;
}
static PetscErrorCode exactSol(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx)
{
PetscReal posx,posy,xx,yy,x0,y0,sigma2,C0,Diffu;
AppCtx *user = (AppCtx *) ctx;
x0 = user->pos_init[0];
y0 = user->pos_init[1];
sigma2 = PetscPowReal(user->sigma,2);
Diffu = user->Diffu;
C0 = 2*sigma2/(2*sigma2 + 4*Diffu*time);
posx = x[0]-0.5;
posy = x[1]-0.5;
xx = posx*PetscCosReal(4.*time) + posy*PetscSinReal(4.*time) - x0;
yy = -posx*PetscSinReal(4.*time) + posy*PetscCosReal(4.*time) - y0;
u[0] = C0*exp(-(PetscPowReal(xx,2) + PetscPowReal(yy,2))/(2*sigma2 + 4*Diffu*time));
return 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[])
{
PetscInt d;
f0[0] = u_t[0];
for (d = 0; d < dim; ++d) f0[0] += a[d] * u_x[d];
}
/* gradU[comp*dim+d] = {u_x, u_y} or {u_x, u_y, u_z} */
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[])
{
PetscInt d;
for (d = 0; d < dim; ++d) f1[d] = constants[0]*u_x[d];
}
/* <psi, u_t> */
static void g0_ut(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[])
{
g0[0] = u_tShift*1.0;
}
/* <psi, v . grad u> */
static void g1_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, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g1[])
{
PetscInt d;
for (d = 0; d < dim; ++d) g1[d] = a[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[])
{
PetscInt d;
for (d = 0; d < dim; ++d) g3[d*dim+d] = constants[0]*1.0;
}
static PetscErrorCode ProcessOptions(MPI_Comm comm, AppCtx *options)
{
PetscInt n;
PetscBool flg;
PetscErrorCode ierr;
PetscFunctionBeginUser;
options->dim = 2;
options->cells[0] = 25;
options->cells[1] = 25;
options->filename[0] = '\0';
options->interpolate = PETSC_FALSE;
ierr = PetscOptionsBegin(comm, "", "Poisson Problem Options", "DMPLEX");CHKERRQ(ierr);
ierr = PetscOptionsInt("-dim", "The topological mesh dimension", "ex12.c", options->dim, &options->dim, NULL);CHKERRQ(ierr);
n = 2;
ierr = PetscOptionsIntArray("-cells", "The initial mesh division", "ex12.c", options->cells, &n, NULL);CHKERRQ(ierr);
ierr = PetscOptionsString("-f", "Mesh filename to read", "ex12.c", options->filename, options->filename, sizeof(options->filename), &flg);CHKERRQ(ierr);
ierr = PetscOptionsBool("-interpolate", "Generate intermediate mesh elements", "ex12.c", options->interpolate, &options->interpolate, NULL);CHKERRQ(ierr);
ierr = PetscOptionsEnd();
ierr = PetscLogEventRegister("CreateMesh", DM_CLASSID, &options->createMeshEvent);CHKERRQ(ierr);
PetscFunctionReturn(0);
}
static PetscErrorCode CreateBCLabel(DM dm, const char name[])
{
DMLabel label;
PetscErrorCode ierr;
PetscFunctionBeginUser;
ierr = DMCreateLabel(dm, name);CHKERRQ(ierr);
ierr = DMGetLabel(dm, name, &label);CHKERRQ(ierr);
ierr = DMPlexMarkBoundaryFaces(dm, label);CHKERRQ(ierr);
ierr = DMPlexLabelComplete(dm, label);CHKERRQ(ierr);
PetscFunctionReturn(0);
}
static PetscErrorCode CreateMesh(MPI_Comm comm, AppCtx *user, DM *dm)
{
PetscInt dim = user->dim;
const char *filename = user->filename;
PetscBool interpolate = user->interpolate;
PetscBool hasLabel;
size_t len;
PetscErrorCode ierr;
PetscFunctionBeginUser;
ierr = PetscLogEventBegin(user->createMeshEvent,0,0,0,0);CHKERRQ(ierr);
ierr = PetscStrlen(filename, &len);CHKERRQ(ierr);
if (!len) {
ierr = DMPlexCreateHexBoxMesh(comm, dim, user->cells, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, dm);CHKERRQ(ierr);
ierr = PetscObjectSetName((PetscObject) *dm, "Mesh");CHKERRQ(ierr);
} else {
ierr = DMPlexCreateFromFile(comm, filename, interpolate, dm);CHKERRQ(ierr);
ierr = DMPlexSetRefinementUniform(*dm, PETSC_FALSE);CHKERRQ(ierr);
}
/* If no boundary marker exists, mark the whole boundary */
ierr = DMHasLabel(*dm, "marker", &hasLabel);CHKERRQ(ierr);
if (!hasLabel) {ierr = CreateBCLabel(*dm, "marker");CHKERRQ(ierr);}
{
PetscPartitioner part;
DM distributedMesh = NULL;
/* Distribute mesh over processes */
ierr = DMPlexGetPartitioner(*dm,&part);CHKERRQ(ierr);
ierr = PetscPartitionerSetFromOptions(part);CHKERRQ(ierr);
ierr = DMPlexDistribute(*dm, 0, NULL, &distributedMesh);CHKERRQ(ierr);
if (distributedMesh) {
ierr = DMDestroy(dm);CHKERRQ(ierr);
*dm = distributedMesh;
}
}
{
DM dmConv;
ierr = DMConvert(*dm,DMP4EST,&dmConv);CHKERRQ(ierr);
if (dmConv) {
ierr = DMDestroy(dm);CHKERRQ(ierr);
*dm = dmConv;
}
}
ierr = DMLocalizeCoordinates(*dm);CHKERRQ(ierr); /* needed for periodic */
ierr = DMSetFromOptions(*dm);CHKERRQ(ierr);
ierr = DMViewFromOptions(*dm, NULL, "-dm_view");CHKERRQ(ierr);
ierr = PetscLogEventEnd(user->createMeshEvent,0,0,0,0);CHKERRQ(ierr);
PetscFunctionReturn(0);
}
static PetscErrorCode SetupProblem(PetscDS prob, AppCtx *user)
{
PetscErrorCode ierr;
const PetscInt id = 1;
PetscScalar myconst[1];
PetscFunctionBeginUser;
ierr = PetscDSSetResidual(prob, 0, f0_u, f1_u);CHKERRQ(ierr);
ierr = PetscDSSetJacobian(prob, 0, 0, g0_ut, g1_u, NULL, g3_uu);CHKERRQ(ierr);
user->exactConc[0] = exactSol;
user->bord_d[0] = exactSol;
/*user->bord_d[0] = dirichlet;*/
myconst[0] = user->Diffu;
ierr = PetscDSSetConstants(prob,1,myconst);
ierr = PetscDSAddBoundary(prob, DM_BC_ESSENTIAL, "wall", "marker", 0, 0, NULL, (void (*)(void)) user->bord_d[0], 1, &id, user);CHKERRQ(ierr);
PetscFunctionReturn(0);
}
static PetscErrorCode SetupVelocity(DM dm, DM dmAux, AppCtx *user)
{
PetscErrorCode (*funcs[1])(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar v[], void *ctx) = {rotate_field};
Vec v;
PetscErrorCode ierr;
PetscFunctionBeginUser;
ierr = DMCreateLocalVector(dmAux, &v);CHKERRQ(ierr);
ierr = DMProjectFunctionLocal(dmAux, 0.0, funcs, NULL, INSERT_ALL_VALUES, v);CHKERRQ(ierr);
ierr = PetscObjectCompose((PetscObject) dm, "A", (PetscObject) v);CHKERRQ(ierr);
ierr = VecDestroy(&v);CHKERRQ(ierr);
PetscFunctionReturn(0);
}
static PetscErrorCode SetupDiscretization(DM dm, AppCtx *user)
{
const PetscInt debug = 0;
DM cdm = dm;
const PetscInt dim = user->dim;
PetscInt order,Ncomp,npoints;
PetscFE fe, feAux;
PetscDS prob, probAux;
PetscQuadrature q;
PetscErrorCode ierr;
PetscFunctionBeginUser;
/* Create finite element */
ierr = PetscFECreateDefault(dm, dim, 1, 0, "conc_", PETSC_DEFAULT, &fe);CHKERRQ(ierr); /*dim 1 : c'est un champ scalaire */
ierr = PetscObjectSetName((PetscObject) fe, "conc");CHKERRQ(ierr);
/* Create velocity */
ierr = PetscFECreateDefault(dm, dim, dim, 0, "vel_", -1, &feAux);CHKERRQ(ierr);
ierr = PetscFEGetQuadrature(fe, &q);CHKERRQ(ierr);
ierr = PetscFESetQuadrature(feAux, q);CHKERRQ(ierr);
ierr = PetscDSCreate(PetscObjectComm((PetscObject) dm), &probAux);CHKERRQ(ierr);
ierr = PetscDSSetDiscretization(probAux, 0, (PetscObject) feAux);CHKERRQ(ierr);
/* Set discretization and boundary conditions for each mesh */
ierr = DMGetDS(dm, &prob);CHKERRQ(ierr);
ierr = PetscDSSetDiscretization(prob, 0, (PetscObject) fe);CHKERRQ(ierr);
ierr = SetupProblem(prob, user);CHKERRQ(ierr);
while (cdm) {
DM dmAux, coordDM;
PetscBool hasLabel;
ierr = DMSetDS(cdm, prob);CHKERRQ(ierr);
ierr = DMGetCoordinateDM(cdm, &coordDM);CHKERRQ(ierr);
ierr = DMClone(cdm, &dmAux);CHKERRQ(ierr);
ierr = DMSetCoordinateDM(dmAux, coordDM);CHKERRQ(ierr);
ierr = DMSetDS(dmAux, probAux);CHKERRQ(ierr);
ierr = PetscObjectCompose((PetscObject) cdm, "dmAux", (PetscObject) dmAux);CHKERRQ(ierr);
ierr = SetupVelocity(cdm, dmAux, user);CHKERRQ(ierr);
ierr = DMDestroy(&dmAux);CHKERRQ(ierr);
ierr = DMHasLabel(cdm, "marker", &hasLabel);CHKERRQ(ierr);
if (!hasLabel) {ierr = CreateBCLabel(cdm, "marker");CHKERRQ(ierr);}
ierr = DMGetCoarseDM(cdm, &cdm);CHKERRQ(ierr);
}
/* On test l ordre de quadrature */
if (debug) {
ierr = PetscQuadratureGetOrder(q,&order);CHKERRQ(ierr);
ierr = PetscPrintf(PETSC_COMM_WORLD, "Ordre de quadrature : %D\n",order);CHKERRQ(ierr);
ierr = PetscQuadratureGetData(q,NULL,&Ncomp,&npoints,NULL,NULL);CHKERRQ(ierr);
ierr = PetscPrintf(PETSC_COMM_WORLD, "Nombre de points d'integration : %D\n\n",npoints);CHKERRQ(ierr);
}
/*pas de condition limite pour le moment */
ierr = PetscFEDestroy(&fe);CHKERRQ(ierr);
ierr = PetscFEDestroy(&feAux);CHKERRQ(ierr);
ierr = PetscDSDestroy(&probAux);CHKERRQ(ierr);
PetscFunctionReturn(0);
}
static PetscErrorCode GetDiscretizationInfo(DM dm) {
DM cdm = dm;
PetscFE fe;
PetscQuadrature quad;
PetscInt field, order, Ncomp, npoints;
PetscErrorCode ierr;
field = 0;
ierr = DMGetField(cdm, field, (PetscObject *) &fe);CHKERRQ(ierr);
ierr = PetscFEGetQuadrature(fe, &quad);CHKERRQ(ierr);
ierr = PetscQuadratureGetOrder(quad,&order);CHKERRQ(ierr);
ierr = PetscPrintf(PETSC_COMM_WORLD, "Ordre de quadrature : %D\n",order);CHKERRQ(ierr);
ierr = PetscQuadratureGetData(quad,NULL,&Ncomp,&npoints,NULL,NULL);CHKERRQ(ierr);
ierr = PetscPrintf(PETSC_COMM_WORLD, "Nombre de points d'integration : %D\n\n",npoints);CHKERRQ(ierr);
PetscFunctionReturn(0);
}
static PetscErrorCode TSMonitorError(TS ts, PetscInt step, PetscReal crtime, Vec u, void *ctx)
{
AppCtx *user = (AppCtx *) ctx;
DM dm;
PetscReal error;
PetscErrorCode ierr;
void *myctx[1]; /* pour envoyer la structure dans DMProjectFunction */
PetscFunctionBeginUser;
if (step >= 0) {
step += user->monitorStepOffset;
}
myctx[0] = (void *) user;
if (step%10 == 0) {
ierr = TSGetDM(ts, &dm);CHKERRQ(ierr);
ierr = DMComputeL2Diff(dm, crtime, user->exactConc, myctx, u, &error);CHKERRQ(ierr);
ierr = PetscPrintf(PETSC_COMM_WORLD, "Timestep: %04d time = %-8.4g \t L_2 Error: %2.5g\n", (int) step, (double) crtime, (double) error);CHKERRQ(ierr);
}
PetscFunctionReturn(0);
}
static PetscErrorCode initializeTS(DM dm, AppCtx *ctx, TS *ts, PetscReal tmax, PetscReal dt)
{
PetscErrorCode ierr;
KSP ksp;
SNES ts_snes; /* nonlinear solver */
PC pc; /* preconditioner context */
PetscFunctionBegin;
/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
Create timestepping solver context
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
ierr = TSCreate(PetscObjectComm((PetscObject)dm), ts);CHKERRQ(ierr);
ierr = TSSetDM(*ts, dm);CHKERRQ(ierr);
/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
Customize nonlinear solver
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
ierr = TSMonitorSet(*ts, TSMonitorError, ctx, NULL);CHKERRQ(ierr);
ierr = DMTSSetBoundaryLocal(dm, DMPlexTSComputeBoundary, ctx);CHKERRQ(ierr);
ierr = DMTSSetIFunctionLocal(dm, DMPlexTSComputeIFunctionFEM, ctx);CHKERRQ(ierr);
ierr = DMTSSetIJacobianLocal(dm, DMPlexTSComputeIJacobianFEM, ctx);CHKERRQ(ierr);
ierr = TSSetExactFinalTime(*ts,TS_EXACTFINALTIME_MATCHSTEP);CHKERRQ(ierr);
ierr = TSSetMaxTime(*ts,tmax);CHKERRQ(ierr);
ierr = TSSetTimeStep(*ts,dt);CHKERRQ(ierr);
ierr = TSSetType(*ts,TSBEULER);CHKERRQ(ierr);
/*
TSGetType(TS,&time_scheme);
PCGetType(PC
TSThetaSetTheta(ts,1.0);
TSSetType(ts,TSCN);
*/
ierr = TSGetSNES(*ts,&ts_snes);CHKERRQ(ierr);
ierr = SNESGetKSP(ts_snes, &ksp);CHKERRQ(ierr);
ierr = KSPGetPC(ksp,&pc);CHKERRQ(ierr);
ierr = PCSetType(pc,"lu");CHKERRQ(ierr);
PetscFunctionReturn(0);
}
int main(int argc, char **argv)
{
MPI_Comm comm;
DM base, postForest; /* Problem specification */
TS ts; /* Pour l'iteration en temps */
PetscReal ftime;
PetscReal tmax,dt; /* temps max, delta t */
PetscInt steps; /* iterations for convergence */
PetscInt adaptInterval;/* TS step entre les remaillage */
Vec u,uNew,u_exact,r; /* solution, exact solution, residual vectors */
AppCtx user; /* user-defined work context */
PetscReal error = 0.0; /* L_2 error in the solution */
PetscReal Peclet,cfl,dx,Lx,Um;
PetscErrorCode ierr;
void *ctx[1]; /* pour envoyer la structure dans DMProjectFunction */
ctx[0] = (void *) &user;
TSType time_scheme;
DMLabel adaptLabel;
PetscBool hasLabel;
PetscDS prob;
ierr = PetscInitialize(&argc, &argv, NULL,help);if (ierr) return ierr;
ierr = ProcessOptions(PETSC_COMM_WORLD, &user);CHKERRQ(ierr);
comm = PETSC_COMM_WORLD;
Lx = 1;
Um = 4*Lx/2.;
dx = 1.*Lx/user.cells[0];
user.dx = dx;
tmax = 0.25*PETSC_PI/2.; /* rotation complete en PETSC_PI/2 */
/* dt = 0.000625; */
dt = 0.0025;
/* On init des parametres de la gaussienne */
user.pos_init[0] = 0.;
user.pos_init[1] = -0.125;
user.sigma = 3*dx;
user.Diffu = 0.02; /* Um*Lx/Diffu */
Peclet = Um*dx/user.Diffu;
cfl = Um*dt/dx;
ierr = PetscPrintf(comm, "Probleme : advection/diffusion\n cas test rotation de la gaussienne\n");CHKERRQ(ierr);
ierr = PetscPrintf(comm, " cas test de la gaussienne en rotation \n");CHKERRQ(ierr);
ierr = PetscPrintf(comm, "Parametres simu : \n");CHKERRQ(ierr);
ierr = PetscPrintf(comm, " CellsX : %D Lx : %g dx : %g dt : %g\n",user.cells[0],(double)Lx,(double)dx,(double)dt);CHKERRQ(ierr);
ierr = PetscPrintf(comm, " Tmax : %g avec dt : %g soit %f iterations\n",(double)tmax,(double)dt,(double)tmax/dt);CHKERRQ(ierr);
ierr = PetscPrintf(comm, "Parametres gaussienne : \n");CHKERRQ(ierr);
ierr = PetscPrintf(comm, " Sigma : %g x0 : %g y0 : %g\n",(double)user.sigma,(double)user.pos_init[0],(double)user.pos_init[1]);CHKERRQ(ierr);
ierr = PetscPrintf(comm, "Parametres physiques : \n");CHKERRQ(ierr);
ierr = PetscPrintf(comm, " Um : %g Diffu : %g\n",(double)Um,(double)user.Diffu);CHKERRQ(ierr);
ierr = PetscPrintf(comm, " Peclet de maille : %g\n",(double)Peclet);CHKERRQ(ierr);
ierr = PetscPrintf(comm, " CFL : %g\n",(double)cfl);CHKERRQ(ierr);
/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
sur le maillage
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
ierr = CreateMesh(comm, &user, &base);CHKERRQ(ierr);
ierr = DMSetApplicationContext(base, &user);CHKERRQ(ierr);
ierr = PetscMalloc2(1, &user.exactConc,1,&user.bord_d);CHKERRQ(ierr);
ierr = SetupDiscretization(base, &user);CHKERRQ(ierr);
ierr = AddIdentityLabel(base);CHKERRQ(ierr);
ierr = DMViewFromOptions(base,NULL,"-dm_pre_adapt");CHKERRQ(ierr);
ierr = DMCreateGlobalVector(base, &u);CHKERRQ(ierr);
ierr = PetscObjectSetName((PetscObject) u, "conc");CHKERRQ(ierr);
ierr = DMProjectFunction(base, 0.0, user.exactConc, ctx, INSERT_ALL_VALUES, u);CHKERRQ(ierr);
ierr = VecViewFromOptions(u,NULL,"-vec_pre_adapt");CHKERRQ(ierr);
/* adapt */
ierr = CreateAdaptivityLabel(base, &adaptLabel, u, 1.);CHKERRQ(ierr);
ierr = DMForestTemplate(base,comm,&postForest);CHKERRQ(ierr);
ierr = DMForestSetMinimumRefinement(postForest,0);CHKERRQ(ierr);
ierr = DMForestSetInitialRefinement(postForest,0);CHKERRQ(ierr);
ierr = DMForestSetAdaptivityLabel(postForest,adaptLabel);CHKERRQ(ierr);
ierr = DMLabelDestroy(&adaptLabel);CHKERRQ(ierr);
ierr = DMSetUp(postForest);CHKERRQ(ierr);
ierr = DMViewFromOptions(postForest,NULL,"-dm_post_view");CHKERRQ(ierr);
/* transfer */
ierr = DMCreateGlobalVector(postForest,&uNew);CHKERRQ(ierr);
ierr = DMForestTransferVec(base,u,postForest,uNew,PETSC_TRUE,0.0);CHKERRQ(ierr);
ierr = VecViewFromOptions(uNew,NULL,"-vec_post_transfer_view");CHKERRQ(ierr);
GetDiscretizationInfo(postForest);
ierr = DMGetDS(base, &prob);CHKERRQ(ierr);
ierr = DMSetDS(postForest, prob);CHKERRQ(ierr);
initializeTS(postForest,&user,&ts,tmax,dt);
ierr = TSSetFromOptions(ts);CHKERRQ(ierr);
ierr = TSSetSolution(ts,uNew);CHKERRQ(ierr);
ierr = VecDuplicate(uNew, &r);CHKERRQ(ierr);
/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
Solve linear system
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
TSSolve(ts,uNew);
ierr = TSGetTime(ts, &ftime);CHKERRQ(ierr);
ierr = TSGetStepNumber(ts, &steps);CHKERRQ(ierr);
ierr = DMComputeL2Diff(postForest, ftime, user.exactConc, ctx, uNew, &error);CHKERRQ(ierr);
if (error < 1.0e-12) {ierr = PetscPrintf(PETSC_COMM_WORLD, "L_2 Error: < 1.0e-12\n");CHKERRQ(ierr);}
else {ierr = PetscPrintf(PETSC_COMM_WORLD, "L_2 Error: %g at time %f\n",(double) error,(double) ftime);CHKERRQ(ierr);}
ierr = VecViewFromOptions(uNew, NULL, "-vec_sol_u");CHKERRQ(ierr);
ierr = VecDuplicate(uNew, &u_exact);CHKERRQ(ierr);
ierr = PetscObjectSetName((PetscObject) u_exact, "conc_ex");CHKERRQ(ierr);
ierr = DMProjectFunction(postForest, ftime, user.exactConc, ctx, INSERT_ALL_VALUES, u_exact);CHKERRQ(ierr);
ierr = VecViewFromOptions(u_exact, NULL, "-vec_sol_exact");CHKERRQ(ierr);
ierr = VecDestroy(&u);CHKERRQ(ierr);
ierr = VecDestroy(&uNew);CHKERRQ(ierr);
ierr = VecDestroy(&u_exact);CHKERRQ(ierr);
ierr = VecDestroy(&r);CHKERRQ(ierr);
ierr = TSDestroy(&ts);CHKERRQ(ierr);
ierr = DMDestroy(&base);CHKERRQ(ierr);
ierr = DMDestroy(&postForest);CHKERRQ(ierr);
ierr = PetscFree2(user.exactConc,user.bord_d);CHKERRQ(ierr);
ierr = PetscFinalize();
return ierr;
}
More information about the petsc-users
mailing list