Actual source code: ex19.c

  1: static char help[] = "Tests mesh adaptation with DMPlex and pragmatic.\n";

  3: #include <petsc/private/dmpleximpl.h>

  5: #include <petscsnes.h>

  7: typedef struct {
  8:   PetscInt  Nr;         /* The number of refinement passes */
  9:   PetscInt  metOpt;     /* Different choices of metric */
 10:   PetscReal hmax, hmin; /* Max and min sizes prescribed by the metric */
 11:   PetscBool doL2;       /* Test L2 projection */
 12: } AppCtx;

 14: /*
 15: Classic hyperbolic sensor function for testing multi-scale anisotropic mesh adaptation:

 17:   f:[-1, 1]x[-1, 1] \to R,
 18:     f(x, y) = sin(50xy)/100 if |xy| > 2\pi/50 else sin(50xy)

 20: (mapped to have domain [0,1] x [0,1] in this case).
 21: */
 22: static PetscErrorCode sensor(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar u[], void *ctx)
 23: {
 24:   const PetscReal xref = 2. * x[0] - 1.;
 25:   const PetscReal yref = 2. * x[1] - 1.;
 26:   const PetscReal xy   = xref * yref;

 29:   u[0] = PetscSinReal(50. * xy);
 30:   if (PetscAbsReal(xy) > 2. * PETSC_PI / 50.) u[0] *= 0.01;
 31:   return 0;
 32: }

 34: static PetscErrorCode ProcessOptions(MPI_Comm comm, AppCtx *options)
 35: {
 36:   options->Nr     = 1;
 37:   options->metOpt = 1;
 38:   options->hmin   = 0.05;
 39:   options->hmax   = 0.5;
 40:   options->doL2   = PETSC_FALSE;

 42:   PetscOptionsBegin(comm, "", "Meshing Adaptation Options", "DMPLEX");
 43:   PetscOptionsBoundedInt("-Nr", "Numberof refinement passes", "ex19.c", options->Nr, &options->Nr, NULL, 1);
 44:   PetscOptionsBoundedInt("-met", "Different choices of metric", "ex19.c", options->metOpt, &options->metOpt, NULL, 0);
 45:   PetscOptionsReal("-hmax", "Max size prescribed by the metric", "ex19.c", options->hmax, &options->hmax, NULL);
 46:   PetscOptionsReal("-hmin", "Min size prescribed by the metric", "ex19.c", options->hmin, &options->hmin, NULL);
 47:   PetscOptionsBool("-do_L2", "Test L2 projection", "ex19.c", options->doL2, &options->doL2, NULL);
 48:   PetscOptionsEnd();

 50:   return 0;
 51: }

 53: static PetscErrorCode CreateMesh(MPI_Comm comm, DM *dm)
 54: {
 55:   DMCreate(comm, dm);
 56:   DMSetType(*dm, DMPLEX);
 57:   DMSetFromOptions(*dm);
 58:   PetscObjectSetName((PetscObject)*dm, "DMinit");
 59:   DMViewFromOptions(*dm, NULL, "-init_dm_view");
 60:   return 0;
 61: }

 63: static PetscErrorCode ComputeMetricSensor(DM dm, AppCtx *user, Vec *metric)
 64: {
 65:   PetscSimplePointFunc funcs[1] = {sensor};
 66:   DM                   dmSensor, dmGrad, dmHess, dmDet;
 67:   PetscFE              fe;
 68:   Vec                  f, g, H, determinant;
 69:   PetscBool            simplex;
 70:   PetscInt             dim;

 72:   DMGetDimension(dm, &dim);
 73:   DMPlexIsSimplex(dm, &simplex);

 75:   DMClone(dm, &dmSensor);
 76:   PetscFECreateLagrange(PETSC_COMM_SELF, dim, 1, simplex, 1, -1, &fe);
 77:   DMSetField(dmSensor, 0, NULL, (PetscObject)fe);
 78:   PetscFEDestroy(&fe);
 79:   DMCreateDS(dmSensor);
 80:   DMCreateLocalVector(dmSensor, &f);
 81:   DMProjectFunctionLocal(dmSensor, 0., funcs, NULL, INSERT_VALUES, f);
 82:   VecViewFromOptions(f, NULL, "-sensor_view");

 84:   // Recover the gradient of the sensor function
 85:   DMClone(dm, &dmGrad);
 86:   PetscFECreateLagrange(PETSC_COMM_SELF, dim, dim, simplex, 1, -1, &fe);
 87:   DMSetField(dmGrad, 0, NULL, (PetscObject)fe);
 88:   PetscFEDestroy(&fe);
 89:   DMCreateDS(dmGrad);
 90:   DMCreateLocalVector(dmGrad, &g);
 91:   DMPlexComputeGradientClementInterpolant(dmSensor, f, g);
 92:   VecDestroy(&f);
 93:   VecViewFromOptions(g, NULL, "-gradient_view");

 95:   // Recover the Hessian of the sensor function
 96:   DMClone(dm, &dmHess);
 97:   DMPlexMetricCreate(dmHess, 0, &H);
 98:   DMPlexComputeGradientClementInterpolant(dmGrad, g, H);
 99:   VecDestroy(&g);
100:   VecViewFromOptions(H, NULL, "-hessian_view");

102:   // Obtain a metric by Lp normalization
103:   DMPlexMetricCreate(dm, 0, metric);
104:   DMPlexMetricDeterminantCreate(dm, 0, &determinant, &dmDet);
105:   DMPlexMetricNormalize(dmHess, H, PETSC_TRUE, PETSC_TRUE, *metric, determinant);
106:   VecDestroy(&determinant);
107:   DMDestroy(&dmDet);
108:   VecDestroy(&H);
109:   DMDestroy(&dmHess);
110:   DMDestroy(&dmGrad);
111:   DMDestroy(&dmSensor);
112:   return 0;
113: }

