1 static char help[] = "Tests mesh adaptation with DMPlex and pragmatic.\n"; 2 3 #include <petsc/private/dmpleximpl.h> 4 5 #include <petscsnes.h> 6 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; 13 14 /* 15 Classic hyperbolic sensor function for testing multi-scale anisotropic mesh adaptation: 16 17 f:[-1, 1]x[-1, 1] \to R, 18 f(x, y) = sin(50xy)/100 if |xy| > 2\pi/50 else sin(50xy) 19 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; 27 28 PetscFunctionBeginUser; 29 u[0] = PetscSinReal(50.*xy); 30 if (PetscAbsReal(xy) > 2.*PETSC_PI/50.) u[0] *= 0.01; 31 PetscFunctionReturn(0); 32 } 33 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; 42 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 51 PetscFunctionReturn(0); 52 } 53 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(0); 63 } 64 65 static PetscErrorCode ComputeMetricSensor(DM dm, AppCtx *user, Vec *metric) 66 { 67 PetscSimplePointFunc funcs[1] = {sensor}; 68 DM dmSensor, dmGrad, dmHess; 69 PetscFE fe; 70 Vec f, g, H; 71 PetscBool simplex; 72 PetscInt dim; 73 74 PetscFunctionBegin; 75 PetscCall(DMGetDimension(dm, &dim)); 76 PetscCall(DMPlexIsSimplex(dm, &simplex)); 77 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")); 86 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")); 97 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")); 104 105 // Obtain a metric by Lp normalization 106 PetscCall(DMPlexMetricNormalize(dmHess, H, PETSC_TRUE, PETSC_TRUE, metric)); 107 PetscCall(VecDestroy(&H)); 108 PetscCall(DMDestroy(&dmHess)); 109 PetscCall(DMDestroy(&dmGrad)); 110 PetscCall(DMDestroy(&dmSensor)); 111 PetscFunctionReturn(0); 112 } 113 114 static PetscErrorCode ComputeMetric(DM dm, AppCtx *user, Vec *metric) 115 { 116 PetscReal lambda = 1/(user->hmax*user->hmax); 117 118 PetscFunctionBeginUser; 119 if (user->metOpt == 0) { 120 /* Specify a uniform, isotropic metric */ 121 PetscCall(DMPlexMetricCreateUniform(dm, 0, lambda, metric)); 122 } else if (user->metOpt == 3) { 123 PetscCall(ComputeMetricSensor(dm, user, metric)); 124 } else { 125 DM cdm; 126 Vec coordinates; 127 const PetscScalar *coords; 128 PetscScalar *met; 129 PetscReal h; 130 PetscInt dim, i, j, vStart, vEnd, v; 131 132 PetscCall(DMPlexMetricCreate(dm, 0, metric)); 133 PetscCall(DMGetDimension(dm, &dim)); 134 PetscCall(DMGetCoordinateDM(dm, &cdm)); 135 PetscCall(DMGetCoordinatesLocal(dm, &coordinates)); 136 PetscCall(VecGetArrayRead(coordinates, &coords)); 137 PetscCall(VecGetArray(*metric, &met)); 138 PetscCall(DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd)); 139 for (v = vStart; v < vEnd; ++v) { 140 PetscScalar *vcoords; 141 PetscScalar *pmet; 142 143 PetscCall(DMPlexPointLocalRead(cdm, v, coords, &vcoords)); 144 switch (user->metOpt) { 145 case 1: 146 h = user->hmax - (user->hmax-user->hmin)*PetscRealPart(vcoords[0]); 147 break; 148 case 2: 149 h = user->hmax*PetscAbsReal(((PetscReal) 1.0)-PetscExpReal(-PetscAbsScalar(vcoords[0]-(PetscReal)0.5))) + user->hmin; 150 break; 151 default: 152 SETERRQ(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_WRONG, "metOpt = 0, 1, 2 or 3, cannot be %d", user->metOpt); 153 } 154 PetscCall(DMPlexPointLocalRef(dm, v, met, &pmet)); 155 for (i = 0; i < dim; ++i) { 156 for (j = 0; j < dim; ++j) { 157 if (i == j) { 158 if (i == 0) pmet[i*dim+j] = 1/(h*h); 159 else pmet[i*dim+j] = lambda; 160 } else pmet[i*dim+j] = 0.0; 161 } 162 } 163 } 164 PetscCall(VecRestoreArray(*metric, &met)); 165 PetscCall(VecRestoreArrayRead(coordinates, &coords)); 166 } 167 PetscFunctionReturn(0); 168 } 169 170 static PetscErrorCode linear(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx) 171 { 172 u[0] = x[0] + x[1]; 173 return 0; 174 } 175 176 static PetscErrorCode TestL2Projection(DM dm, DM dma, AppCtx *user) 177 { 178 PetscErrorCode (*funcs[1])(PetscInt, PetscReal, const PetscReal [], PetscInt, PetscScalar *, void *) = {linear}; 179 DM dmProj, dmaProj; 180 PetscFE fe; 181 KSP ksp; 182 Mat Interp, mass, mass2; 183 Vec u, ua, scaling, rhs, uproj; 184 PetscReal error; 185 PetscBool simplex; 186 PetscInt dim; 187 188 PetscFunctionBeginUser; 189 PetscCall(DMGetDimension(dm, &dim)); 190 PetscCall(DMPlexIsSimplex(dm, &simplex)); 191 192 PetscCall(DMClone(dm, &dmProj)); 193 PetscCall(PetscFECreateDefault(PETSC_COMM_SELF, dim, 1, simplex, NULL, -1, &fe)); 194 PetscCall(DMSetField(dmProj, 0, NULL, (PetscObject) fe)); 195 PetscCall(PetscFEDestroy(&fe)); 196 PetscCall(DMCreateDS(dmProj)); 197 198 PetscCall(DMClone(dma, &dmaProj)); 199 PetscCall(PetscFECreateDefault(PETSC_COMM_SELF, dim, 1, simplex, NULL, -1, &fe)); 200 PetscCall(DMSetField(dmaProj, 0, NULL, (PetscObject) fe)); 201 PetscCall(PetscFEDestroy(&fe)); 202 PetscCall(DMCreateDS(dmaProj)); 203 204 PetscCall(DMGetGlobalVector(dmProj, &u)); 205 PetscCall(DMGetGlobalVector(dmaProj, &ua)); 206 PetscCall(DMGetGlobalVector(dmaProj, &rhs)); 207 PetscCall(DMGetGlobalVector(dmaProj, &uproj)); 208 209 // Interpolate onto original mesh using dual basis 210 PetscCall(DMProjectFunction(dmProj, 0.0, funcs, NULL, INSERT_VALUES, u)); 211 PetscCall(PetscObjectSetName((PetscObject) u, "Original")); 212 PetscCall(VecViewFromOptions(u, NULL, "-orig_vec_view")); 213 PetscCall(DMComputeL2Diff(dmProj, 0.0, funcs, NULL, u, &error)); 214 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Original L2 Error: %g\n", (double) error)); 215 // Interpolate onto NEW mesh using dual basis 216 PetscCall(DMProjectFunction(dmaProj, 0.0, funcs, NULL, INSERT_VALUES, ua)); 217 PetscCall(PetscObjectSetName((PetscObject) ua, "Adapted")); 218 PetscCall(VecViewFromOptions(ua, NULL, "-adapt_vec_view")); 219 PetscCall(DMComputeL2Diff(dmaProj, 0.0, funcs, NULL, ua, &error)); 220 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Adapted L2 Error: %g\n", (double) error)); 221 // Interpolate between meshes using interpolation matrix 222 PetscCall(DMCreateInterpolation(dmProj, dmaProj, &Interp, &scaling)); 223 PetscCall(MatInterpolate(Interp, u, ua)); 224 PetscCall(MatDestroy(&Interp)); 225 PetscCall(VecDestroy(&scaling)); 226 PetscCall(PetscObjectSetName((PetscObject) ua, "Interpolation")); 227 PetscCall(VecViewFromOptions(ua, NULL, "-interp_vec_view")); 228 PetscCall(DMComputeL2Diff(dmaProj, 0.0, funcs, NULL, ua, &error)); 229 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Interpolated L2 Error: %g\n", (double) error)); 230 // L2 projection 231 PetscCall(DMCreateMassMatrix(dmaProj, dmaProj, &mass)); 232 PetscCall(MatViewFromOptions(mass, NULL, "-mass_mat_view")); 233 PetscCall(KSPCreate(PETSC_COMM_WORLD, &ksp)); 234 PetscCall(KSPSetOperators(ksp, mass, mass)); 235 PetscCall(KSPSetFromOptions(ksp)); 236 // Compute rhs as M f, could also direclty project the analytic function but we might not have it 237 PetscCall(DMCreateMassMatrix(dmProj, dmaProj, &mass2)); 238 PetscCall(MatMult(mass2, u, rhs)); 239 PetscCall(MatDestroy(&mass2)); 240 PetscCall(KSPSolve(ksp, rhs, uproj)); 241 PetscCall(PetscObjectSetName((PetscObject) uproj, "L_2 Projection")); 242 PetscCall(VecViewFromOptions(uproj, NULL, "-proj_vec_view")); 243 PetscCall(DMComputeL2Diff(dmaProj, 0.0, funcs, NULL, uproj, &error)); 244 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Projected L2 Error: %g\n", (double) error)); 245 PetscCall(KSPDestroy(&ksp)); 246 PetscCall(MatDestroy(&mass)); 247 PetscCall(DMRestoreGlobalVector(dmProj, &u)); 248 PetscCall(DMRestoreGlobalVector(dmaProj, &ua)); 249 PetscCall(DMRestoreGlobalVector(dmaProj, &rhs)); 250 PetscCall(DMRestoreGlobalVector(dmaProj, &uproj)); 251 PetscCall(DMDestroy(&dmProj)); 252 PetscCall(DMDestroy(&dmaProj)); 253 PetscFunctionReturn(0); 254 } 255 256 int main (int argc, char * argv[]) { 257 DM dm; 258 AppCtx user; /* user-defined work context */ 259 MPI_Comm comm; 260 DM dma, odm; 261 Vec metric; 262 PetscInt r; 263 264 PetscCall(PetscInitialize(&argc, &argv, NULL, help)); 265 comm = PETSC_COMM_WORLD; 266 PetscCall(ProcessOptions(comm, &user)); 267 PetscCall(CreateMesh(comm, &dm)); 268 269 odm = dm; 270 PetscCall(DMPlexDistributeOverlap(odm, 1, NULL, &dm)); 271 if (!dm) {dm = odm;} 272 else PetscCall(DMDestroy(&odm)); 273 274 for (r = 0; r < user.Nr; ++r) { 275 DMLabel label; 276 277 PetscCall(ComputeMetric(dm, &user, &metric)); 278 PetscCall(DMGetLabel(dm, "marker", &label)); 279 PetscCall(DMAdaptMetric(dm, metric, label, NULL, &dma)); 280 PetscCall(VecDestroy(&metric)); 281 PetscCall(PetscObjectSetName((PetscObject) dma, "DMadapt")); 282 PetscCall(PetscObjectSetOptionsPrefix((PetscObject) dma, "adapt_")); 283 PetscCall(DMViewFromOptions(dma, NULL, "-dm_view")); 284 if (user.doL2) PetscCall(TestL2Projection(dm, dma, &user)); 285 PetscCall(DMDestroy(&dm)); 286 dm = dma; 287 } 288 PetscCall(PetscObjectSetOptionsPrefix((PetscObject) dm, "final_")); 289 PetscCall(DMViewFromOptions(dm, NULL, "-dm_view")); 290 PetscCall(DMDestroy(&dm)); 291 PetscCall(PetscFinalize()); 292 return 0; 293 } 294 295 /*TEST 296 297 build: 298 requires: pragmatic 299 300 testset: 301 args: -dm_plex_box_faces 4,4,4 -dm_adaptor pragmatic -met 2 -init_dm_view -adapt_dm_view -dm_adaptor pragmatic 302 303 test: 304 suffix: 0 305 args: -dm_plex_separate_marker 0 306 test: 307 suffix: 1 308 args: -dm_plex_separate_marker 1 309 test: 310 suffix: 2 311 args: -dm_plex_dim 3 312 test: 313 suffix: 3 314 args: -dm_plex_dim 3 315 316 # Pragmatic hangs for simple partitioner 317 testset: 318 requires: parmetis 319 args: -dm_plex_box_faces 2,2 -dm_adaptor pragmatic -petscpartitioner_type parmetis -met 2 -init_dm_view -adapt_dm_view -dm_adaptor pragmatic 320 321 test: 322 suffix: 4 323 nsize: 2 324 test: 325 suffix: 5 326 nsize: 4 327 328 test: 329 requires: parmetis 330 suffix: 6 331 nsize: 2 332 args: -dm_plex_dim 3 -dm_plex_box_faces 9,9,9 -dm_adaptor pragmatic -petscpartitioner_type parmetis \ 333 -met 0 -hmin 0.01 -hmax 0.03 -init_dm_view -adapt_dm_view -dm_adaptor pragmatic 334 test: 335 requires: parmetis 336 suffix: 7 337 nsize: 2 338 args: -dm_plex_box_faces 19,19 -dm_adaptor pragmatic -petscpartitioner_type parmetis \ 339 -met 2 -hmax 0.5 -hmin 0.001 -init_dm_view -adapt_dm_view -dm_adaptor pragmatic 340 test: 341 suffix: proj_0 342 args: -dm_plex_box_faces 2,2 -dm_plex_hash_location -dm_adaptor pragmatic -init_dm_view -adapt_dm_view -do_L2 \ 343 -petscspace_degree 1 -petscfe_default_quadrature_order 1 -pc_type lu -dm_adaptor pragmatic 344 test: 345 suffix: proj_1 346 args: -dm_plex_box_faces 4,4 -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 349 test: 350 suffix: sensor 351 args: -dm_plex_box_faces 9,9 -met 3 -dm_adaptor pragmatic -init_dm_view -adapt_dm_view \ 352 -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 \ 353 -dm_plex_metric_target_complexity 10000.0 -dm_adaptor pragmatic 354 355 TEST*/ 356