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;

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

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

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

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

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

 73:   PetscFunctionBegin;
 74:   PetscCall(DMGetDimension(dm, &dim));
 75:   PetscCall(DMPlexIsSimplex(dm, &simplex));

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

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

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

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

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

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

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

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

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

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

191:   PetscFunctionBeginUser;
192:   PetscCall(DMGetDimension(dm, &dim));
193:   PetscCall(DMPlexIsSimplex(dm, &simplex));

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

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

207:   PetscCall(DMGetGlobalVector(dmProj, &u));
208:   PetscCall(DMGetGlobalVector(dmaProj, &ua));
209:   PetscCall(DMGetGlobalVector(dmaProj, &rhs));
210:   PetscCall(DMGetGlobalVector(dmaProj, &uproj));

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

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

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

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

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

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

301: /*TEST

303:   build:
304:     requires: pragmatic

306:   testset:
307:     args: -dm_plex_box_faces 4,4,4 -dm_adaptor pragmatic -met 2 -init_dm_view -adapt_dm_view -dm_adaptor pragmatic

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

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

324:     test:
325:       suffix: 2d_parmetis_np2
326:       nsize: 2
327:     test:
328:       suffix: 2d_parmetis_np4
329:       nsize: 4

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

352:   test:
353:     suffix: 2d_met3
354:     args: -dm_plex_box_faces 9,9 -met 3 -dm_adaptor pragmatic -init_dm_view -adapt_dm_view \
355:           -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 \
356:             -dm_plex_metric_target_complexity 10000.0 -dm_adaptor pragmatic

358: TEST*/