Actual source code: ex12.c

petsc-3.6.1 2015-07-22
Report Typos and Errors
  1: static char help[] = "Poisson Problem in 2d and 3d with simplicial finite elements.\n\
  2: We solve the Poisson problem in a rectangular\n\
  3: domain, using a parallel unstructured mesh (DMPLEX) to discretize it.\n\
  4: This example supports discretized auxiliary fields (conductivity) as well as\n\
  5: multilevel nonlinear solvers.\n\n\n";

  7: #include <petscdmplex.h>
  8: #include <petscsnes.h>
  9: #include <petscds.h>
 10: #include <petscviewerhdf5.h>

 12: typedef enum {NEUMANN, DIRICHLET, NONE} BCType;
 13: typedef enum {RUN_FULL, RUN_TEST, RUN_PERF} RunType;
 14: typedef enum {COEFF_NONE, COEFF_ANALYTIC, COEFF_FIELD, COEFF_NONLINEAR} CoeffType;

 16: typedef struct {
 17:   PetscInt      debug;             /* The debugging level */
 18:   RunType       runType;           /* Whether to run tests, or solve the full problem */
 19:   PetscBool     jacobianMF;        /* Whether to calculate the Jacobian action on the fly */
 20:   PetscLogEvent createMeshEvent;
 21:   PetscBool     showInitial, showSolution, restart, check;
 22:   PetscViewer   checkpoint;
 23:   /* Domain and mesh definition */
 24:   PetscInt      dim;               /* The topological mesh dimension */
 25:   char          filename[2048];    /* The optional ExodusII file */
 26:   PetscBool     interpolate;       /* Generate intermediate mesh elements */
 27:   PetscReal     refinementLimit;   /* The largest allowable cell volume */
 28:   /* Problem definition */
 29:   BCType        bcType;
 30:   CoeffType     variableCoefficient;
 31:   PetscErrorCode (**exactFuncs)(PetscInt dim, const PetscReal x[], PetscInt Nf, PetscScalar *u, void *ctx);
 32: } AppCtx;

 34: PetscErrorCode zero(PetscInt dim, const PetscReal x[], PetscInt Nf, PetscScalar *u, void *ctx)
 35: {
 36:   u[0] = 0.0;
 37:   return 0;
 38: }

 40: /*
 41:   In 2D for Dirichlet conditions, we use exact solution:

 43:     u = x^2 + y^2
 44:     f = 4

 46:   so that

 48:     -\Delta u + f = -4 + 4 = 0

 50:   For Neumann conditions, we have

 52:     -\nabla u \cdot -\hat y |_{y=0} =  (2y)|_{y=0} =  0 (bottom)
 53:     -\nabla u \cdot  \hat y |_{y=1} = -(2y)|_{y=1} = -2 (top)
 54:     -\nabla u \cdot -\hat x |_{x=0} =  (2x)|_{x=0} =  0 (left)
 55:     -\nabla u \cdot  \hat x |_{x=1} = -(2x)|_{x=1} = -2 (right)

 57:   Which we can express as

 59:     \nabla u \cdot  \hat n|_\Gamma = {2 x, 2 y} \cdot \hat n = 2 (x + y)
 60: */
 61: PetscErrorCode quadratic_u_2d(PetscInt dim, const PetscReal x[], PetscInt Nf, PetscScalar *u, void *ctx)
 62: {
 63:   *u = x[0]*x[0] + x[1]*x[1];
 64:   return 0;
 65: }

 67: void f0_u(PetscInt dim, PetscInt Nf, PetscInt NfAux,
 68:           const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
 69:           const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
 70:           PetscReal t, const PetscReal x[], PetscScalar f0[])
 71: {
 72:   f0[0] = 4.0;
 73: }

 75: void f0_bd_u(PetscInt dim, PetscInt Nf, PetscInt NfAux,
 76:              const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
 77:              const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
 78:              PetscReal t, const PetscReal x[], const PetscReal n[], PetscScalar f0[])
 79: {
 80:   PetscInt d;
 81:   for (d = 0, f0[0] = 0.0; d < dim; ++d) f0[0] += -n[d]*2.0*x[d];
 82: }

 84: void f0_bd_zero(PetscInt dim, PetscInt Nf, PetscInt NfAux,
 85:                 const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
 86:                 const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
 87:                 PetscReal t, const PetscReal x[], const PetscReal n[], PetscScalar f0[])
 88: {
 89:   f0[0] = 0.0;
 90: }

 92: void f1_bd_zero(PetscInt dim, PetscInt Nf, PetscInt NfAux,
 93:                 const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
 94:                 const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
 95:                 PetscReal t, const PetscReal x[], const PetscReal n[], PetscScalar f1[])
 96: {
 97:   PetscInt comp;
 98:   for (comp = 0; comp < dim; ++comp) f1[comp] = 0.0;
 99: }