115: static PetscErrorCode ComputeMetric(DM dm, AppCtx *user, Vec *metric)
116: {
117:   PetscReal lambda = 1 / (user->hmax * user->hmax);

120:   if (user->metOpt == 0) {
121:     /* Specify a uniform, isotropic metric */
122:     DMPlexMetricCreateUniform(dm, 0, lambda, metric);
123:   } else if (user->metOpt == 3) {
124:     ComputeMetricSensor(dm, user, metric);
125:   } else {
126:     DM                 cdm;
127:     Vec                coordinates;
128:     const PetscScalar *coords;
129:     PetscScalar       *met;
130:     PetscReal          h;
131:     PetscInt           dim, i, j, vStart, vEnd, v;

133:     DMPlexMetricCreate(dm, 0, metric);
134:     DMGetDimension(dm, &dim);
135:     DMGetCoordinateDM(dm, &cdm);
136:     DMGetCoordinatesLocal(dm, &coordinates);
137:     VecGetArrayRead(coordinates, &coords);
138:     VecGetArray(*metric, &met);
139:     DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);
140:     for (v = vStart; v < vEnd; ++v) {
141:       PetscScalar *vcoords;
142:       PetscScalar *pmet;

144:       DMPlexPointLocalRead(cdm, v, coords, &vcoords);
145:       switch (user->metOpt) {
146:       case 1:
147:         h = user->hmax - (user->hmax - user->hmin) * PetscRealPart(vcoords[0]);
148:         break;
149:       case 2:
150:         h = user->hmax * PetscAbsReal(((PetscReal)1.0) - PetscExpReal(-PetscAbsScalar(vcoords[0] - (PetscReal)0.5))) + user->hmin;
151:         break;
152:       default:
153:         SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "metOpt = 0, 1, 2 or 3, cannot be %d", user->metOpt);
154:       }
155:       DMPlexPointLocalRef(dm, v, met, &pmet);
156:       for (i = 0; i < dim; ++i) {
157:         for (j = 0; j < dim; ++j) {
158:           if (i == j) {
159:             if (i == 0) pmet[i * dim + j] = 1 / (h * h);
160:             else pmet[i * dim + j] = lambda;
161:           } else pmet[i * dim + j] = 0.0;
162:         }
163:       }
164:     }
165:     VecRestoreArray(*metric, &met);
166:     VecRestoreArrayRead(coordinates, &coords);
167:   }
168:   return 0;
169: }

171: static PetscErrorCode linear(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx)
172: {
173:   u[0] = x[0] + x[1];
174:   return 0;
175: }

