Actual source code: ex52.c

  1: static char help[] = "Testbed for FEM operations on the GPU\n\n";

  3: #include<petscdmmesh.h>
  4: #include<petscsnes.h>

  6: #ifndef PETSC_HAVE_SIEVE
  7: #error This example requires Sieve. Reconfigure using --with-sieve
  8: #endif

 10: typedef struct {
 11:   DM            dm;                /* REQUIRED in order to use SNES evaluation functions */
 12:   PetscInt      debug;             /* The debugging level */
 13:   PetscMPIInt   rank;              /* The process rank */
 14:   PetscMPIInt   numProcs;          /* The number of processes */
 15:   PetscInt      dim;               /* The topological mesh dimension */
 16:   PetscBool     interpolate;       /* Generate intermediate mesh elements */
 17:   PetscReal     refinementLimit;   /* The largest allowable cell volume */
 18:   char          partitioner[2048]; /* The graph partitioner */
 19:   PetscBool     computeFunction;   /* The flag for computing a residual */
 20:   PetscBool     computeJacobian;   /* The flag for computing a Jacobian */
 21:   PetscBool     batch;             /* The flag for batch assembly */
 22:   /* Element quadrature */
 23:   PetscQuadrature q;
 24: } AppCtx;

 26: /*------------------------------------------------------------------------------
 27:   This code can be generated using 'bin/pythonscripts/PetscGenerateFEMQuadrature.py 2 1 src/snes/examples/tutorials/ex52.h'
 28:  -----------------------------------------------------------------------------*/
 29: #include "ex52.h"

 33: PetscErrorCode ProcessOptions(MPI_Comm comm, AppCtx *options) {

 37:   options->debug           = 0;
 38:   options->dim             = 2;
 39:   options->interpolate     = PETSC_FALSE;
 40:   options->refinementLimit = 0.0;
 41:   options->computeFunction = PETSC_FALSE;
 42:   options->computeJacobian = PETSC_FALSE;
 43:   options->batch           = PETSC_FALSE;

 45:   MPI_Comm_size(comm, &options->numProcs);
 46:   MPI_Comm_rank(comm, &options->rank);
 47:   PetscOptionsBegin(comm, "", "Bratu Problem Options", "DMMESH");
 48:   PetscOptionsInt("-debug", "The debugging level", "ex52.c", options->debug, &options->debug, PETSC_NULL);
 49:   PetscOptionsInt("-dim", "The topological mesh dimension", "ex52.c", options->dim, &options->dim, PETSC_NULL);
 50:   PetscOptionsBool("-interpolate", "Generate intermediate mesh elements", "ex52.c", options->interpolate, &options->interpolate, PETSC_NULL);
 51:   PetscOptionsReal("-refinement_limit", "The largest allowable cell volume", "ex52.c", options->refinementLimit, &options->refinementLimit, PETSC_NULL);
 52:   PetscStrcpy(options->partitioner, "chaco");
 53:   PetscOptionsString("-partitioner", "The graph partitioner", "pflotran.cxx", options->partitioner, options->partitioner, 2048, PETSC_NULL);
 54:   PetscOptionsBool("-compute_function", "Compute the residual", "ex52.c", options->computeFunction, &options->computeFunction, PETSC_NULL);
 55:   PetscOptionsBool("-compute_jacobian", "Compute the Jacobian", "ex52.c", options->computeJacobian, &options->computeJacobian, PETSC_NULL);
 56:   PetscOptionsBool("-batch", "Use the batch assembly method", "ex52.c", options->batch, &options->batch, PETSC_NULL);
 57:   PetscOptionsEnd();
 58:   return(0);
 59: };

 63: PetscErrorCode CreateMesh(MPI_Comm comm, AppCtx *user, DM *dm)
 64: {
 65:   PetscInt       dim             = user->dim;
 66:   PetscBool      interpolate     = user->interpolate;
 67:   PetscReal      refinementLimit = user->refinementLimit;
 68:   const char    *partitioner     = user->partitioner;

 72:   DMMeshCreateBoxMesh(comm, dim, interpolate, dm);
 73:   {
 74:     DM refinedMesh     = PETSC_NULL;
 75:     DM distributedMesh = PETSC_NULL;

 77:     /* Refine mesh using a volume constraint */
 78:     DMMeshRefine(*dm, refinementLimit, interpolate, &refinedMesh);
 79:     if (refinedMesh) {
 80:       DMDestroy(dm);
 81:       *dm  = refinedMesh;
 82:     }
 83:     /* Distribute mesh over processes */
 84:     DMMeshDistribute(*dm, partitioner, &distributedMesh);
 85:     if (distributedMesh) {
 86:       DMDestroy(dm);
 87:       *dm  = distributedMesh;
 88:     }
 89:   }
 90:   DMSetFromOptions(*dm);
 91:   user->dm = *dm;
 92:   return(0);
 93: }

 97: PetscErrorCode SetupQuadrature(AppCtx *user) {
 99:   switch(user->dim) {
100:   case 2:
101:     user->q.numQuadPoints = NUM_QUADRATURE_POINTS_0;
102:     user->q.quadPoints    = points_0;
103:     user->q.quadWeights   = weights_0;
104:     user->q.numBasisFuncs = NUM_BASIS_FUNCTIONS_0;
105:     user->q.basis         = Basis_0;
106:     user->q.basisDer      = BasisDerivatives_0;
107:     break;
108:   default:
109:     SETERRQ1(PETSC_COMM_WORLD, PETSC_ERR_ARG_OUTOFRANGE, "Invalid dimension %d", user->dim);
110:   }
111:   return(0);
112: }

