1 #include <petscdmadaptor.h> /*I "petscdmadaptor.h" I*/
2 #include <petscdmplex.h>
3 #include <petscdmforest.h>
4 #include <petscds.h>
5 #include <petscblaslapack.h>
6
7 #include <petsc/private/dmadaptorimpl.h>
8 #include <petsc/private/dmpleximpl.h>
9 #include <petsc/private/petscfeimpl.h>
10
11 static PetscErrorCode DMAdaptorSimpleErrorIndicator_Private(DMAdaptor, PetscInt, PetscInt, const PetscScalar *, const PetscScalar *, const PetscFVCellGeom *, PetscReal *, void *);
12
DMAdaptorTransferSolution_Exact_Private(DMAdaptor adaptor,DM dm,Vec u,DM adm,Vec au,void * ctx)13 static PetscErrorCode DMAdaptorTransferSolution_Exact_Private(DMAdaptor adaptor, DM dm, Vec u, DM adm, Vec au, void *ctx)
14 {
15 PetscErrorCode ierr;
16
17 PetscFunctionBeginUser;
18 ierr = DMProjectFunction(adm, 0.0, adaptor->exactSol, adaptor->exactCtx, INSERT_ALL_VALUES, au);CHKERRQ(ierr);
19 PetscFunctionReturn(0);
20 }
21
22 /*@
23 DMAdaptorCreate - Create a DMAdaptor object. Its purpose is to construct a adaptation DMLabel or metric Vec that can be used to modify the DM.
24
25 Collective
26
27 Input Parameter:
28 . comm - The communicator for the DMAdaptor object
29
30 Output Parameter:
31 . adaptor - The DMAdaptor object
32
33 Level: beginner
34
35 .seealso: DMAdaptorDestroy(), DMAdaptorAdapt()
36 @*/
DMAdaptorCreate(MPI_Comm comm,DMAdaptor * adaptor)37 PetscErrorCode DMAdaptorCreate(MPI_Comm comm, DMAdaptor *adaptor)
38 {
39 VecTaggerBox refineBox, coarsenBox;
40 PetscErrorCode ierr;
41
42 PetscFunctionBegin;
43 PetscValidPointer(adaptor, 2);
44 ierr = PetscSysInitializePackage();CHKERRQ(ierr);
45 ierr = PetscHeaderCreate(*adaptor, DM_CLASSID, "DMAdaptor", "DM Adaptor", "SNES", comm, DMAdaptorDestroy, DMAdaptorView);CHKERRQ(ierr);
46
47 (*adaptor)->monitor = PETSC_FALSE;
48 (*adaptor)->adaptCriterion = DM_ADAPTATION_NONE;
49 (*adaptor)->numSeq = 1;
50 (*adaptor)->Nadapt = -1;
51 (*adaptor)->refinementFactor = 2.0;
52 (*adaptor)->h_min = 1.;
53 (*adaptor)->h_max = 10000.;
54 (*adaptor)->ops->computeerrorindicator = DMAdaptorSimpleErrorIndicator_Private;
55 refineBox.min = refineBox.max = PETSC_MAX_REAL;
56 ierr = VecTaggerCreate(PetscObjectComm((PetscObject) *adaptor), &(*adaptor)->refineTag);CHKERRQ(ierr);
57 ierr = PetscObjectSetOptionsPrefix((PetscObject) (*adaptor)->refineTag, "refine_");CHKERRQ(ierr);
58 ierr = VecTaggerSetType((*adaptor)->refineTag, VECTAGGERABSOLUTE);CHKERRQ(ierr);
59 ierr = VecTaggerAbsoluteSetBox((*adaptor)->refineTag, &refineBox);CHKERRQ(ierr);
60 coarsenBox.min = coarsenBox.max = PETSC_MAX_REAL;
61 ierr = VecTaggerCreate(PetscObjectComm((PetscObject) *adaptor), &(*adaptor)->coarsenTag);CHKERRQ(ierr);
62 ierr = PetscObjectSetOptionsPrefix((PetscObject) (*adaptor)->coarsenTag, "coarsen_");CHKERRQ(ierr);
63 ierr = VecTaggerSetType((*adaptor)->coarsenTag, VECTAGGERABSOLUTE);CHKERRQ(ierr);
64 ierr = VecTaggerAbsoluteSetBox((*adaptor)->coarsenTag, &coarsenBox);CHKERRQ(ierr);
65 PetscFunctionReturn(0);
66 }
67
68 /*@
69 DMAdaptorDestroy - Destroys a DMAdaptor object
70
71 Collective on DMAdaptor
72
73 Input Parameter:
74 . adaptor - The DMAdaptor object
75
76 Level: beginner
77
78 .seealso: DMAdaptorCreate(), DMAdaptorAdapt()
79 @*/
DMAdaptorDestroy(DMAdaptor * adaptor)80 PetscErrorCode DMAdaptorDestroy(DMAdaptor *adaptor)
81 {
82 PetscErrorCode ierr;
83
84 PetscFunctionBegin;
85 if (!*adaptor) PetscFunctionReturn(0);
86 PetscValidHeaderSpecific((*adaptor), DM_CLASSID, 1);
87 if (--((PetscObject)(*adaptor))->refct > 0) {
88 *adaptor = NULL;
89 PetscFunctionReturn(0);
90 }
91 ierr = VecTaggerDestroy(&(*adaptor)->refineTag);CHKERRQ(ierr);
92 ierr = VecTaggerDestroy(&(*adaptor)->coarsenTag);CHKERRQ(ierr);
93 ierr = PetscFree2((*adaptor)->exactSol, (*adaptor)->exactCtx);CHKERRQ(ierr);
94 ierr = PetscHeaderDestroy(adaptor);CHKERRQ(ierr);
95 PetscFunctionReturn(0);
96 }
97
98 /*@
99 DMAdaptorSetFromOptions - Sets a DMAdaptor object from options
100
101 Collective on DMAdaptor
102
103 Input Parameters:
104 . adaptor - The DMAdaptor object
105
106 Options Database Keys:
107 + -adaptor_monitor <bool> : Monitor the adaptation process
108 . -adaptor_sequence_num <num> : Number of adaptations to generate an optimal grid
109 . -adaptor_target_num <num> : Set the target number of vertices N_adapt, -1 for automatic determination
110 . -adaptor_refinement_factor <r> : Set r such that N_adapt = r^dim N_orig
111 . -adaptor_metric_h_min <min> : Set the minimum eigenvalue of Hessian (sqr max edge length)
112 - -adaptor_metric_h_max <max> : Set the maximum eigenvalue of Hessian (sqr min edge length)
113
114 Level: beginner
115
116 .seealso: DMAdaptorCreate(), DMAdaptorAdapt()
117 @*/
DMAdaptorSetFromOptions(DMAdaptor adaptor)118 PetscErrorCode DMAdaptorSetFromOptions(DMAdaptor adaptor)
119 {
120 PetscErrorCode ierr;
121
122 PetscFunctionBegin;
123 ierr = PetscOptionsBegin(PetscObjectComm((PetscObject) adaptor), "", "DM Adaptor Options", "DMAdaptor");CHKERRQ(ierr);
124 ierr = PetscOptionsBool("-adaptor_monitor", "Monitor the adaptation process", "DMAdaptorMonitor", adaptor->monitor, &adaptor->monitor, NULL);CHKERRQ(ierr);
125 ierr = PetscOptionsInt("-adaptor_sequence_num", "Number of adaptations to generate an optimal grid", "DMAdaptorSetSequenceLength", adaptor->numSeq, &adaptor->numSeq, NULL);CHKERRQ(ierr);
126 ierr = PetscOptionsInt("-adaptor_target_num", "Set the target number of vertices N_adapt, -1 for automatic determination", "DMAdaptor", adaptor->Nadapt, &adaptor->Nadapt, NULL);CHKERRQ(ierr);
127 ierr = PetscOptionsReal("-adaptor_refinement_factor", "Set r such that N_adapt = r^dim N_orig", "DMAdaptor", adaptor->refinementFactor, &adaptor->refinementFactor, NULL);CHKERRQ(ierr);
128 ierr = PetscOptionsReal("-adaptor_metric_h_min", "Set the minimum eigenvalue of Hessian (sqr max edge length)", "DMAdaptor", adaptor->h_min, &adaptor->h_min, NULL);CHKERRQ(ierr);
129 ierr = PetscOptionsReal("-adaptor_metric_h_max", "Set the maximum eigenvalue of Hessian (sqr min edge length)", "DMAdaptor", adaptor->h_max, &adaptor->h_max, NULL);CHKERRQ(ierr);
130 ierr = PetscOptionsEnd();
131 ierr = VecTaggerSetFromOptions(adaptor->refineTag);CHKERRQ(ierr);
132 ierr = VecTaggerSetFromOptions(adaptor->coarsenTag);CHKERRQ(ierr);
133 PetscFunctionReturn(0);
134 }
135
136 /*@
137 DMAdaptorView - Views a DMAdaptor object
138
139 Collective on DMAdaptor
140
141 Input Parameters:
142 + adaptor - The DMAdaptor object
143 - viewer - The PetscViewer object
144
145 Level: beginner
146
147 .seealso: DMAdaptorCreate(), DMAdaptorAdapt()
148 @*/
DMAdaptorView(DMAdaptor adaptor,PetscViewer viewer)149 PetscErrorCode DMAdaptorView(DMAdaptor adaptor, PetscViewer viewer)
150 {
151 PetscErrorCode ierr;
152
153 PetscFunctionBegin;
154 ierr = PetscObjectPrintClassNamePrefixType((PetscObject) adaptor, viewer);CHKERRQ(ierr);
155 ierr = PetscViewerASCIIPrintf(viewer, "DM Adaptor\n");CHKERRQ(ierr);
156 ierr = PetscViewerASCIIPrintf(viewer, " sequence length: %D\n", adaptor->numSeq);CHKERRQ(ierr);
157 ierr = VecTaggerView(adaptor->refineTag, viewer);CHKERRQ(ierr);
158 ierr = VecTaggerView(adaptor->coarsenTag, viewer);CHKERRQ(ierr);
159 PetscFunctionReturn(0);
160 }
161
162 /*@
163 DMAdaptorGetSolver - Gets the solver used to produce discrete solutions
164
165 Not collective
166
167 Input Parameter:
168 . adaptor - The DMAdaptor object
169
170 Output Parameter:
171 . snes - The solver
172
173 Level: intermediate
174
175 .seealso: DMAdaptorSetSolver(), DMAdaptorCreate(), DMAdaptorAdapt()
176 @*/
DMAdaptorGetSolver(DMAdaptor adaptor,SNES * snes)177 PetscErrorCode DMAdaptorGetSolver(DMAdaptor adaptor, SNES *snes)
178 {
179 PetscFunctionBegin;
180 PetscValidHeaderSpecific(adaptor, DM_CLASSID, 1);
181 PetscValidPointer(snes, 2);
182 *snes = adaptor->snes;
183 PetscFunctionReturn(0);
184 }
185
186 /*@
187 DMAdaptorSetSolver - Sets the solver used to produce discrete solutions
188
189 Not collective
190
191 Input Parameters:
192 + adaptor - The DMAdaptor object
193 - snes - The solver
194
195 Level: intermediate
196
197 Note: The solver MUST have an attached DM/DS, so that we know the exact solution
198
199 .seealso: DMAdaptorGetSolver(), DMAdaptorCreate(), DMAdaptorAdapt()
200 @*/
DMAdaptorSetSolver(DMAdaptor adaptor,SNES snes)201 PetscErrorCode DMAdaptorSetSolver(DMAdaptor adaptor, SNES snes)
202 {
203 PetscErrorCode ierr;
204
205 PetscFunctionBegin;
206 PetscValidHeaderSpecific(adaptor, DM_CLASSID, 1);
207 PetscValidHeaderSpecific(snes, SNES_CLASSID, 2);
208 adaptor->snes = snes;
209 ierr = SNESGetDM(adaptor->snes, &adaptor->idm);CHKERRQ(ierr);
210 PetscFunctionReturn(0);
211 }
212
213 /*@
214 DMAdaptorGetSequenceLength - Gets the number of sequential adaptations
215
216 Not collective
217
218 Input Parameter:
219 . adaptor - The DMAdaptor object
220
221 Output Parameter:
222 . num - The number of adaptations
223
224 Level: intermediate
225
226 .seealso: DMAdaptorSetSequenceLength(), DMAdaptorCreate(), DMAdaptorAdapt()
227 @*/
DMAdaptorGetSequenceLength(DMAdaptor adaptor,PetscInt * num)228 PetscErrorCode DMAdaptorGetSequenceLength(DMAdaptor adaptor, PetscInt *num)
229 {
230 PetscFunctionBegin;
231 PetscValidHeaderSpecific(adaptor, DM_CLASSID, 1);
232 PetscValidPointer(num, 2);
233 *num = adaptor->numSeq;
234 PetscFunctionReturn(0);
235 }
236
237 /*@
238 DMAdaptorSetSequenceLength - Sets the number of sequential adaptations
239
240 Not collective
241
242 Input Parameters:
243 + adaptor - The DMAdaptor object
244 - num - The number of adaptations
245
246 Level: intermediate
247
248 .seealso: DMAdaptorGetSequenceLength(), DMAdaptorCreate(), DMAdaptorAdapt()
249 @*/
DMAdaptorSetSequenceLength(DMAdaptor adaptor,PetscInt num)250 PetscErrorCode DMAdaptorSetSequenceLength(DMAdaptor adaptor, PetscInt num)
251 {
252 PetscFunctionBegin;
253 PetscValidHeaderSpecific(adaptor, DM_CLASSID, 1);
254 adaptor->numSeq = num;
255 PetscFunctionReturn(0);
256 }
257
258 /*@
259 DMAdaptorSetUp - After the solver is specified, we create structures for controlling adaptivity
260
261 Collective on DMAdaptor
262
263 Input Parameters:
264 . adaptor - The DMAdaptor object
265
266 Level: beginner
267
268 .seealso: DMAdaptorCreate(), DMAdaptorAdapt()
269 @*/
DMAdaptorSetUp(DMAdaptor adaptor)270 PetscErrorCode DMAdaptorSetUp(DMAdaptor adaptor)
271 {
272 PetscDS prob;
273 PetscInt Nf, f;
274 PetscErrorCode ierr;
275
276 PetscFunctionBegin;
277 ierr = DMGetDS(adaptor->idm, &prob);CHKERRQ(ierr);
278 ierr = VecTaggerSetUp(adaptor->refineTag);CHKERRQ(ierr);
279 ierr = VecTaggerSetUp(adaptor->coarsenTag);CHKERRQ(ierr);
280 ierr = PetscDSGetNumFields(prob, &Nf);CHKERRQ(ierr);
281 ierr = PetscMalloc2(Nf, &adaptor->exactSol, Nf, &adaptor->exactCtx);CHKERRQ(ierr);
282 for (f = 0; f < Nf; ++f) {
283 ierr = PetscDSGetExactSolution(prob, f, &adaptor->exactSol[f], &adaptor->exactCtx[f]);CHKERRQ(ierr);
284 /* TODO Have a flag that forces projection rather than using the exact solution */
285 if (adaptor->exactSol[0]) {ierr = DMAdaptorSetTransferFunction(adaptor, DMAdaptorTransferSolution_Exact_Private);CHKERRQ(ierr);}
286 }
287 PetscFunctionReturn(0);
288 }
289
DMAdaptorGetTransferFunction(DMAdaptor adaptor,PetscErrorCode (** tfunc)(DMAdaptor,DM,Vec,DM,Vec,void *))290 PetscErrorCode DMAdaptorGetTransferFunction(DMAdaptor adaptor, PetscErrorCode (**tfunc)(DMAdaptor, DM, Vec, DM, Vec, void *))
291 {
292 PetscFunctionBegin;
293 *tfunc = adaptor->ops->transfersolution;
294 PetscFunctionReturn(0);
295 }
296
DMAdaptorSetTransferFunction(DMAdaptor adaptor,PetscErrorCode (* tfunc)(DMAdaptor,DM,Vec,DM,Vec,void *))297 PetscErrorCode DMAdaptorSetTransferFunction(DMAdaptor adaptor, PetscErrorCode (*tfunc)(DMAdaptor, DM, Vec, DM, Vec, void *))
298 {
299 PetscFunctionBegin;
300 adaptor->ops->transfersolution = tfunc;
301 PetscFunctionReturn(0);
302 }
303
DMAdaptorPreAdapt(DMAdaptor adaptor,Vec locX)304 PetscErrorCode DMAdaptorPreAdapt(DMAdaptor adaptor, Vec locX)
305 {
306 DM plex;
307 PetscDS prob;
308 PetscObject obj;
309 PetscClassId id;
310 PetscBool isForest;
311 PetscErrorCode ierr;
312
313 PetscFunctionBegin;
314 ierr = DMConvert(adaptor->idm, DMPLEX, &plex);CHKERRQ(ierr);
315 ierr = DMGetDS(adaptor->idm, &prob);CHKERRQ(ierr);
316 ierr = PetscDSGetDiscretization(prob, 0, &obj);CHKERRQ(ierr);
317 ierr = PetscObjectGetClassId(obj, &id);CHKERRQ(ierr);
318 ierr = DMIsForest(adaptor->idm, &isForest);CHKERRQ(ierr);
319 if (adaptor->adaptCriterion == DM_ADAPTATION_NONE) {
320 if (isForest) {adaptor->adaptCriterion = DM_ADAPTATION_LABEL;}
321 #if defined(PETSC_HAVE_PRAGMATIC)
322 else {adaptor->adaptCriterion = DM_ADAPTATION_METRIC;}
323 #else
324 else {adaptor->adaptCriterion = DM_ADAPTATION_REFINE;}
325 #endif
326 }
327 if (id == PETSCFV_CLASSID) {adaptor->femType = PETSC_FALSE;}
328 else {adaptor->femType = PETSC_TRUE;}
329 if (adaptor->femType) {
330 /* Compute local solution bc */
331 ierr = DMPlexInsertBoundaryValues(plex, PETSC_TRUE, locX, 0.0, adaptor->faceGeom, adaptor->cellGeom, NULL);CHKERRQ(ierr);
332 } else {
333 PetscFV fvm = (PetscFV) obj;
334 PetscLimiter noneLimiter;
335 Vec grad;
336
337 ierr = PetscFVGetComputeGradients(fvm, &adaptor->computeGradient);CHKERRQ(ierr);
338 ierr = PetscFVSetComputeGradients(fvm, PETSC_TRUE);CHKERRQ(ierr);
339 /* Use no limiting when reconstructing gradients for adaptivity */
340 ierr = PetscFVGetLimiter(fvm, &adaptor->limiter);CHKERRQ(ierr);
341 ierr = PetscObjectReference((PetscObject) adaptor->limiter);CHKERRQ(ierr);
342 ierr = PetscLimiterCreate(PetscObjectComm((PetscObject) fvm), &noneLimiter);CHKERRQ(ierr);
343 ierr = PetscLimiterSetType(noneLimiter, PETSCLIMITERNONE);CHKERRQ(ierr);
344 ierr = PetscFVSetLimiter(fvm, noneLimiter);CHKERRQ(ierr);
345 /* Get FVM data */
346 ierr = DMPlexGetDataFVM(plex, fvm, &adaptor->cellGeom, &adaptor->faceGeom, &adaptor->gradDM);CHKERRQ(ierr);
347 ierr = VecGetDM(adaptor->cellGeom, &adaptor->cellDM);CHKERRQ(ierr);
348 ierr = VecGetArrayRead(adaptor->cellGeom, &adaptor->cellGeomArray);CHKERRQ(ierr);
349 /* Compute local solution bc */
350 ierr = DMPlexInsertBoundaryValues(plex, PETSC_TRUE, locX, 0.0, adaptor->faceGeom, adaptor->cellGeom, NULL);CHKERRQ(ierr);
351 /* Compute gradients */
352 ierr = DMCreateGlobalVector(adaptor->gradDM, &grad);CHKERRQ(ierr);
353 ierr = DMPlexReconstructGradientsFVM(plex, locX, grad);CHKERRQ(ierr);
354 ierr = DMGetLocalVector(adaptor->gradDM, &adaptor->cellGrad);CHKERRQ(ierr);
355 ierr = DMGlobalToLocalBegin(adaptor->gradDM, grad, INSERT_VALUES, adaptor->cellGrad);CHKERRQ(ierr);
356 ierr = DMGlobalToLocalEnd(adaptor->gradDM, grad, INSERT_VALUES, adaptor->cellGrad);CHKERRQ(ierr);
357 ierr = VecDestroy(&grad);CHKERRQ(ierr);
358 ierr = VecGetArrayRead(adaptor->cellGrad, &adaptor->cellGradArray);CHKERRQ(ierr);
359 }
360 ierr = DMDestroy(&plex);CHKERRQ(ierr);
361 PetscFunctionReturn(0);
362 }
363
DMAdaptorTransferSolution(DMAdaptor adaptor,DM dm,Vec x,DM adm,Vec ax)364 PetscErrorCode DMAdaptorTransferSolution(DMAdaptor adaptor, DM dm, Vec x, DM adm, Vec ax)
365 {
366 PetscReal time = 0.0;
367 Mat interp;
368 void *ctx;
369 PetscErrorCode ierr;
370
371 PetscFunctionBegin;
372 ierr = DMGetApplicationContext(dm, &ctx);CHKERRQ(ierr);
373 if (adaptor->ops->transfersolution) {
374 ierr = (*adaptor->ops->transfersolution)(adaptor, dm, x, adm, ax, ctx);CHKERRQ(ierr);
375 } else {
376 switch (adaptor->adaptCriterion) {
377 case DM_ADAPTATION_LABEL:
378 ierr = DMForestTransferVec(dm, x, adm, ax, PETSC_TRUE, time);CHKERRQ(ierr);
379 break;
380 case DM_ADAPTATION_REFINE:
381 case DM_ADAPTATION_METRIC:
382 ierr = DMCreateInterpolation(dm, adm, &interp, NULL);CHKERRQ(ierr);
383 ierr = MatInterpolate(interp, x, ax);CHKERRQ(ierr);
384 ierr = DMInterpolate(dm, interp, adm);CHKERRQ(ierr);
385 ierr = MatDestroy(&interp);CHKERRQ(ierr);
386 break;
387 default: SETERRQ1(PetscObjectComm((PetscObject) adaptor), PETSC_ERR_SUP, "No built-in projection for this adaptation criterion: %D", adaptor->adaptCriterion);
388 }
389 }
390 PetscFunctionReturn(0);
391 }
392
DMAdaptorPostAdapt(DMAdaptor adaptor)393 PetscErrorCode DMAdaptorPostAdapt(DMAdaptor adaptor)
394 {
395 PetscDS prob;
396 PetscObject obj;
397 PetscClassId id;
398 PetscErrorCode ierr;
399
400 PetscFunctionBegin;
401 ierr = DMGetDS(adaptor->idm, &prob);CHKERRQ(ierr);
402 ierr = PetscDSGetDiscretization(prob, 0, &obj);CHKERRQ(ierr);
403 ierr = PetscObjectGetClassId(obj, &id);CHKERRQ(ierr);
404 if (id == PETSCFV_CLASSID) {
405 PetscFV fvm = (PetscFV) obj;
406
407 ierr = PetscFVSetComputeGradients(fvm, adaptor->computeGradient);CHKERRQ(ierr);
408 /* Restore original limiter */
409 ierr = PetscFVSetLimiter(fvm, adaptor->limiter);CHKERRQ(ierr);
410
411 ierr = VecRestoreArrayRead(adaptor->cellGeom, &adaptor->cellGeomArray);CHKERRQ(ierr);
412 ierr = VecRestoreArrayRead(adaptor->cellGrad, &adaptor->cellGradArray);CHKERRQ(ierr);
413 ierr = DMRestoreLocalVector(adaptor->gradDM, &adaptor->cellGrad);CHKERRQ(ierr);
414 }
415 PetscFunctionReturn(0);
416 }
417
DMAdaptorModifyHessian_Private(PetscInt dim,PetscReal h_min,PetscReal h_max,PetscScalar Hp[])418 static PetscErrorCode DMAdaptorModifyHessian_Private(PetscInt dim, PetscReal h_min, PetscReal h_max, PetscScalar Hp[])
419 {
420 PetscScalar *Hpos;
421 PetscReal *eigs;
422 PetscInt i, j, k;
423 PetscErrorCode ierr;
424
425 PetscFunctionBegin;
426 ierr = PetscMalloc2(dim*dim, &Hpos, dim, &eigs);CHKERRQ(ierr);
427 #if 0
428 ierr = PetscPrintf(PETSC_COMM_SELF, "H = [");CHKERRQ(ierr);
429 for (i = 0; i < dim; ++i) {
430 if (i > 0) {ierr = PetscPrintf(PETSC_COMM_SELF, " ");CHKERRQ(ierr);}
431 for (j = 0; j < dim; ++j) {
432 if (j > 0) {ierr = PetscPrintf(PETSC_COMM_SELF, ", ");CHKERRQ(ierr);}
433 ierr = PetscPrintf(PETSC_COMM_SELF, "%g", Hp[i*dim+j]);CHKERRQ(ierr);
434 }
435 ierr = PetscPrintf(PETSC_COMM_SELF, "\n");CHKERRQ(ierr);
436 }
437 ierr = PetscPrintf(PETSC_COMM_SELF, "]\n");CHKERRQ(ierr);
438 #endif
439 /* Symmetrize */
440 for (i = 0; i < dim; ++i) {
441 Hpos[i*dim+i] = Hp[i*dim+i];
442 for (j = i+1; j < dim; ++j) {
443 Hpos[i*dim+j] = 0.5*(Hp[i*dim+j] + Hp[j*dim+i]);
444 Hpos[j*dim+i] = Hpos[i*dim+j];
445 }
446 }
447 #if 0
448 ierr = PetscPrintf(PETSC_COMM_SELF, "Hs = [");CHKERRQ(ierr);
449 for (i = 0; i < dim; ++i) {
450 if (i > 0) {ierr = PetscPrintf(PETSC_COMM_SELF, " ");CHKERRQ(ierr);}
451 for (j = 0; j < dim; ++j) {
452 if (j > 0) {ierr = PetscPrintf(PETSC_COMM_SELF, ", ");CHKERRQ(ierr);}
453 ierr = PetscPrintf(PETSC_COMM_SELF, "%g", Hpos[i*dim+j]);CHKERRQ(ierr);
454 }
455 ierr = PetscPrintf(PETSC_COMM_SELF, "\n");CHKERRQ(ierr);
456 }
457 ierr = PetscPrintf(PETSC_COMM_SELF, "]\n");CHKERRQ(ierr);
458 #endif
459 /* Compute eigendecomposition */
460 {
461 PetscScalar *work;
462 PetscBLASInt lwork;
463
464 lwork = 5*dim;
465 ierr = PetscMalloc1(5*dim, &work);CHKERRQ(ierr);
466 {
467 PetscBLASInt lierr;
468 PetscBLASInt nb;
469
470 ierr = PetscBLASIntCast(dim, &nb);CHKERRQ(ierr);
471 ierr = PetscFPTrapPush(PETSC_FP_TRAP_OFF);CHKERRQ(ierr);
472 #if defined(PETSC_USE_COMPLEX)
473 {
474 PetscReal *rwork;
475 ierr = PetscMalloc1(3*dim, &rwork);CHKERRQ(ierr);
476 PetscStackCallBLAS("LAPACKsyev",LAPACKsyev_("V","U",&nb,Hpos,&nb,eigs,work,&lwork,rwork,&lierr));
477 ierr = PetscFree(rwork);CHKERRQ(ierr);
478 }
479 #else
480 PetscStackCallBLAS("LAPACKsyev",LAPACKsyev_("V","U",&nb,Hpos,&nb,eigs,work,&lwork,&lierr));
481 #endif
482 if (lierr) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_LIB, "Error in LAPACK routine %d", (int) lierr);
483 ierr = PetscFPTrapPop();CHKERRQ(ierr);
484 }
485 ierr = PetscFree(work);CHKERRQ(ierr);
486 }
487 #if 0
488 ierr = PetscPrintf(PETSC_COMM_SELF, "L = [");CHKERRQ(ierr);
489 for (i = 0; i < dim; ++i) {
490 if (i > 0) {ierr = PetscPrintf(PETSC_COMM_SELF, ", ");CHKERRQ(ierr);}
491 ierr = PetscPrintf(PETSC_COMM_SELF, "%g", eigs[i]);CHKERRQ(ierr);
492 }
493 ierr = PetscPrintf(PETSC_COMM_SELF, "]\n");CHKERRQ(ierr);
494 ierr = PetscPrintf(PETSC_COMM_SELF, "Q = [");CHKERRQ(ierr);
495 for (i = 0; i < dim; ++i) {
496 if (i > 0) {ierr = PetscPrintf(PETSC_COMM_SELF, " ");CHKERRQ(ierr);}
497 for (j = 0; j < dim; ++j) {
498 if (j > 0) {ierr = PetscPrintf(PETSC_COMM_SELF, ", ");CHKERRQ(ierr);}
499 ierr = PetscPrintf(PETSC_COMM_SELF, "%g", Hpos[i*dim+j]);CHKERRQ(ierr);
500 }
501 ierr = PetscPrintf(PETSC_COMM_SELF, "\n");CHKERRQ(ierr);
502 }
503 ierr = PetscPrintf(PETSC_COMM_SELF, "]\n");CHKERRQ(ierr);
504 #endif
505 /* Reflect to positive orthant, enforce maximum and minimum size, \lambda \propto 1/h^2
506 TODO get domain bounding box */
507 for (i = 0; i < dim; ++i) eigs[i] = PetscMin(h_max, PetscMax(h_min, PetscAbsReal(eigs[i])));
508 /* Reconstruct Hessian */
509 for (i = 0; i < dim; ++i) {
510 for (j = 0; j < dim; ++j) {
511 Hp[i*dim+j] = 0.0;
512 for (k = 0; k < dim; ++k) {
513 Hp[i*dim+j] += Hpos[k*dim+i] * eigs[k] * Hpos[k*dim+j];
514 }
515 }
516 }
517 ierr = PetscFree2(Hpos, eigs);CHKERRQ(ierr);
518 #if 0
519 ierr = PetscPrintf(PETSC_COMM_SELF, "H+ = [");CHKERRQ(ierr);
520 for (i = 0; i < dim; ++i) {
521 if (i > 0) {ierr = PetscPrintf(PETSC_COMM_SELF, " ");CHKERRQ(ierr);}
522 for (j = 0; j < dim; ++j) {
523 if (j > 0) {ierr = PetscPrintf(PETSC_COMM_SELF, ", ");CHKERRQ(ierr);}
524 ierr = PetscPrintf(PETSC_COMM_SELF, "%g", Hp[i*dim+j]);CHKERRQ(ierr);
525 }
526 ierr = PetscPrintf(PETSC_COMM_SELF, "\n");CHKERRQ(ierr);
527 }
528 ierr = PetscPrintf(PETSC_COMM_SELF, "]\n");CHKERRQ(ierr);
529 #endif
530 PetscFunctionReturn(0);
531 }
532
detHFunc(PetscInt dim,PetscInt Nf,PetscInt NfAux,const PetscInt uOff[],const PetscInt uOff_x[],const PetscScalar u[],const PetscScalar u_t[],const PetscScalar u_x[],const PetscInt aOff[],const PetscInt aOff_x[],const PetscScalar a[],const PetscScalar a_t[],const PetscScalar a_x[],PetscReal t,const PetscReal x[],PetscInt numConstants,const PetscScalar constants[],PetscScalar f0[])533 static void detHFunc(PetscInt dim, PetscInt Nf, PetscInt NfAux,
534 const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
535 const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
536 PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[])
537 {
538 const PetscInt p = 1;
539 PetscReal detH = 0.0;
540
541 if (dim == 2) DMPlex_Det2D_Scalar_Internal(&detH, u);
542 else if (dim == 3) DMPlex_Det3D_Scalar_Internal(&detH, u);
543 f0[0] = PetscPowReal(detH, p/(2.*p + dim));
544 }
545
546 /*
547 DMAdaptorSimpleErrorIndicator - Just use the integrated gradient as an error indicator
548
549 Input Parameters:
550 + adaptor - The DMAdaptor object
551 . dim - The topological dimension
552 . cell - The cell
553 . field - The field integrated over the cell
554 . gradient - The gradient integrated over the cell
555 . cg - A PetscFVCellGeom struct
556 - ctx - A user context
557
558 Output Parameter:
559 . errInd - The error indicator
560
561 .seealso: DMAdaptorComputeErrorIndicator()
562 */
DMAdaptorSimpleErrorIndicator_Private(DMAdaptor adaptor,PetscInt dim,PetscInt Nc,const PetscScalar * field,const PetscScalar * gradient,const PetscFVCellGeom * cg,PetscReal * errInd,void * ctx)563 static PetscErrorCode DMAdaptorSimpleErrorIndicator_Private(DMAdaptor adaptor, PetscInt dim, PetscInt Nc, const PetscScalar *field, const PetscScalar *gradient, const PetscFVCellGeom *cg, PetscReal *errInd, void *ctx)
564 {
565 PetscReal err = 0.;
566 PetscInt c, d;
567
568 PetscFunctionBeginHot;
569 for (c = 0; c < Nc; c++) {
570 for (d = 0; d < dim; ++d) {
571 err += PetscSqr(PetscRealPart(gradient[c*dim+d]));
572 }
573 }
574 *errInd = cg->volume * err;
575 PetscFunctionReturn(0);
576 }
577
DMAdaptorComputeErrorIndicator_Private(DMAdaptor adaptor,DM plex,PetscInt cell,Vec locX,PetscReal * errInd)578 static PetscErrorCode DMAdaptorComputeErrorIndicator_Private(DMAdaptor adaptor, DM plex, PetscInt cell, Vec locX, PetscReal *errInd)
579 {
580 PetscDS prob;
581 PetscObject obj;
582 PetscClassId id;
583 void *ctx;
584 PetscQuadrature quad;
585 PetscInt dim, d, cdim, Nc;
586 PetscErrorCode ierr;
587
588 PetscFunctionBegin;
589 *errInd = 0.;
590 ierr = DMGetDimension(plex, &dim);CHKERRQ(ierr);
591 ierr = DMGetCoordinateDim(plex, &cdim);CHKERRQ(ierr);
592 ierr = DMGetApplicationContext(plex, &ctx);CHKERRQ(ierr);
593 ierr = DMGetDS(plex, &prob);CHKERRQ(ierr);
594 ierr = PetscDSGetDiscretization(prob, 0, &obj);CHKERRQ(ierr);
595 ierr = PetscObjectGetClassId(obj, &id);CHKERRQ(ierr);
596 if (id == PETSCFV_CLASSID) {
597 const PetscScalar *pointSols;
598 const PetscScalar *pointSol;
599 const PetscScalar *pointGrad;
600 PetscFVCellGeom *cg;
601
602 ierr = PetscFVGetNumComponents((PetscFV) obj, &Nc);CHKERRQ(ierr);
603 ierr = VecGetArrayRead(locX, &pointSols);CHKERRQ(ierr);
604 ierr = DMPlexPointLocalRead(plex, cell, pointSols, (void *) &pointSol);CHKERRQ(ierr);
605 ierr = DMPlexPointLocalRead(adaptor->gradDM, cell, adaptor->cellGradArray, (void *) &pointGrad);CHKERRQ(ierr);
606 ierr = DMPlexPointLocalRead(adaptor->cellDM, cell, adaptor->cellGeomArray, &cg);CHKERRQ(ierr);
607 ierr = (*adaptor->ops->computeerrorindicator)(adaptor, dim, Nc, pointSol, pointGrad, cg, errInd, ctx);CHKERRQ(ierr);
608 ierr = VecRestoreArrayRead(locX, &pointSols);CHKERRQ(ierr);
609 } else {
610 PetscScalar *x = NULL, *field, *gradient, *interpolant, *interpolantGrad;
611 PetscFVCellGeom cg;
612 PetscFEGeom fegeom;
613 const PetscReal *quadWeights;
614 PetscReal *coords;
615 PetscInt Nb, fc, Nq, qNc, Nf, f, fieldOffset;
616
617 fegeom.dim = dim;
618 fegeom.dimEmbed = cdim;
619 ierr = PetscDSGetNumFields(prob, &Nf);CHKERRQ(ierr);
620 ierr = PetscFEGetQuadrature((PetscFE) obj, &quad);CHKERRQ(ierr);
621 ierr = DMPlexVecGetClosure(plex, NULL, locX, cell, NULL, &x);CHKERRQ(ierr);
622 ierr = PetscFEGetDimension((PetscFE) obj, &Nb);CHKERRQ(ierr);
623 ierr = PetscFEGetNumComponents((PetscFE) obj, &Nc);CHKERRQ(ierr);
624 ierr = PetscQuadratureGetData(quad, NULL, &qNc, &Nq, NULL, &quadWeights);CHKERRQ(ierr);
625 ierr = PetscMalloc6(Nc,&field,cdim*Nc,&gradient,cdim*Nq,&coords,Nq,&fegeom.detJ,cdim*cdim*Nq,&fegeom.J,cdim*cdim*Nq,&fegeom.invJ);CHKERRQ(ierr);
626 ierr = PetscMalloc2(Nc, &interpolant, cdim*Nc, &interpolantGrad);CHKERRQ(ierr);
627 ierr = DMPlexComputeCellGeometryFEM(plex, cell, quad, coords, fegeom.J, fegeom.invJ, fegeom.detJ);CHKERRQ(ierr);
628 ierr = DMPlexComputeCellGeometryFVM(plex, cell, &cg.volume, NULL, NULL);CHKERRQ(ierr);
629 ierr = PetscArrayzero(gradient, cdim*Nc);CHKERRQ(ierr);
630 for (f = 0, fieldOffset = 0; f < Nf; ++f) {
631 PetscInt qc = 0, q;
632
633 ierr = PetscDSGetDiscretization(prob, f, &obj);CHKERRQ(ierr);
634 ierr = PetscArrayzero(interpolant,Nc);CHKERRQ(ierr);
635 ierr = PetscArrayzero(interpolantGrad, cdim*Nc);CHKERRQ(ierr);
636 for (q = 0; q < Nq; ++q) {
637 ierr = PetscFEInterpolateFieldAndGradient_Static((PetscFE) obj, x, &fegeom, q, interpolant, interpolantGrad);CHKERRQ(ierr);
638 for (fc = 0; fc < Nc; ++fc) {
639 const PetscReal wt = quadWeights[q*qNc+qc+fc];
640
641 field[fc] += interpolant[fc]*wt*fegeom.detJ[q];
642 for (d = 0; d < cdim; ++d) gradient[fc*cdim+d] += interpolantGrad[fc*dim+d]*wt*fegeom.detJ[q];
643 }
644 }
645 fieldOffset += Nb;
646 qc += Nc;
647 }
648 ierr = PetscFree2(interpolant, interpolantGrad);CHKERRQ(ierr);
649 ierr = DMPlexVecRestoreClosure(plex, NULL, locX, cell, NULL, &x);CHKERRQ(ierr);
650 for (fc = 0; fc < Nc; ++fc) {
651 field[fc] /= cg.volume;
652 for (d = 0; d < cdim; ++d) gradient[fc*cdim+d] /= cg.volume;
653 }
654 ierr = (*adaptor->ops->computeerrorindicator)(adaptor, dim, Nc, field, gradient, &cg, errInd, ctx);CHKERRQ(ierr);
655 ierr = PetscFree6(field,gradient,coords,fegeom.detJ,fegeom.J,fegeom.invJ);CHKERRQ(ierr);
656 }
657 PetscFunctionReturn(0);
658 }
659
DMAdaptorAdapt_Sequence_Private(DMAdaptor adaptor,Vec inx,PetscBool doSolve,DM * adm,Vec * ax)660 static PetscErrorCode DMAdaptorAdapt_Sequence_Private(DMAdaptor adaptor, Vec inx, PetscBool doSolve, DM *adm, Vec *ax)
661 {
662 PetscDS prob;
663 void *ctx;
664 MPI_Comm comm;
665 PetscInt numAdapt = adaptor->numSeq, adaptIter;
666 PetscInt dim, coordDim, numFields, cStart, cEnd, c;
667 PetscErrorCode ierr;
668
669 PetscFunctionBegin;
670 ierr = DMViewFromOptions(adaptor->idm, NULL, "-dm_adapt_pre_view");CHKERRQ(ierr);
671 ierr = VecViewFromOptions(inx, NULL, "-sol_adapt_pre_view");CHKERRQ(ierr);
672 ierr = PetscObjectGetComm((PetscObject) adaptor, &comm);CHKERRQ(ierr);
673 ierr = DMGetDimension(adaptor->idm, &dim);CHKERRQ(ierr);
674 ierr = DMGetCoordinateDim(adaptor->idm, &coordDim);CHKERRQ(ierr);
675 ierr = DMGetApplicationContext(adaptor->idm, &ctx);CHKERRQ(ierr);
676 ierr = DMGetDS(adaptor->idm, &prob);CHKERRQ(ierr);
677 ierr = PetscDSGetNumFields(prob, &numFields);CHKERRQ(ierr);
678
679 /* Adapt until nothing changes */
680 /* Adapt for a specified number of iterates */
681 for (adaptIter = 0; adaptIter < numAdapt-1; ++adaptIter) {ierr = PetscViewerASCIIPushTab(PETSC_VIEWER_STDOUT_(comm));CHKERRQ(ierr);}
682 for (adaptIter = 0; adaptIter < numAdapt; ++adaptIter) {
683 PetscBool adapted = PETSC_FALSE;
684 DM dm = adaptIter ? *adm : adaptor->idm, odm;
685 Vec x = adaptIter ? *ax : inx, locX, ox;
686
687 ierr = DMGetLocalVector(dm, &locX);CHKERRQ(ierr);
688 ierr = DMGlobalToLocalBegin(dm, adaptIter ? *ax : x, INSERT_VALUES, locX);CHKERRQ(ierr);
689 ierr = DMGlobalToLocalEnd(dm, adaptIter ? *ax : x, INSERT_VALUES, locX);CHKERRQ(ierr);
690 ierr = DMAdaptorPreAdapt(adaptor, locX);CHKERRQ(ierr);
691 if (doSolve) {
692 SNES snes;
693
694 ierr = DMAdaptorGetSolver(adaptor, &snes);CHKERRQ(ierr);
695 ierr = SNESSolve(snes, NULL, adaptIter ? *ax : x);CHKERRQ(ierr);
696 }
697 /* ierr = DMAdaptorMonitor(adaptor);CHKERRQ(ierr);
698 Print iterate, memory used, DM, solution */
699 switch (adaptor->adaptCriterion) {
700 case DM_ADAPTATION_REFINE:
701 ierr = DMRefine(dm, comm, &odm);CHKERRQ(ierr);
702 if (!odm) SETERRQ(comm, PETSC_ERR_ARG_INCOMP, "DMRefine() did not perform any refinement, cannot continue grid sequencing");
703 adapted = PETSC_TRUE;
704 break;
705 case DM_ADAPTATION_LABEL:
706 {
707 /* Adapt DM
708 Create local solution
709 Reconstruct gradients (FVM) or solve adjoint equation (FEM)
710 Produce cellwise error indicator */
711 DM plex;
712 DMLabel adaptLabel;
713 IS refineIS, coarsenIS;
714 Vec errVec;
715 PetscScalar *errArray;
716 const PetscScalar *pointSols;
717 PetscReal minMaxInd[2] = {PETSC_MAX_REAL, PETSC_MIN_REAL}, minMaxIndGlobal[2];
718 PetscInt nRefine, nCoarsen;
719
720 ierr = DMConvert(dm, DMPLEX, &plex);CHKERRQ(ierr);
721 ierr = DMLabelCreate(PETSC_COMM_SELF, "adapt", &adaptLabel);CHKERRQ(ierr);
722 ierr = DMPlexGetSimplexOrBoxCells(plex, 0, &cStart, &cEnd);CHKERRQ(ierr);
723
724 ierr = VecCreateMPI(PetscObjectComm((PetscObject) adaptor), cEnd-cStart, PETSC_DETERMINE, &errVec);CHKERRQ(ierr);
725 ierr = VecSetUp(errVec);CHKERRQ(ierr);
726 ierr = VecGetArray(errVec, &errArray);CHKERRQ(ierr);
727 ierr = VecGetArrayRead(locX, &pointSols);CHKERRQ(ierr);
728 for (c = cStart; c < cEnd; ++c) {
729 PetscReal errInd;
730
731 ierr = DMAdaptorComputeErrorIndicator_Private(adaptor, plex, c, locX, &errInd);CHKERRQ(ierr);
732 errArray[c-cStart] = errInd;
733 minMaxInd[0] = PetscMin(minMaxInd[0], errInd);
734 minMaxInd[1] = PetscMax(minMaxInd[1], errInd);
735 }
736 ierr = VecRestoreArrayRead(locX, &pointSols);CHKERRQ(ierr);
737 ierr = VecRestoreArray(errVec, &errArray);CHKERRQ(ierr);
738 ierr = PetscGlobalMinMaxReal(PetscObjectComm((PetscObject) adaptor), minMaxInd, minMaxIndGlobal);CHKERRQ(ierr);
739 ierr = PetscInfo2(adaptor, "DMAdaptor: error indicator range (%E, %E)\n", minMaxIndGlobal[0], minMaxIndGlobal[1]);CHKERRQ(ierr);
740 /* Compute IS from VecTagger */
741 ierr = VecTaggerComputeIS(adaptor->refineTag, errVec, &refineIS);CHKERRQ(ierr);
742 ierr = VecTaggerComputeIS(adaptor->coarsenTag, errVec, &coarsenIS);CHKERRQ(ierr);
743 ierr = ISGetSize(refineIS, &nRefine);CHKERRQ(ierr);
744 ierr = ISGetSize(coarsenIS, &nCoarsen);CHKERRQ(ierr);
745 ierr = PetscInfo2(adaptor, "DMAdaptor: numRefine %D, numCoarsen %D\n", nRefine, nCoarsen);CHKERRQ(ierr);
746 if (nRefine) {ierr = DMLabelSetStratumIS(adaptLabel, DM_ADAPT_REFINE, refineIS);CHKERRQ(ierr);}
747 if (nCoarsen) {ierr = DMLabelSetStratumIS(adaptLabel, DM_ADAPT_COARSEN, coarsenIS);CHKERRQ(ierr);}
748 ierr = ISDestroy(&coarsenIS);CHKERRQ(ierr);
749 ierr = ISDestroy(&refineIS);CHKERRQ(ierr);
750 ierr = VecDestroy(&errVec);CHKERRQ(ierr);
751 /* Adapt DM from label */
752 if (nRefine || nCoarsen) {
753 ierr = DMAdaptLabel(dm, adaptLabel, &odm);CHKERRQ(ierr);
754 adapted = PETSC_TRUE;
755 }
756 ierr = DMLabelDestroy(&adaptLabel);CHKERRQ(ierr);
757 ierr = DMDestroy(&plex);CHKERRQ(ierr);
758 }
759 break;
760 case DM_ADAPTATION_METRIC:
761 {
762 DM dmGrad, dmHess, dmMetric;
763 PetscDS probGrad, probHess;
764 Vec xGrad, xHess, metric;
765 PetscSection sec, msec;
766 PetscScalar *H, *M, integral;
767 PetscReal N;
768 DMLabel bdLabel;
769 PetscInt Nd = coordDim*coordDim, f, vStart, vEnd, v;
770
771 /* Compute vertexwise gradients from cellwise gradients */
772 ierr = DMClone(dm, &dmGrad);CHKERRQ(ierr);
773 ierr = DMClone(dm, &dmHess);CHKERRQ(ierr);
774 ierr = DMGetDS(dmGrad, &probGrad);CHKERRQ(ierr);
775 ierr = DMGetDS(dmHess, &probHess);CHKERRQ(ierr);
776 for (f = 0; f < numFields; ++f) {
777 PetscFE fe, feGrad, feHess;
778 PetscDualSpace Q;
779 DM K;
780 PetscQuadrature q;
781 PetscInt Nc, qorder;
782 const char *prefix;
783
784 ierr = PetscDSGetDiscretization(prob, f, (PetscObject *) &fe);CHKERRQ(ierr);
785 ierr = PetscFEGetNumComponents(fe, &Nc);CHKERRQ(ierr);
786 ierr = PetscFEGetDualSpace(fe, &Q);CHKERRQ(ierr);
787 ierr = PetscDualSpaceGetDM(Q, &K);CHKERRQ(ierr);
788 ierr = DMPlexGetDepthStratum(K, 0, &vStart, &vEnd);CHKERRQ(ierr);
789 ierr = PetscFEGetQuadrature(fe, &q);CHKERRQ(ierr);
790 ierr = PetscQuadratureGetOrder(q, &qorder);CHKERRQ(ierr);
791 ierr = PetscObjectGetOptionsPrefix((PetscObject) fe, &prefix);CHKERRQ(ierr);
792 ierr = PetscFECreateDefault(PetscObjectComm((PetscObject) dmGrad), dim, Nc*coordDim, (vEnd-vStart) == dim+1 ? PETSC_TRUE : PETSC_FALSE, prefix, qorder, &feGrad);CHKERRQ(ierr);
793 ierr = PetscDSSetDiscretization(probGrad, f, (PetscObject) feGrad);CHKERRQ(ierr);
794 ierr = PetscFECreateDefault(PetscObjectComm((PetscObject) dmHess), dim, Nc*Nd, (vEnd-vStart) == dim+1 ? PETSC_TRUE : PETSC_FALSE, prefix, qorder, &feHess);CHKERRQ(ierr);
795 ierr = PetscDSSetDiscretization(probHess, f, (PetscObject) feHess);CHKERRQ(ierr);
796 ierr = PetscFEDestroy(&feGrad);CHKERRQ(ierr);
797 ierr = PetscFEDestroy(&feHess);CHKERRQ(ierr);
798 }
799 ierr = DMGetGlobalVector(dmGrad, &xGrad);CHKERRQ(ierr);
800 ierr = VecViewFromOptions(x, NULL, "-sol_adapt_loc_pre_view");CHKERRQ(ierr);
801 ierr = DMPlexComputeGradientClementInterpolant(dm, locX, xGrad);CHKERRQ(ierr);
802 ierr = VecViewFromOptions(xGrad, NULL, "-adapt_gradient_view");CHKERRQ(ierr);
803 /* Compute vertexwise Hessians from cellwise Hessians */
804 ierr = DMGetGlobalVector(dmHess, &xHess);CHKERRQ(ierr);
805 ierr = DMPlexComputeGradientClementInterpolant(dmGrad, xGrad, xHess);CHKERRQ(ierr);
806 ierr = VecViewFromOptions(xHess, NULL, "-adapt_hessian_view");CHKERRQ(ierr);
807 /* Compute metric */
808 ierr = DMClone(dm, &dmMetric);CHKERRQ(ierr);
809 ierr = DMGetLocalSection(dm, &sec);CHKERRQ(ierr);
810 ierr = DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);CHKERRQ(ierr);
811 ierr = PetscSectionCreate(PetscObjectComm((PetscObject) dm), &msec);CHKERRQ(ierr);
812 ierr = PetscSectionSetNumFields(msec, 1);CHKERRQ(ierr);
813 ierr = PetscSectionSetFieldComponents(msec, 0, Nd);CHKERRQ(ierr);
814 ierr = PetscSectionSetChart(msec, vStart, vEnd);CHKERRQ(ierr);
815 for (v = vStart; v < vEnd; ++v) {
816 ierr = PetscSectionSetDof(msec, v, Nd);CHKERRQ(ierr);
817 ierr = PetscSectionSetFieldDof(msec, v, 0, Nd);CHKERRQ(ierr);
818 }
819 ierr = PetscSectionSetUp(msec);CHKERRQ(ierr);
820 ierr = DMSetLocalSection(dmMetric, msec);CHKERRQ(ierr);
821 ierr = PetscSectionDestroy(&msec);CHKERRQ(ierr);
822 ierr = DMGetLocalVector(dmMetric, &metric);CHKERRQ(ierr);
823 /* N is the target size */
824 N = adaptor->Nadapt >= 0 ? adaptor->Nadapt : PetscPowRealInt(adaptor->refinementFactor, dim)*((PetscReal) (vEnd - vStart));
825 if (adaptor->monitor) {ierr = PetscPrintf(PETSC_COMM_SELF, "N_orig: %D N_adapt: %g\n", vEnd - vStart, N);CHKERRQ(ierr);}
826 /* |H| means take the absolute value of eigenvalues */
827 ierr = VecGetArray(xHess, &H);CHKERRQ(ierr);
828 ierr = VecGetArray(metric, &M);CHKERRQ(ierr);
829 for (v = vStart; v < vEnd; ++v) {
830 PetscScalar *Hp;
831
832 ierr = DMPlexPointLocalRef(dmHess, v, H, &Hp);CHKERRQ(ierr);
833 ierr = DMAdaptorModifyHessian_Private(coordDim, adaptor->h_min, adaptor->h_max, Hp);CHKERRQ(ierr);
834 }
835 /* Pointwise on vertices M(x) = N^{2/d} (\int_\Omega det(|H|)^{p/(2p+d)})^{-2/d} det(|H|)^{-1/(2p+d)} |H| for L_p */
836 ierr = PetscDSSetObjective(probHess, 0, detHFunc);CHKERRQ(ierr);
837 ierr = DMPlexComputeIntegralFEM(dmHess, xHess, &integral, NULL);CHKERRQ(ierr);
838 for (v = vStart; v < vEnd; ++v) {
839 const PetscInt p = 1;
840 const PetscScalar *Hp;
841 PetscScalar *Mp;
842 PetscReal detH, fact;
843 PetscInt i;
844
845 ierr = DMPlexPointLocalRead(dmHess, v, H, (void *) &Hp);CHKERRQ(ierr);
846 ierr = DMPlexPointLocalRef(dmMetric, v, M, &Mp);CHKERRQ(ierr);
847 if (dim == 2) DMPlex_Det2D_Scalar_Internal(&detH, Hp);
848 else if (dim == 3) DMPlex_Det3D_Scalar_Internal(&detH, Hp);
849 else SETERRQ1(PetscObjectComm((PetscObject) adaptor), PETSC_ERR_SUP, "Dimension %d not supported", dim);
850 fact = PetscPowReal(N, 2.0/dim) * PetscPowReal(PetscRealPart(integral), -2.0/dim) * PetscPowReal(PetscAbsReal(detH), -1.0/(2*p+dim));
851 #if 0
852 ierr = PetscPrintf(PETSC_COMM_SELF, "fact: %g integral: %g |detH|: %g termA: %g termB: %g\n", fact, integral, PetscAbsReal(detH), PetscPowReal(integral, -2.0/dim), PetscPowReal(PetscAbsReal(detH), -1.0/(2*p+dim)));CHKERRQ(ierr);
853 ierr = DMPrintCellMatrix(v, "H", coordDim, coordDim, Hp);CHKERRQ(ierr);
854 #endif
855 for (i = 0; i < Nd; ++i) {
856 Mp[i] = fact * Hp[i];
857 }
858 }
859 ierr = VecRestoreArray(xHess, &H);CHKERRQ(ierr);
860 ierr = VecRestoreArray(metric, &M);CHKERRQ(ierr);
861 /* Adapt DM from metric */
862 ierr = DMGetLabel(dm, "marker", &bdLabel);CHKERRQ(ierr);
863 ierr = DMAdaptMetric(dm, metric, bdLabel, &odm);CHKERRQ(ierr);
864 adapted = PETSC_TRUE;
865 /* Cleanup */
866 ierr = DMRestoreLocalVector(dmMetric, &metric);CHKERRQ(ierr);
867 ierr = DMDestroy(&dmMetric);CHKERRQ(ierr);
868 ierr = DMRestoreGlobalVector(dmHess, &xHess);CHKERRQ(ierr);
869 ierr = DMRestoreGlobalVector(dmGrad, &xGrad);CHKERRQ(ierr);
870 ierr = DMDestroy(&dmGrad);CHKERRQ(ierr);
871 ierr = DMDestroy(&dmHess);CHKERRQ(ierr);
872 }
873 break;
874 default: SETERRQ1(comm, PETSC_ERR_ARG_WRONG, "Invalid adaptation type: %D", adaptor->adaptCriterion);
875 }
876 ierr = DMAdaptorPostAdapt(adaptor);CHKERRQ(ierr);
877 ierr = DMRestoreLocalVector(dm, &locX);CHKERRQ(ierr);
878 /* If DM was adapted, replace objects and recreate solution */
879 if (adapted) {
880 const char *name;
881
882 ierr = PetscObjectGetName((PetscObject) dm, &name);CHKERRQ(ierr);
883 ierr = PetscObjectSetName((PetscObject) odm, name);CHKERRQ(ierr);
884 /* Reconfigure solver */
885 ierr = SNESReset(adaptor->snes);CHKERRQ(ierr);
886 ierr = SNESSetDM(adaptor->snes, odm);CHKERRQ(ierr);
887 ierr = DMAdaptorSetSolver(adaptor, adaptor->snes);CHKERRQ(ierr);
888 ierr = DMPlexSetSNESLocalFEM(odm, ctx, ctx, ctx);CHKERRQ(ierr);
889 ierr = SNESSetFromOptions(adaptor->snes);CHKERRQ(ierr);
890 /* Transfer system */
891 ierr = DMCopyDisc(adaptor->idm, odm);CHKERRQ(ierr);
892 /* Transfer solution */
893 ierr = DMCreateGlobalVector(odm, &ox);CHKERRQ(ierr);
894 ierr = PetscObjectGetName((PetscObject) x, &name);CHKERRQ(ierr);
895 ierr = PetscObjectSetName((PetscObject) ox, name);CHKERRQ(ierr);
896 ierr = DMAdaptorTransferSolution(adaptor, dm, x, odm, ox);CHKERRQ(ierr);
897 /* Cleanup adaptivity info */
898 if (adaptIter > 0) {ierr = PetscViewerASCIIPopTab(PETSC_VIEWER_STDOUT_(comm));CHKERRQ(ierr);}
899 ierr = DMForestSetAdaptivityForest(dm, NULL);CHKERRQ(ierr); /* clear internal references to the previous dm */
900 ierr = DMDestroy(&dm);CHKERRQ(ierr);
901 ierr = VecDestroy(&x);CHKERRQ(ierr);
902 *adm = odm;
903 *ax = ox;
904 } else {
905 *adm = dm;
906 *ax = x;
907 adaptIter = numAdapt;
908 }
909 if (adaptIter < numAdapt-1) {
910 ierr = DMViewFromOptions(odm, NULL, "-dm_adapt_iter_view");CHKERRQ(ierr);
911 ierr = VecViewFromOptions(ox, NULL, "-sol_adapt_iter_view");CHKERRQ(ierr);
912 }
913 }
914 ierr = DMViewFromOptions(*adm, NULL, "-dm_adapt_view");CHKERRQ(ierr);
915 ierr = VecViewFromOptions(*ax, NULL, "-sol_adapt_view");CHKERRQ(ierr);
916 PetscFunctionReturn(0);
917 }
918
919 /*@
920 DMAdaptorAdapt - Creates a new DM that is adapted to the problem
921
922 Not collective
923
924 Input Parameter:
925 + adaptor - The DMAdaptor object
926 . x - The global approximate solution
927 - strategy - The adaptation strategy
928
929 Output Parameters:
930 + adm - The adapted DM
931 - ax - The adapted solution
932
933 Options database keys:
934 + -snes_adapt <strategy> : initial, sequential, multigrid
935 . -adapt_gradient_view : View the Clement interpolant of the solution gradient
936 . -adapt_hessian_view : View the Clement interpolant of the solution Hessian
937 - -adapt_metric_view : View the metric tensor for adaptive mesh refinement
938
939 Note: The available adaptation strategies are:
940 $ 1) Adapt the initial mesh until a quality metric, e.g., a priori error bound, is satisfied
941 $ 2) Solve the problem on a series of adapted meshes until a quality metric, e.g. a posteriori error bound, is satisfied
942 $ 3) Solve the problem on a hierarchy of adapted meshes generated to satisfy a quality metric using multigrid
943
944 Level: intermediate
945
946 .seealso: DMAdaptorSetSolver(), DMAdaptorCreate(), DMAdaptorAdapt()
947 @*/
DMAdaptorAdapt(DMAdaptor adaptor,Vec x,DMAdaptationStrategy strategy,DM * adm,Vec * ax)948 PetscErrorCode DMAdaptorAdapt(DMAdaptor adaptor, Vec x, DMAdaptationStrategy strategy, DM *adm, Vec *ax)
949 {
950 PetscErrorCode ierr;
951
952 PetscFunctionBegin;
953 switch (strategy)
954 {
955 case DM_ADAPTATION_INITIAL:
956 ierr = DMAdaptorAdapt_Sequence_Private(adaptor, x, PETSC_FALSE, adm, ax);CHKERRQ(ierr);
957 break;
958 case DM_ADAPTATION_SEQUENTIAL:
959 ierr = DMAdaptorAdapt_Sequence_Private(adaptor, x, PETSC_TRUE, adm, ax);CHKERRQ(ierr);
960 break;
961 default: SETERRQ1(PetscObjectComm((PetscObject) adaptor), PETSC_ERR_ARG_WRONG, "Unrecognized adaptation strategy %d", strategy);
962 }
963 PetscFunctionReturn(0);
964 }
965