177: static PetscErrorCode TestL2Projection(DM dm, DM dma, AppCtx *user)
178: {
179:   PetscErrorCode (*funcs[1])(PetscInt, PetscReal, const PetscReal[], PetscInt, PetscScalar *, void *) = {linear};
180:   DM        dmProj, dmaProj;
181:   PetscFE   fe;
182:   KSP       ksp;
183:   Mat       Interp, mass, mass2;
184:   Vec       u, ua, scaling, rhs, uproj;
185:   PetscReal error;
186:   PetscBool simplex;
187:   PetscInt  dim;

190:   DMGetDimension(dm, &dim);
191:   DMPlexIsSimplex(dm, &simplex);

193:   DMClone(dm, &dmProj);
194:   PetscFECreateDefault(PETSC_COMM_SELF, dim, 1, simplex, NULL, -1, &fe);
195:   DMSetField(dmProj, 0, NULL, (PetscObject)fe);
196:   PetscFEDestroy(&fe);
197:   DMCreateDS(dmProj);

199:   DMClone(dma, &dmaProj);
200:   PetscFECreateDefault(PETSC_COMM_SELF, dim, 1, simplex, NULL, -1, &fe);
201:   DMSetField(dmaProj, 0, NULL, (PetscObject)fe);
202:   PetscFEDestroy(&fe);
203:   DMCreateDS(dmaProj);

205:   DMGetGlobalVector(dmProj, &u);
206:   DMGetGlobalVector(dmaProj, &ua);
207:   DMGetGlobalVector(dmaProj, &rhs);
208:   DMGetGlobalVector(dmaProj, &uproj);

210:   // Interpolate onto original mesh using dual basis
211:   DMProjectFunction(dmProj, 0.0, funcs, NULL, INSERT_VALUES, u);
212:   PetscObjectSetName((PetscObject)u, "Original");
213:   VecViewFromOptions(u, NULL, "-orig_vec_view");
214:   DMComputeL2Diff(dmProj, 0.0, funcs, NULL, u, &error);
215:   PetscPrintf(PETSC_COMM_WORLD, "Original L2 Error: %g\n", (double)error);
216:   // Interpolate onto NEW mesh using dual basis
217:   DMProjectFunction(dmaProj, 0.0, funcs, NULL, INSERT_VALUES, ua);
218:   PetscObjectSetName((PetscObject)ua, "Adapted");
219:   VecViewFromOptions(ua, NULL, "-adapt_vec_view");
220:   DMComputeL2Diff(dmaProj, 0.0, funcs, NULL, ua, &error);
221:   PetscPrintf(PETSC_COMM_WORLD, "Adapted L2 Error: %g\n", (double)error);
222:   // Interpolate between meshes using interpolation matrix
223:   DMCreateInterpolation(dmProj, dmaProj, &Interp, &scaling);
224:   MatInterpolate(Interp, u, ua);
225:   MatDestroy(&Interp);
226:   VecDestroy(&scaling);
227:   PetscObjectSetName((PetscObject)ua, "Interpolation");
228:   VecViewFromOptions(ua, NULL, "-interp_vec_view");
229:   DMComputeL2Diff(dmaProj, 0.0, funcs, NULL, ua, &error);
230:   PetscPrintf(PETSC_COMM_WORLD, "Interpolated L2 Error: %g\n", (double)error);
231:   // L2 projection
232:   DMCreateMassMatrix(dmaProj, dmaProj, &mass);
233:   MatViewFromOptions(mass, NULL, "-mass_mat_view");
234:   KSPCreate(PETSC_COMM_WORLD, &ksp);
235:   KSPSetOperators(ksp, mass, mass);
236:   KSPSetFromOptions(ksp);
237:   //   Compute rhs as M f, could also directly project the analytic function but we might not have it
238:   DMCreateMassMatrix(dmProj, dmaProj, &mass2);
239:   MatMult(mass2, u, rhs);
240:   MatDestroy(&mass2);
241:   KSPSolve(ksp, rhs, uproj);
242:   PetscObjectSetName((PetscObject)uproj, "L_2 Projection");
243:   VecViewFromOptions(uproj, NULL, "-proj_vec_view");
244:   DMComputeL2Diff(dmaProj, 0.0, funcs, NULL, uproj, &error);
245:   PetscPrintf(PETSC_COMM_WORLD, "Projected L2 Error: %g\n", (double)error);
246:   KSPDestroy(&ksp);
247:   MatDestroy(&mass);
248:   DMRestoreGlobalVector(dmProj, &u);
249:   DMRestoreGlobalVector(dmaProj, &ua);
250:   DMRestoreGlobalVector(dmaProj, &rhs);
251:   DMRestoreGlobalVector(dmaProj, &uproj);
252:   DMDestroy(&dmProj);
253:   DMDestroy(&dmaProj);
254:   return 0;
255: }