101: /* gradU[comp*dim+d] = {u_x, u_y} or {u_x, u_y, u_z} */
102: void f1_u(PetscInt dim, PetscInt Nf, PetscInt NfAux,
103:           const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
104:           const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
105:           PetscReal t, const PetscReal x[], PetscScalar f1[])
106: {
107:   PetscInt d;
108:   for (d = 0; d < dim; ++d) f1[d] = u_x[d];
109: }

111: /* < \nabla v, \nabla u + {\nabla u}^T >
112:    This just gives \nabla u, give the perdiagonal for the transpose */
113: void g3_uu(PetscInt dim, PetscInt Nf, PetscInt NfAux,
114:            const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
115:            const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
116:            PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscScalar g3[])
117: {
118:   PetscInt d;
119:   for (d = 0; d < dim; ++d) g3[d*dim+d] = 1.0;
120: }

122: /*
123:   In 2D for Dirichlet conditions with a variable coefficient, we use exact solution:

125:     u  = x^2 + y^2
126:     f  = 6 (x + y)
127:     nu = (x + y)

129:   so that

131:     -\div \nu \grad u + f = -6 (x + y) + 6 (x + y) = 0
132: */
133: PetscErrorCode nu_2d(PetscInt dim, const PetscReal x[], PetscInt Nf, PetscScalar *u, void *ctx)
134: {
135:   *u = x[0] + x[1];
136:   return 0;
137: }

139: void f0_analytic_u(PetscInt dim, PetscInt Nf, PetscInt NfAux,
140:                    const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
141:                    const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
142:                    PetscReal t, const PetscReal x[], PetscScalar f0[])
143: {
144:   f0[0] = 6.0*(x[0] + x[1]);
145: }

147: /* gradU[comp*dim+d] = {u_x, u_y} or {u_x, u_y, u_z} */
148: void f1_analytic_u(PetscInt dim, PetscInt Nf, PetscInt NfAux,
149:                    const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
150:                    const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
151:                    PetscReal t, const PetscReal x[], PetscScalar f1[])
152: {
153:   PetscInt d;
154:   for (d = 0; d < dim; ++d) f1[d] = (x[0] + x[1])*u_x[d];
155: }

157: void f1_field_u(PetscInt dim, PetscInt Nf, PetscInt NfAux,
158:                 const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
159:                 const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
160:                 PetscReal t, const PetscReal x[], PetscScalar f1[])
161: {
162:   PetscInt d;
163:   for (d = 0; d < dim; ++d) f1[d] = a[0]*u_x[d];
164: }

166: /* < \nabla v, \nabla u + {\nabla u}^T >
167:    This just gives \nabla u, give the perdiagonal for the transpose */
168: void g3_analytic_uu(PetscInt dim, PetscInt Nf, PetscInt NfAux,
169:                     const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
170:                     const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
171:                     PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscScalar g3[])
172: {
173:   PetscInt d;
174:   for (d = 0; d < dim; ++d) g3[d*dim+d] = x[0] + x[1];
175: }

177: void g3_field_uu(PetscInt dim, PetscInt Nf, PetscInt NfAux,
178:                  const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
179:                  const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
180:                  PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscScalar g3[])
181: {
182:   PetscInt d;
183:   for (d = 0; d < dim; ++d) g3[d*dim+d] = a[0];
184: }

