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();

 51:   PetscFunctionReturn(PETSC_SUCCESS);
 52: }

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

302: /*TEST

304:   build:
305:     requires: pragmatic

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

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

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

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

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

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

359: TEST*/