257: int main(int argc, char *argv[])
258: {
259:   DM       dm;
260:   AppCtx   user; /* user-defined work context */
261:   MPI_Comm comm;
262:   DM       dma, odm;
263:   Vec      metric;
264:   PetscInt r;

267:   PetscInitialize(&argc, &argv, NULL, help);
268:   comm = PETSC_COMM_WORLD;
269:   ProcessOptions(comm, &user);
270:   CreateMesh(comm, &dm);

272:   odm = dm;
273:   DMPlexDistributeOverlap(odm, 1, NULL, &dm);
274:   if (!dm) {
275:     dm = odm;
276:   } else DMDestroy(&odm);

278:   for (r = 0; r < user.Nr; ++r) {
279:     DMLabel label;

281:     ComputeMetric(dm, &user, &metric);
282:     DMGetLabel(dm, "marker", &label);
283:     DMAdaptMetric(dm, metric, label, NULL, &dma);
284:     VecDestroy(&metric);
285:     PetscObjectSetName((PetscObject)dma, "DMadapt");
286:     PetscObjectSetOptionsPrefix((PetscObject)dma, "adapt_");
287:     DMViewFromOptions(dma, NULL, "-dm_view");
288:     if (user.doL2) TestL2Projection(dm, dma, &user);
289:     DMDestroy(&dm);
290:     dm = dma;
291:   }
292:   PetscObjectSetOptionsPrefix((PetscObject)dm, "final_");
293:   DMViewFromOptions(dm, NULL, "-dm_view");
294:   DMDestroy(&dm);
295:   PetscFinalize();
296:   return 0;
297: }

299: /*TEST

301:   build:
302:     requires: pragmatic

304:   testset:
305:     args: -dm_plex_box_faces 4,4,4 -dm_adaptor pragmatic -met 2 -init_dm_view -adapt_dm_view -dm_adaptor pragmatic

307:     test:
308:       suffix: 2d
309:       args: -dm_plex_separate_marker 0
310:     test:
311:       suffix: 2d_sep
312:       args: -dm_plex_separate_marker 1
313:     test:
314:       suffix: 3d
315:       args: -dm_plex_dim 3

317:   # Pragmatic hangs for simple partitioner
318:   testset:
319:     requires: parmetis
320:     args: -dm_plex_box_faces 2,2 -petscpartitioner_type parmetis -met 2 -init_dm_view -adapt_dm_view -dm_adaptor pragmatic

322:     test:
323:       suffix: 2d_parmetis_np2
324:       nsize: 2
325:     test:
326:       suffix: 2d_parmetis_np4
327:       nsize: 4

329:   test:
330:     requires: parmetis
331:     suffix: 3d_parmetis_met0
332:     nsize: 2
333:     args: -dm_plex_dim 3 -dm_plex_box_faces 9,9,9 -dm_adaptor pragmatic -petscpartitioner_type parmetis \
334:           -met 0 -hmin 0.01 -hmax 0.03 -init_dm_view -adapt_dm_view -dm_adaptor pragmatic
335:   test:
336:     requires: parmetis
337:     suffix: 3d_parmetis_met2
338:     nsize: 2
339:     args: -dm_plex_box_faces 19,19 -dm_adaptor pragmatic -petscpartitioner_type parmetis \
340:           -met 2 -hmax 0.5 -hmin 0.001 -init_dm_view -adapt_dm_view -dm_adaptor pragmatic
341:   test:
342:     suffix: proj2
343:     args: -dm_plex_box_faces 2,2 -dm_plex_hash_location -dm_adaptor pragmatic -init_dm_view -adapt_dm_view -do_L2 \
344:           -petscspace_degree 1 -petscfe_default_quadrature_order 1 -pc_type lu -dm_adaptor pragmatic
345:   test:
346:     suffix: proj4
347:     args: -dm_plex_box_faces 4,4 -dm_plex_hash_location -dm_adaptor pragmatic -init_dm_view -adapt_dm_view -do_L2 \
348:           -petscspace_degree 1 -petscfe_default_quadrature_order 1 -pc_type lu -dm_adaptor pragmatic

350:   test:
351:     suffix: 2d_met3
352:     args: -dm_plex_box_faces 9,9 -met 3 -dm_adaptor pragmatic -init_dm_view -adapt_dm_view \
353:           -dm_plex_metric_h_min 1.e-10 -dm_plex_metric_h_max 1.0e-01 -dm_plex_metric_a_max 1.0e+05 -dm_plex_metric_p 1.0 \
354:             -dm_plex_metric_target_complexity 10000.0 -dm_adaptor pragmatic

356: TEST*/