186: /*
187:   In 2D for Dirichlet conditions with a nonlinear coefficient (p-Laplacian with p = 4), we use exact solution:

189:     u  = x^2 + y^2
190:     f  = 16 (x^2 + y^2)
191:     nu = 1/2 |grad u|^2

193:   so that

195:     -\div \nu \grad u + f = -16 (x^2 + y^2) + 16 (x^2 + y^2) = 0
196: */
197: void f0_analytic_nonlinear_u(PetscInt dim, PetscInt Nf, PetscInt NfAux,
198:                              const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
199:                              const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
200:                              PetscReal t, const PetscReal x[], PetscScalar f0[])
201: {
202:   f0[0] = 16.0*(x[0]*x[0] + x[1]*x[1]);
203: }

205: /* gradU[comp*dim+d] = {u_x, u_y} or {u_x, u_y, u_z} */
206: void f1_analytic_nonlinear_u(PetscInt dim, PetscInt Nf, PetscInt NfAux,
207:                              const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
208:                              const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
209:                              PetscReal t, const PetscReal x[], PetscScalar f1[])
210: {
211:   PetscScalar nu = 0.0;
212:   PetscInt    d;
213:   for (d = 0; d < dim; ++d) nu += u_x[d]*u_x[d];
214:   for (d = 0; d < dim; ++d) f1[d] = 0.5*nu*u_x[d];
215: }

217: /*
218:   grad (u + eps w) - grad u = eps grad w

220:   1/2 |grad (u + eps w)|^2 grad (u + eps w) - 1/2 |grad u|^2 grad u
221: = 1/2 (|grad u|^2 + 2 eps <grad u,grad w>) (grad u + eps grad w) - 1/2 |grad u|^2 grad u
222: = 1/2 (eps |grad u|^2 grad w + 2 eps <grad u,grad w> grad u)
223: = eps (1/2 |grad u|^2 grad w + grad u <grad u,grad w>)
224: */
225: void g3_analytic_nonlinear_uu(PetscInt dim, PetscInt Nf, PetscInt NfAux,
226:                               const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
227:                               const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
228:                               PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscScalar g3[])
229: {
230:   PetscScalar nu = 0.0;
231:   PetscInt    d, e;
232:   for (d = 0; d < dim; ++d) nu += u_x[d]*u_x[d];
233:   for (d = 0; d < dim; ++d) {
234:     g3[d*dim+d] = 0.5*nu;
235:     for (e = 0; e < dim; ++e) {
236:       g3[d*dim+e] += u_x[d]*u_x[e];
237:     }
238:   }
239: }

241: /*
242:   In 3D for Dirichlet conditions we use exact solution:

244:     u = x^2 + y^2 + z^2
245:     f = 6

247:   so that

249:     -\Delta u + f = -6 + 6 = 0

251:   For Neumann conditions, we have

253:     -\nabla u \cdot -\hat z |_{z=0} =  (2z)|_{z=0} =  0 (bottom)
254:     -\nabla u \cdot  \hat z |_{z=1} = -(2z)|_{z=1} = -2 (top)
255:     -\nabla u \cdot -\hat y |_{y=0} =  (2y)|_{y=0} =  0 (front)
256:     -\nabla u \cdot  \hat y |_{y=1} = -(2y)|_{y=1} = -2 (back)
257:     -\nabla u \cdot -\hat x |_{x=0} =  (2x)|_{x=0} =  0 (left)
258:     -\nabla u \cdot  \hat x |_{x=1} = -(2x)|_{x=1} = -2 (right)

260:   Which we can express as

262:     \nabla u \cdot  \hat n|_\Gamma = {2 x, 2 y, 2z} \cdot \hat n = 2 (x + y + z)
263: */
264: PetscErrorCode quadratic_u_3d(PetscInt dim, const PetscReal x[], PetscInt Nf, PetscScalar *u, void *ctx)
265: {
266:   *u = x[0]*x[0] + x[1]*x[1] + x[2]*x[2];
267:   return 0;
268: }