116: PetscErrorCode SetupSection(DM dm, AppCtx *user) {
117:   PetscSection   section;
118:   /* These can be generated using config/PETSc/FEM.py */
119:   PetscInt       numDof_0[2] = {1, 0};
120:   PetscInt       numDof_1[3] = {1, 0, 0};
121:   PetscInt       numDof_2[4] = {1, 0, 0, 0};
122:   PetscInt       dim         = user->dim;
123:   PetscInt      *numDof;
124:   const char    *bcLabel = PETSC_NULL;

128:   switch(user->dim) {
129:   case 1:
130:     numDof = numDof_0;
131:     break;
132:   case 2:
133:     numDof = numDof_1;
134:     break;
135:   case 3:
136:     numDof = numDof_2;
137:     break;
138:   default:
139:     SETERRQ1(PETSC_COMM_WORLD, PETSC_ERR_ARG_OUTOFRANGE, "Invalid spatial dimension %d", user->dim);
140:   }
141:   //if (user->bcType == DIRICHLET) {
142:   //  bcLabel = "marker";
143:   //}
144:   DMMeshCreateSection(dm, dim, numDof, bcLabel, 1, &section);
145:   DMMeshSetSection(dm, "default", section);
146:   return(0);
147: }

151: /*
152:  dm - The mesh
153:  X  - Local intput vector
154:  F  - Local output vector
155:  */
156: PetscErrorCode FormFunctionLocal(DM dm, Vec X, Vec F, AppCtx *user)
157: {
158:   //PetscScalar    (*rhsFunc)(const PetscReal []) = user->rhsFunc;
159:   const PetscInt   debug         = user->debug;
160:   const PetscInt   dim           = user->dim;
161:   const PetscInt   numQuadPoints = user->q.numQuadPoints;
162:   const PetscReal *quadPoints    = user->q.quadPoints;
163:   const PetscReal *quadWeights   = user->q.quadWeights;
164:   const PetscInt   numBasisFuncs = user->q.numBasisFuncs;
165:   const PetscReal *basis         = user->q.basis;
166:   const PetscReal *basisDer      = user->q.basisDer;
167:   PetscReal       *coords, *v0, *J, *invJ, detJ;
168:   PetscScalar     *realSpaceDer, *fieldGrad, *elemVec;
169:   PetscInt         cStart, cEnd;
170:   PetscErrorCode   ierr;

173:   VecSet(F, 0.0);
174:   PetscMalloc3(dim,PetscScalar,&realSpaceDer,dim,PetscScalar,&fieldGrad,numBasisFuncs,PetscScalar,&elemVec);
175:   PetscMalloc4(dim,PetscReal,&coords,dim,PetscReal,&v0,dim*dim,PetscReal,&J,dim*dim,PetscReal,&invJ);
176:   DMMeshGetHeightStratum(dm, 0, &cStart, &cEnd);
177:   for(PetscInt c = cStart; c < cEnd; ++c) {
178:     const PetscScalar *x;

180:     PetscMemzero(elemVec, numBasisFuncs * sizeof(PetscScalar));
181:     DMMeshComputeCellGeometry(dm, c, v0, J, invJ, &detJ);
182:     if (detJ <= 0.0) {SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid determinant %g for element %d", detJ, c);}
183:     DMMeshVecGetClosure(dm, X, c, &x);
184:     if (debug) {DMMeshPrintCellVector(c, "Solution", numBasisFuncs, x);}

186:     for(int q = 0; q < numQuadPoints; ++q) {
187:       PetscScalar fieldVal = 0.0;

189:       if (debug) {PetscPrintf(PETSC_COMM_SELF, "  quad point %d\n", q);}
190:       for(int d = 0; d < dim; ++d) {
191:         fieldGrad[d] = 0.0;
192:         coords[d] = v0[d];
193:         for(int e = 0; e < dim; ++e) {
194:           coords[d] += J[d*dim+e]*(quadPoints[q*dim+e] + 1.0);
195:         }
196:         if (debug) {PetscPrintf(PETSC_COMM_SELF, "    coords[%d] %g\n", d, coords[d]);}
197:       }
198:       for(int f = 0; f < numBasisFuncs; ++f) {
199:         fieldVal += x[f]*basis[q*numBasisFuncs+f];

201:         for(int d = 0; d < dim; ++d) {
202:           realSpaceDer[d] = 0.0;
203:           for(int e = 0; e < dim; ++e) {
204:             realSpaceDer[d] += invJ[e*dim+d]*basisDer[(q*numBasisFuncs+f)*dim+e];
205:           }
206:           fieldGrad[d] += realSpaceDer[d]*x[f];
207:         }
208:       }
209:       if (debug) {
210:         for(int d = 0; d < dim; ++d) {
211:           PetscPrintf(PETSC_COMM_SELF, "    fieldGrad[%d] %g\n", d, fieldGrad[d]);
212:         }
213:       }
214:       const PetscScalar funcVal = 0.0; //(*rhsFunc)(coords);
215:       for(int f = 0; f < numBasisFuncs; ++f) {
216:         /* Constant term: -f(x) */
217:         elemVec[f] -= basis[q*numBasisFuncs+f]*funcVal*quadWeights[q]*detJ;
218:         /* Linear term: -\Delta u */
219:         PetscScalar product = 0.0;
220:         for(int d = 0; d < dim; ++d) {
221:           realSpaceDer[d] = 0.0;
222:           for(int e = 0; e < dim; ++e) {
223:             realSpaceDer[d] += invJ[e*dim+d]*basisDer[(q*numBasisFuncs+f)*dim+e];
224:           }
225:           product += realSpaceDer[d]*fieldGrad[d];
226:         }
227:         elemVec[f] += product*quadWeights[q]*detJ;
228:         /* Nonlinear term: -\lambda e^{u} */
229:         elemVec[f] -= basis[q*numBasisFuncs+f]*0.0*PetscExpScalar(fieldVal)*quadWeights[q]*detJ;
230:       }
231:     }
232:     if (debug) {DMMeshPrintCellVector(c, "Residual", numBasisFuncs, elemVec);}
233:     DMMeshVecSetClosure(dm, F, c, elemVec, ADD_VALUES);
234:   }
235:   PetscLogFlops((cEnd-cStart)*numQuadPoints*numBasisFuncs*(dim*(dim*5+4)+14));
236:   PetscFree3(realSpaceDer,fieldGrad,elemVec);
237:   PetscFree4(coords,v0,J,invJ);

239:   PetscPrintf(PETSC_COMM_WORLD, "Residual:\n");
240:   for(int p = 0; p < user->numProcs; ++p) {
241:     if (p == user->rank) {VecView(F, PETSC_VIEWER_STDOUT_SELF);}
242:     PetscBarrier((PetscObject) dm);
243:   }
244:   return(0);
245: }


251: /*
252:  dm - The mesh
253:  X  - Local intput vector
254:  F  - Local output vector
255:  */
256: PetscErrorCode FormFunctionLocalBatch(DM dm, Vec X, Vec F, AppCtx *user)
257: {
258:   //PetscScalar    (*rhsFunc)(const PetscReal []) = user->rhsFunc;
259:   const PetscInt   debug         = user->debug;
260:   const PetscInt   dim           = user->dim;
261:   const PetscInt   numQuadPoints = user->q.numQuadPoints;
262:   const PetscReal *quadPoints    = user->q.quadPoints;
263:   const PetscReal *quadWeights   = user->q.quadWeights;
264:   const PetscInt   numBasisFuncs = user->q.numBasisFuncs;
265:   const PetscReal *basis         = user->q.basis;
266:   const PetscReal *basisDer      = user->q.basisDer;
267:   PetscReal       *coords, *v0, *J, *invJ, *detJ;
268:   PetscScalar     *realSpaceDer, *fieldGrad, *elemVec;
269:   PetscInt         cStart, cEnd;
270:   PetscErrorCode   ierr;

273:   VecSet(F, 0.0);
274:   //PetscMalloc3(dim,PetscScalar,&realSpaceDer,dim,PetscScalar,&fieldGrad,numBasisFuncs,PetscScalar,&elemVec);
275:   PetscMalloc4(dim,PetscReal,&coords,dim,PetscReal,&v0,dim*dim,PetscReal,&J,dim*dim,PetscReal,&invJ);
276:   DMMeshGetHeightStratum(dm, 0, &cStart, &cEnd);
277:   const PetscInt numCells = cEnd - cStart;
278:   PetscScalar   *u;

280:   PetscMalloc3(numCells*numBasisFuncs,PetscScalar,&u,numCells,PetscReal,&detJ,numCells*numBasisFuncs,PetscScalar,&elemVec);
281:   for(PetscInt c = cStart; c < cEnd; ++c) {
282:     const PetscScalar *x;

284:     DMMeshComputeCellGeometry(dm, c, v0, J, invJ, &detJ[c]);
285:     DMMeshVecGetClosure(dm, X, c, &x);

287:     for(int f = 0; f < numBasisFuncs; ++f) {
288:       u[c*numBasisFuncs+f] = x[f];
289:     }
290:   }
291:   IntegrateElementBatch(numCells, numBasisFuncs, u, detJ, numQuadPoints, quadPoints, quadWeights, basis, basisDer, elemVec);
292:   for(PetscInt c = cStart; c < cEnd; ++c) {
293:     DMMeshVecSetClosure(dm, F, c, &elemVec[c*numBasisFuncs], ADD_VALUES);
294:   }
295:   PetscFree3(u,detJ,elemVec);

297:   PetscLogFlops(0);
298:   //PetscFree3(realSpaceDer,fieldGrad,elemVec);
299:   PetscFree4(coords,v0,J,invJ);

301:   PetscPrintf(PETSC_COMM_WORLD, "Residual:\n");
302:   for(int p = 0; p < user->numProcs; ++p) {
303:     if (p == user->rank) {VecView(F, PETSC_VIEWER_STDOUT_SELF);}
304:     PetscBarrier((PetscObject) dm);
305:   }
306:   return(0);
307: }

311: /*
312:   dm  - The mesh
313:   X   - The local input vector
314:   Jac - The output matrix
315: */
316: PetscErrorCode FormJacobianLocal(DM dm, Vec X, Mat Jac, AppCtx *user)
317: {
318:   const PetscInt   debug         = user->debug;
319:   const PetscInt   dim           = user->dim;
320:   const PetscInt   numQuadPoints = user->q.numQuadPoints;
321:   const PetscReal *quadWeights   = user->q.quadWeights;
322:   const PetscInt   numBasisFuncs = user->q.numBasisFuncs;
323:   const PetscReal *basis         = user->q.basis;
324:   const PetscReal *basisDer      = user->q.basisDer;
325:   PetscReal       *v0, *J, *invJ, detJ;
326:   PetscScalar     *realSpaceTestDer, *realSpaceBasisDer, *elemMat;
327:   PetscInt         cStart, cEnd;
328:   PetscErrorCode   ierr;

331:   MatZeroEntries(Jac);
332:   PetscMalloc3(dim,PetscScalar,&realSpaceTestDer,dim,PetscScalar,&realSpaceBasisDer,numBasisFuncs*numBasisFuncs,PetscScalar,&elemMat);
333:   PetscMalloc3(dim,PetscReal,&v0,dim*dim,PetscReal,&J,dim*dim,PetscReal,&invJ);
334:   DMMeshGetHeightStratum(dm, 0, &cStart, &cEnd);
335:   for(PetscInt c = cStart; c < cEnd; ++c) {
336:     const PetscScalar *x;

338:     PetscMemzero(elemMat, numBasisFuncs*numBasisFuncs * sizeof(PetscScalar));
339:     DMMeshComputeCellGeometry(dm, c, v0, J, invJ, &detJ);
340:     if (detJ <= 0.0) {SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid determinant %g for element %d", detJ, c);}
341:     DMMeshVecGetClosure(dm, X, c, &x);

343:     for(int q = 0; q < numQuadPoints; ++q) {
344:       PetscScalar fieldVal = 0.0;

346:       for(int f = 0; f < numBasisFuncs; ++f) {
347:         fieldVal += x[f]*basis[q*numBasisFuncs+f];
348:       }
349:       for(int f = 0; f < numBasisFuncs; ++f) {
350:         for(int d = 0; d < dim; ++d) {
351:           realSpaceTestDer[d] = 0.0;
352:           for(int e = 0; e < dim; ++e) {
353:             realSpaceTestDer[d] += invJ[e*dim+d]*basisDer[(q*numBasisFuncs+f)*dim+e];
354:           }
355:         }
356:         for(int g = 0; g < numBasisFuncs; ++g) {
357:           for(int d = 0; d < dim; ++d) {
358:             realSpaceBasisDer[d] = 0.0;
359:             for(int e = 0; e < dim; ++e) {
360:               realSpaceBasisDer[d] += invJ[e*dim+d]*basisDer[(q*numBasisFuncs+g)*dim+e];
361:             }
362:           }
363:           /* Linear term: -\Delta u */
364:           PetscScalar product = 0.0;
365:           for(int d = 0; d < dim; ++d) product += realSpaceTestDer[d]*realSpaceBasisDer[d];
366:           elemMat[f*numBasisFuncs+g] += product*quadWeights[q]*detJ;
367:           /* Nonlinear term: -\lambda e^{u} */
368:           elemMat[f*numBasisFuncs+g] -= basis[q*numBasisFuncs+f]*basis[q*numBasisFuncs+g]*0.0*PetscExpScalar(fieldVal)*quadWeights[q]*detJ;
369:         }
370:       }
371:     }
372:     if (debug) {DMMeshPrintCellMatrix(c, "Jacobian", numBasisFuncs, numBasisFuncs, elemMat);}
373:     DMMeshMatSetClosure(dm, Jac, c, elemMat, ADD_VALUES);
374:   }
375:   PetscLogFlops((cEnd-cStart)*numQuadPoints*numBasisFuncs*(dim*(dim*5+4)+14));
376:   PetscFree3(realSpaceTestDer,realSpaceBasisDer,elemMat);
377:   PetscFree3(v0,J,invJ);

379:   /* Assemble matrix, using the 2-step process:
380:        MatAssemblyBegin(), MatAssemblyEnd(). */
381:   MatAssemblyBegin(Jac, MAT_FINAL_ASSEMBLY);
382:   MatAssemblyEnd(Jac, MAT_FINAL_ASSEMBLY);
383:   /* Tell the matrix we will never add a new nonzero location to the
384:      matrix. If we do, it will generate an error. */
385:   MatSetOption(Jac, MAT_NEW_NONZERO_LOCATION_ERR, PETSC_TRUE);
386:   return(0);
387: }

391: /*
392:   dm  - The mesh
393:   X   - The local input vector
394:   Jac - The output matrix
395: */
396: PetscErrorCode FormJacobianLocalBatch(DM dm, Vec X, Mat Jac, AppCtx *user)
397: {
398:   const PetscInt   debug         = user->debug;
399:   const PetscInt   dim           = user->dim;
400:   const PetscInt   numQuadPoints = user->q.numQuadPoints;
401:   const PetscReal *quadWeights   = user->q.quadWeights;
402:   const PetscInt   numBasisFuncs = user->q.numBasisFuncs;
403:   const PetscReal *basis         = user->q.basis;
404:   const PetscReal *basisDer      = user->q.basisDer;
405:   PetscReal       *v0, *J, *invJ, detJ;
406:   PetscScalar     *realSpaceTestDer, *realSpaceBasisDer, *elemMat;
407:   PetscInt         cStart, cEnd;
408:   PetscErrorCode   ierr;

411:   MatZeroEntries(Jac);
412:   PetscMalloc3(dim,PetscScalar,&realSpaceTestDer,dim,PetscScalar,&realSpaceBasisDer,numBasisFuncs*numBasisFuncs,PetscScalar,&elemMat);
413:   PetscMalloc3(dim,PetscReal,&v0,dim*dim,PetscReal,&J,dim*dim,PetscReal,&invJ);
414:   DMMeshGetHeightStratum(dm, 0, &cStart, &cEnd);
415:   for(PetscInt c = cStart; c < cEnd; ++c) {
416:     const PetscScalar *x;

418:     PetscMemzero(elemMat, numBasisFuncs*numBasisFuncs * sizeof(PetscScalar));
419:     DMMeshComputeCellGeometry(dm, c, v0, J, invJ, &detJ);
420:     if (detJ <= 0.0) {SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid determinant %g for element %d", detJ, c);}
421:     DMMeshVecGetClosure(dm, X, c, &x);

423:     for(int q = 0; q < numQuadPoints; ++q) {
424:       PetscScalar fieldVal = 0.0;

426:       for(int f = 0; f < numBasisFuncs; ++f) {
427:         fieldVal += x[f]*basis[q*numBasisFuncs+f];
428:       }
429:       for(int f = 0; f < numBasisFuncs; ++f) {
430:         for(int d = 0; d < dim; ++d) {
431:           realSpaceTestDer[d] = 0.0;
432:           for(int e = 0; e < dim; ++e) {
433:             realSpaceTestDer[d] += invJ[e*dim+d]*basisDer[(q*numBasisFuncs+f)*dim+e];
434:           }
435:         }
436:         for(int g = 0; g < numBasisFuncs; ++g) {
437:           for(int d = 0; d < dim; ++d) {
438:             realSpaceBasisDer[d] = 0.0;
439:             for(int e = 0; e < dim; ++e) {
440:               realSpaceBasisDer[d] += invJ[e*dim+d]*basisDer[(q*numBasisFuncs+g)*dim+e];
441:             }
442:           }
443:           /* Linear term: -\Delta u */
444:           PetscScalar product = 0.0;
445:           for(int d = 0; d < dim; ++d) product += realSpaceTestDer[d]*realSpaceBasisDer[d];
446:           elemMat[f*numBasisFuncs+g] += product*quadWeights[q]*detJ;
447:           /* Nonlinear term: -\lambda e^{u} */
448:           elemMat[f*numBasisFuncs+g] -= basis[q*numBasisFuncs+f]*basis[q*numBasisFuncs+g]*0.0*PetscExpScalar(fieldVal)*quadWeights[q]*detJ;
449:         }
450:       }
451:     }
452:     if (debug) {DMMeshPrintCellMatrix(c, "Jacobian", numBasisFuncs, numBasisFuncs, elemMat);}
453:     DMMeshMatSetClosure(dm, Jac, c, elemMat, ADD_VALUES);
454:   }
455:   PetscLogFlops((cEnd-cStart)*numQuadPoints*numBasisFuncs*(dim*(dim*5+4)+14));
456:   PetscFree3(realSpaceTestDer,realSpaceBasisDer,elemMat);
457:   PetscFree3(v0,J,invJ);

459:   /* Assemble matrix, using the 2-step process:
460:        MatAssemblyBegin(), MatAssemblyEnd(). */
461:   MatAssemblyBegin(Jac, MAT_FINAL_ASSEMBLY);
462:   MatAssemblyEnd(Jac, MAT_FINAL_ASSEMBLY);
463:   /* Tell the matrix we will never add a new nonzero location to the
464:      matrix. If we do, it will generate an error. */
465:   MatSetOption(Jac, MAT_NEW_NONZERO_LOCATION_ERR, PETSC_TRUE);
466:   return(0);
467: }

471: int main(int argc, char **argv)
472: {
473:   DM             dm;
474:   SNES           snes;
475:   AppCtx         user;

478:   PetscInitialize(&argc, &argv, PETSC_NULL, help);
479: #ifndef PETSC_HAVE_CUDA
480:   SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_SUP, "This example requires CUDA support.");
481: #endif
482:   ProcessOptions(PETSC_COMM_WORLD, &user);
483:   CreateMesh(PETSC_COMM_WORLD, &user, &dm);
484:   SetupQuadrature(&user);
485:   SetupSection(dm, &user);

487:   SNESCreate(PETSC_COMM_WORLD, &snes);
488:   if (user.computeFunction) {
489:     Vec            X, F;

491:     DMGetGlobalVector(dm, &X);
492:     DMGetGlobalVector(dm, &F);
493:     if (user.batch) {
494:       DMMeshSetLocalFunction(dm, (PetscErrorCode (*)(DM, Vec, Vec, void*)) FormFunctionLocalBatch);
495:     } else {
496:       DMMeshSetLocalFunction(dm, (PetscErrorCode (*)(DM, Vec, Vec, void*)) FormFunctionLocal);
497:     }
498:     SNESMeshFormFunction(snes, X, F, &user);
499:     DMRestoreGlobalVector(dm, &X);
500:     DMRestoreGlobalVector(dm, &F);
501:   }
502:   if (user.computeJacobian) {
503:     Vec            X;
504:     Mat            J;
505:     MatStructure   flag;

507:     DMGetGlobalVector(dm, &X);
508:     DMGetMatrix(dm, MATAIJ, &J);
509:     if (user.batch) {
510:       DMMeshSetLocalJacobian(dm, (PetscErrorCode (*)(DM, Vec, Mat, void*)) FormJacobianLocalBatch);
511:     } else {
512:       DMMeshSetLocalJacobian(dm, (PetscErrorCode (*)(DM, Vec, Mat, void*)) FormJacobianLocal);
513:     }
514:     SNESMeshFormJacobian(snes, X, &J, &J, &flag, &user);
515:     MatDestroy(&J);
516:     DMRestoreGlobalVector(dm, &X);
517:   }
518:   SNESDestroy(&snes);
519:   DMDestroy(&dm);
520:   PetscFinalize();
521:   return 0;
522: }