272: PetscErrorCode ProcessOptions(MPI_Comm comm, AppCtx *options)
273: {
274:   const char    *bcTypes[3]  = {"neumann", "dirichlet", "none"};
275:   const char    *runTypes[3] = {"full", "test", "perf"};
276:   const char    *coeffTypes[4] = {"none", "analytic", "field", "nonlinear"};
277:   PetscInt       bc, run, coeff;
278:   PetscBool      flg;

282:   options->debug               = 0;
283:   options->runType             = RUN_FULL;
284:   options->dim                 = 2;
285:   options->filename[0]         = '\0';
286:   options->interpolate         = PETSC_FALSE;
287:   options->refinementLimit     = 0.0;
288:   options->bcType              = DIRICHLET;
289:   options->variableCoefficient = COEFF_NONE;
290:   options->jacobianMF          = PETSC_FALSE;
291:   options->showInitial         = PETSC_FALSE;
292:   options->showSolution        = PETSC_FALSE;
293:   options->restart             = PETSC_FALSE;
294:   options->check               = PETSC_FALSE;
295:   options->checkpoint          = NULL;

297:   PetscOptionsBegin(comm, "", "Poisson Problem Options", "DMPLEX");
298:   PetscOptionsInt("-debug", "The debugging level", "ex12.c", options->debug, &options->debug, NULL);
299:   run  = options->runType;
300:   PetscOptionsEList("-run_type", "The run type", "ex12.c", runTypes, 3, runTypes[options->runType], &run, NULL);

302:   options->runType = (RunType) run;

304:   PetscOptionsInt("-dim", "The topological mesh dimension", "ex12.c", options->dim, &options->dim, NULL);
305:   PetscOptionsString("-f", "Exodus.II filename to read", "ex12.c", options->filename, options->filename, sizeof(options->filename), &flg);
306: #if !defined(PETSC_HAVE_EXODUSII)
307:   if (flg) SETERRQ(comm, PETSC_ERR_ARG_WRONG, "This option requires ExodusII support. Reconfigure using --download-exodusii");
308: #endif
309:   PetscOptionsBool("-interpolate", "Generate intermediate mesh elements", "ex12.c", options->interpolate, &options->interpolate, NULL);
310:   PetscOptionsReal("-refinement_limit", "The largest allowable cell volume", "ex12.c", options->refinementLimit, &options->refinementLimit, NULL);
311:   bc   = options->bcType;
312:   PetscOptionsEList("-bc_type","Type of boundary condition","ex12.c",bcTypes,3,bcTypes[options->bcType],&bc,NULL);
313:   options->bcType = (BCType) bc;
314:   coeff = options->variableCoefficient;
315:   PetscOptionsEList("-variable_coefficient","Type of variable coefficent","ex12.c",coeffTypes,4,coeffTypes[options->variableCoefficient],&coeff,NULL);
316:   options->variableCoefficient = (CoeffType) coeff;

318:   PetscOptionsBool("-jacobian_mf", "Calculate the action of the Jacobian on the fly", "ex12.c", options->jacobianMF, &options->jacobianMF, NULL);
319:   PetscOptionsBool("-show_initial", "Output the initial guess for verification", "ex12.c", options->showInitial, &options->showInitial, NULL);
320:   PetscOptionsBool("-show_solution", "Output the solution for verification", "ex12.c", options->showSolution, &options->showSolution, NULL);
321:   PetscOptionsBool("-restart", "Read in the mesh and solution from a file", "ex12.c", options->restart, &options->restart, NULL);
322:   PetscOptionsBool("-check", "Compare with default integration routines", "ex12.c", options->check, &options->check, NULL);
323:   PetscOptionsEnd();

325:   PetscLogEventRegister("CreateMesh", DM_CLASSID, &options->createMeshEvent);

327:   if (options->restart) {
328:     PetscViewerCreate(comm, &options->checkpoint);
329:     PetscViewerSetType(options->checkpoint, PETSCVIEWERHDF5);
330:     PetscViewerFileSetMode(options->checkpoint, FILE_MODE_READ);
331:     PetscViewerFileSetName(options->checkpoint, options->filename);
332:   }
333:   return(0);
334: }

338: PetscErrorCode CreateMesh(MPI_Comm comm, AppCtx *user, DM *dm)
339: {
340:   PetscInt       dim             = user->dim;
341:   const char    *filename        = user->filename;
342:   PetscBool      interpolate     = user->interpolate;
343:   PetscReal      refinementLimit = user->refinementLimit;
344:   size_t         len;

348:   PetscLogEventBegin(user->createMeshEvent,0,0,0,0);
349:   PetscStrlen(filename, &len);
350:   if (!len) {
351:     DMPlexCreateBoxMesh(comm, dim, interpolate, dm);
352:     PetscObjectSetName((PetscObject) *dm, "Mesh");
353:   } else if (user->checkpoint) {
354:     DMCreate(comm, dm);
355:     DMSetType(*dm, DMPLEX);
356:     DMLoad(*dm, user->checkpoint);
357:     DMPlexSetRefinementUniform(*dm, PETSC_FALSE);
358:   } else {
359:     PetscMPIInt rank;

361:     MPI_Comm_rank(comm, &rank);
362:     DMPlexCreateFromFile(comm, filename, interpolate, dm);
363:     DMPlexSetRefinementUniform(*dm, PETSC_FALSE);
364:     /* Must have boundary marker for Dirichlet conditions */
365:   }
366:   {
367:     DM refinedMesh     = NULL;
368:     DM distributedMesh = NULL;

370:     /* Refine mesh using a volume constraint */
371:     DMPlexSetRefinementLimit(*dm, refinementLimit);
372:     DMRefine(*dm, comm, &refinedMesh);
373:     if (refinedMesh) {
374:       const char *name;

376:       PetscObjectGetName((PetscObject) *dm,         &name);
377:       PetscObjectSetName((PetscObject) refinedMesh,  name);
378:       DMDestroy(dm);
379:       *dm  = refinedMesh;
380:     }
381:     /* Distribute mesh over processes */
382:     DMPlexDistribute(*dm, 0, NULL, &distributedMesh);
383:     if (distributedMesh) {
384:       DMDestroy(dm);
385:       *dm  = distributedMesh;
386:     }
387:   }
388:   if (user->bcType == NEUMANN) {
389:     DMLabel label;

391:     DMPlexCreateLabel(*dm, "boundary");
392:     DMPlexGetLabel(*dm, "boundary", &label);
393:     DMPlexMarkBoundaryFaces(*dm, label);
394:   }
395:   DMSetFromOptions(*dm);
396:   DMViewFromOptions(*dm, NULL, "-dm_view");
397:   PetscLogEventEnd(user->createMeshEvent,0,0,0,0);
398:   return(0);
399: }

403: PetscErrorCode SetupProblem(DM dm, AppCtx *user)
404: {
405:   PetscDS        prob;

409:   DMGetDS(dm, &prob);
410:   switch (user->variableCoefficient) {
411:   case COEFF_NONE:
412:     PetscDSSetResidual(prob, 0, f0_u, f1_u);
413:     PetscDSSetJacobian(prob, 0, 0, NULL, NULL, NULL, g3_uu);
414:     break;
415:   case COEFF_ANALYTIC:
416:     PetscDSSetResidual(prob, 0, f0_analytic_u, f1_analytic_u);
417:     PetscDSSetJacobian(prob, 0, 0, NULL, NULL, NULL, g3_analytic_uu);
418:     break;
419:   case COEFF_FIELD:
420:     PetscDSSetResidual(prob, 0, f0_analytic_u, f1_field_u);
421:     PetscDSSetJacobian(prob, 0, 0, NULL, NULL, NULL, g3_field_uu);
422:     break;
423:   case COEFF_NONLINEAR:
424:     PetscDSSetResidual(prob, 0, f0_analytic_nonlinear_u, f1_analytic_nonlinear_u);
425:     PetscDSSetJacobian(prob, 0, 0, NULL, NULL, NULL, g3_analytic_nonlinear_uu);
426:     break;
427:   default: SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Invalid variable coefficient type %d", user->variableCoefficient);
428:   }
429:   switch (user->dim) {
430:   case 2:
431:     user->exactFuncs[0] = quadratic_u_2d;
432:     if (user->bcType == NEUMANN) {PetscDSSetBdResidual(prob, 0, f0_bd_u, f1_bd_zero);}
433:     break;
434:   case 3:
435:     user->exactFuncs[0] = quadratic_u_3d;
436:     if (user->bcType == NEUMANN) {PetscDSSetBdResidual(prob, 0, f0_bd_u, f1_bd_zero);}
437:     break;
438:   default:
439:     SETERRQ1(PETSC_COMM_WORLD, PETSC_ERR_ARG_OUTOFRANGE, "Invalid dimension %d", user->dim);
440:   }
441:   return(0);
442: }

446: PetscErrorCode SetupMaterial(DM dm, DM dmAux, AppCtx *user)
447: {
448:   PetscErrorCode (*matFuncs[1])(PetscInt dim, const PetscReal x[], PetscInt Nf, PetscScalar u[], void *ctx) = {nu_2d};
449:   Vec            nu;

453:   DMCreateLocalVector(dmAux, &nu);
454:   DMPlexProjectFunctionLocal(dmAux, matFuncs, NULL, INSERT_ALL_VALUES, nu);
455:   PetscObjectCompose((PetscObject) dm, "A", (PetscObject) nu);
456:   VecDestroy(&nu);
457:   return(0);
458: }

462: PetscErrorCode SetupDiscretization(DM dm, AppCtx *user)
463: {
464:   DM             cdm   = dm;
465:   const PetscInt dim   = user->dim;
466:   const PetscInt id    = 1;
467:   PetscFE        feAux = NULL;
468:   PetscFE        feBd  = NULL;
469:   PetscFE        feCh  = NULL;
470:   PetscFE        fe;
471:   PetscDS        prob;

475:   /* Create finite element */
476:   PetscFECreateDefault(dm, dim, 1, PETSC_TRUE, NULL, -1, &fe);
477:   PetscObjectSetName((PetscObject) fe, "potential");
478:   if (user->bcType == NEUMANN) {
479:     PetscFECreateDefault(dm, dim-1, 1, PETSC_TRUE, "bd_", -1, &feBd);
480:     PetscObjectSetName((PetscObject) feBd, "potential");
481:   }
482:   if (user->variableCoefficient == COEFF_FIELD) {
483:     PetscQuadrature q;

485:     PetscFECreateDefault(dm, dim, 1, PETSC_TRUE, "mat_", -1, &feAux);
486:     PetscFEGetQuadrature(fe, &q);
487:     PetscFESetQuadrature(feAux, q);
488:   }
489:   if (user->check) {PetscFECreateDefault(dm, dim, 1, PETSC_TRUE, "ch_", -1, &feCh);}
490:   /* Set discretization and boundary conditions for each mesh */
491:   while (cdm) {
492:     DMGetDS(cdm, &prob);
493:     PetscDSSetDiscretization(prob, 0, (PetscObject) fe);
494:     PetscDSSetBdDiscretization(prob, 0, (PetscObject) feBd);
495:     if (feAux) {
496:       DM      dmAux;
497:       PetscDS probAux;

499:       DMClone(cdm, &dmAux);
500:       DMPlexCopyCoordinates(cdm, dmAux);
501:       DMGetDS(dmAux, &probAux);
502:       PetscDSSetDiscretization(probAux, 0, (PetscObject) feAux);
503:       PetscObjectCompose((PetscObject) dm, "dmAux", (PetscObject) dmAux);
504:       SetupMaterial(cdm, dmAux, user);
505:       DMDestroy(&dmAux);
506:     }
507:     if (feCh) {
508:       DM      dmCh;
509:       PetscDS probCh;

511:       DMClone(cdm, &dmCh);
512:       DMPlexCopyCoordinates(cdm, dmCh);
513:       DMGetDS(dmCh, &probCh);
514:       PetscDSSetDiscretization(probCh, 0, (PetscObject) feCh);
515:       PetscObjectCompose((PetscObject) dm, "dmCh", (PetscObject) dmCh);
516:       DMDestroy(&dmCh);
517:     }
518:     SetupProblem(cdm, user);
519:     DMPlexAddBoundary(cdm, user->bcType == DIRICHLET ? PETSC_TRUE : PETSC_FALSE, "wall", user->bcType == NEUMANN ? "boundary" : "marker", 0, 0, NULL, (void (*)()) user->exactFuncs[0], 1, &id, user);
520:     DMPlexGetCoarseDM(cdm, &cdm);
521:   }
522:   PetscFEDestroy(&fe);
523:   PetscFEDestroy(&feBd);
524:   PetscFEDestroy(&feAux);
525:   PetscFEDestroy(&feCh);
526:   return(0);
527: }

531: int main(int argc, char **argv)
532: {
533:   DM             dm;          /* Problem specification */
534:   SNES           snes;        /* nonlinear solver */
535:   Vec            u,r;         /* solution, residual vectors */
536:   Mat            A,J;         /* Jacobian matrix */
537:   MatNullSpace   nullSpace;   /* May be necessary for Neumann conditions */
538:   AppCtx         user;        /* user-defined work context */
539:   JacActionCtx   userJ;       /* context for Jacobian MF action */
540:   PetscInt       its;         /* iterations for convergence */
541:   PetscReal      error = 0.0; /* L_2 error in the solution */

544:   PetscInitialize(&argc, &argv, NULL, help);
545:   ProcessOptions(PETSC_COMM_WORLD, &user);
546:   SNESCreate(PETSC_COMM_WORLD, &snes);
547:   CreateMesh(PETSC_COMM_WORLD, &user, &dm);
548:   SNESSetDM(snes, dm);
549:   DMSetApplicationContext(dm, &user);

551:   PetscMalloc(1 * sizeof(void (*)(const PetscReal[], PetscScalar *, void *)), &user.exactFuncs);
552:   SetupDiscretization(dm, &user);

554:   DMCreateGlobalVector(dm, &u);
555:   PetscObjectSetName((PetscObject) u, "potential");
556:   VecDuplicate(u, &r);

558:   DMSetMatType(dm,MATAIJ);
559:   DMCreateMatrix(dm, &J);
560:   if (user.jacobianMF) {
561:     PetscInt M, m, N, n;

563:     MatGetSize(J, &M, &N);
564:     MatGetLocalSize(J, &m, &n);
565:     MatCreate(PETSC_COMM_WORLD, &A);
566:     MatSetSizes(A, m, n, M, N);
567:     MatSetType(A, MATSHELL);
568:     MatSetUp(A);
569: #if 0
570:     MatShellSetOperation(A, MATOP_MULT, (void (*)(void))FormJacobianAction);
571: #endif

573:     userJ.dm   = dm;
574:     userJ.J    = J;
575:     userJ.user = &user;

577:     DMCreateLocalVector(dm, &userJ.u);
578:     DMPlexProjectFunctionLocal(dm, user.exactFuncs, NULL, INSERT_BC_VALUES, userJ.u);
579:     MatShellSetContext(A, &userJ);
580:   } else {
581:     A = J;
582:   }
583:   if (user.bcType == NEUMANN) {
584:     MatNullSpaceCreate(PetscObjectComm((PetscObject) dm), PETSC_TRUE, 0, NULL, &nullSpace);
585:     MatSetNullSpace(A, nullSpace);
586:   }

588:   DMSNESSetFunctionLocal(dm,  (PetscErrorCode (*)(DM,Vec,Vec,void*)) DMPlexSNESComputeResidualFEM, &user);
589:   DMSNESSetJacobianLocal(dm,  (PetscErrorCode (*)(DM,Vec,Mat,Mat,void*)) DMPlexSNESComputeJacobianFEM, &user);
590:   SNESSetJacobian(snes, A, J, NULL, NULL);

592:   SNESSetFromOptions(snes);

594:   DMPlexProjectFunction(dm, user.exactFuncs, NULL, INSERT_ALL_VALUES, u);
595:   if (user.checkpoint) {
596: #if defined(PETSC_HAVE_HDF5)
597:     PetscViewerHDF5PushGroup(user.checkpoint, "/fields");
598:     VecLoad(u, user.checkpoint);
599:     PetscViewerHDF5PopGroup(user.checkpoint);
600: #endif
601:   }
602:   PetscViewerDestroy(&user.checkpoint);
603:   if (user.showInitial) {
604:     Vec lv;
605:     DMGetLocalVector(dm, &lv);
606:     DMGlobalToLocalBegin(dm, u, INSERT_VALUES, lv);
607:     DMGlobalToLocalEnd(dm, u, INSERT_VALUES, lv);
608:     DMPrintLocalVec(dm, "Local function", 1.0e-10, lv);
609:     DMRestoreLocalVector(dm, &lv);
610:   }
611:   if (user.runType == RUN_FULL) {
612:     PetscErrorCode (*initialGuess[1])(PetscInt dim, const PetscReal x[], PetscInt Nf, PetscScalar u[], void *ctx) = {zero};

614:     DMPlexProjectFunction(dm, initialGuess, NULL, INSERT_VALUES, u);
615:     if (user.debug) {
616:       PetscPrintf(PETSC_COMM_WORLD, "Initial guess\n");
617:       VecView(u, PETSC_VIEWER_STDOUT_WORLD);
618:     }
619:     SNESSolve(snes, NULL, u);
620:     SNESGetIterationNumber(snes, &its);
621:     PetscPrintf(PETSC_COMM_WORLD, "Number of SNES iterations = %D\n", its);
622:     DMPlexComputeL2Diff(dm, user.exactFuncs, NULL, u, &error);
623:     if (error < 1.0e-11) {PetscPrintf(PETSC_COMM_WORLD, "L_2 Error: < 1.0e-11\n");}
624:     else                 {PetscPrintf(PETSC_COMM_WORLD, "L_2 Error: %g\n", error);}
625:     if (user.showSolution) {
626:       PetscPrintf(PETSC_COMM_WORLD, "Solution\n");
627:       VecChop(u, 3.0e-9);
628:       VecView(u, PETSC_VIEWER_STDOUT_WORLD);
629:     }
630:   } else if (user.runType == RUN_PERF) {
631:     PetscReal res = 0.0;

633:     SNESComputeFunction(snes, u, r);
634:     PetscPrintf(PETSC_COMM_WORLD, "Initial Residual\n");
635:     VecChop(r, 1.0e-10);
636:     VecNorm(r, NORM_2, &res);
637:     PetscPrintf(PETSC_COMM_WORLD, "L_2 Residual: %g\n", res);
638:   } else {
639:     PetscReal res = 0.0;

641:     /* Check discretization error */
642:     PetscPrintf(PETSC_COMM_WORLD, "Initial guess\n");
643:     VecView(u, PETSC_VIEWER_STDOUT_WORLD);
644:     DMPlexComputeL2Diff(dm, user.exactFuncs, NULL, u, &error);
645:     if (error < 1.0e-11) {PetscPrintf(PETSC_COMM_WORLD, "L_2 Error: < 1.0e-11\n");}
646:     else                 {PetscPrintf(PETSC_COMM_WORLD, "L_2 Error: %g\n", error);}
647:     /* Check residual */
648:     SNESComputeFunction(snes, u, r);
649:     PetscPrintf(PETSC_COMM_WORLD, "Initial Residual\n");
650:     VecChop(r, 1.0e-10);
651:     VecView(r, PETSC_VIEWER_STDOUT_WORLD);
652:     VecNorm(r, NORM_2, &res);
653:     PetscPrintf(PETSC_COMM_WORLD, "L_2 Residual: %g\n", res);
654:     /* Check Jacobian */
655:     {
656:       Vec          b;

658:       SNESComputeJacobian(snes, u, A, A);
659:       VecDuplicate(u, &b);
660:       VecSet(r, 0.0);
661:       SNESComputeFunction(snes, r, b);
662:       MatMult(A, u, r);
663:       VecAXPY(r, 1.0, b);
664:       VecDestroy(&b);
665:       PetscPrintf(PETSC_COMM_WORLD, "Au - b = Au + F(0)\n");
666:       VecChop(r, 1.0e-10);
667:       VecView(r, PETSC_VIEWER_STDOUT_WORLD);
668:       VecNorm(r, NORM_2, &res);
669:       PetscPrintf(PETSC_COMM_WORLD, "Linear L_2 Residual: %g\n", res);
670:     }
671:   }
672:   VecViewFromOptions(u, NULL, "-vec_view");

674:   if (user.bcType == NEUMANN) {MatNullSpaceDestroy(&nullSpace);}
675:   if (user.jacobianMF) {VecDestroy(&userJ.u);}
676:   if (A != J) {MatDestroy(&A);}
677:   MatDestroy(&J);
678:   VecDestroy(&u);
679:   VecDestroy(&r);
680:   SNESDestroy(&snes);
681:   DMDestroy(&dm);
682:   PetscFree(user.exactFuncs);
683:   PetscFinalize();
684:   return 0;
685: }