1 #include <petsc/private/dmpleximpl.h> /*I "petscdmplex.h" I*/
2 #include <petsc/private/snesimpl.h> /*I "petscsnes.h" I*/
3 #include <petscds.h>
4 #include <petscblaslapack.h>
5 #include <petsc/private/petscimpl.h>
6 #include <petsc/private/petscfeimpl.h>
7
8 /************************** Interpolation *******************************/
9
DMSNESConvertPlex(DM dm,DM * plex,PetscBool copy)10 static PetscErrorCode DMSNESConvertPlex(DM dm, DM *plex, PetscBool copy)
11 {
12 PetscBool isPlex;
13 PetscErrorCode ierr;
14
15 PetscFunctionBegin;
16 ierr = PetscObjectTypeCompare((PetscObject) dm, DMPLEX, &isPlex);CHKERRQ(ierr);
17 if (isPlex) {
18 *plex = dm;
19 ierr = PetscObjectReference((PetscObject) dm);CHKERRQ(ierr);
20 } else {
21 ierr = PetscObjectQuery((PetscObject) dm, "dm_plex", (PetscObject *) plex);CHKERRQ(ierr);
22 if (!*plex) {
23 ierr = DMConvert(dm,DMPLEX,plex);CHKERRQ(ierr);
24 ierr = PetscObjectCompose((PetscObject) dm, "dm_plex", (PetscObject) *plex);CHKERRQ(ierr);
25 if (copy) {
26 PetscInt i;
27 PetscObject obj;
28 const char *comps[3] = {"A","dmAux","dmCh"};
29
30 ierr = DMCopyDMSNES(dm, *plex);CHKERRQ(ierr);
31 for (i = 0; i < 3; i++) {
32 ierr = PetscObjectQuery((PetscObject) dm, comps[i], &obj);CHKERRQ(ierr);
33 ierr = PetscObjectCompose((PetscObject) *plex, comps[i], obj);CHKERRQ(ierr);
34 }
35 }
36 } else {
37 ierr = PetscObjectReference((PetscObject) *plex);CHKERRQ(ierr);
38 }
39 }
40 PetscFunctionReturn(0);
41 }
42
43 /*@C
44 DMInterpolationCreate - Creates a DMInterpolationInfo context
45
46 Collective
47
48 Input Parameter:
49 . comm - the communicator
50
51 Output Parameter:
52 . ctx - the context
53
54 Level: beginner
55
56 .seealso: DMInterpolationEvaluate(), DMInterpolationAddPoints(), DMInterpolationDestroy()
57 @*/
DMInterpolationCreate(MPI_Comm comm,DMInterpolationInfo * ctx)58 PetscErrorCode DMInterpolationCreate(MPI_Comm comm, DMInterpolationInfo *ctx)
59 {
60 PetscErrorCode ierr;
61
62 PetscFunctionBegin;
63 PetscValidPointer(ctx, 2);
64 ierr = PetscNew(ctx);CHKERRQ(ierr);
65
66 (*ctx)->comm = comm;
67 (*ctx)->dim = -1;
68 (*ctx)->nInput = 0;
69 (*ctx)->points = NULL;
70 (*ctx)->cells = NULL;
71 (*ctx)->n = -1;
72 (*ctx)->coords = NULL;
73 PetscFunctionReturn(0);
74 }
75
76 /*@C
77 DMInterpolationSetDim - Sets the spatial dimension for the interpolation context
78
79 Not collective
80
81 Input Parameters:
82 + ctx - the context
83 - dim - the spatial dimension
84
85 Level: intermediate
86
87 .seealso: DMInterpolationGetDim(), DMInterpolationEvaluate(), DMInterpolationAddPoints()
88 @*/
DMInterpolationSetDim(DMInterpolationInfo ctx,PetscInt dim)89 PetscErrorCode DMInterpolationSetDim(DMInterpolationInfo ctx, PetscInt dim)
90 {
91 PetscFunctionBegin;
92 if ((dim < 1) || (dim > 3)) SETERRQ1(ctx->comm, PETSC_ERR_ARG_OUTOFRANGE, "Invalid dimension for points: %D", dim);
93 ctx->dim = dim;
94 PetscFunctionReturn(0);
95 }
96
97 /*@C
98 DMInterpolationGetDim - Gets the spatial dimension for the interpolation context
99
100 Not collective
101
102 Input Parameter:
103 . ctx - the context
104
105 Output Parameter:
106 . dim - the spatial dimension
107
108 Level: intermediate
109
110 .seealso: DMInterpolationSetDim(), DMInterpolationEvaluate(), DMInterpolationAddPoints()
111 @*/
DMInterpolationGetDim(DMInterpolationInfo ctx,PetscInt * dim)112 PetscErrorCode DMInterpolationGetDim(DMInterpolationInfo ctx, PetscInt *dim)
113 {
114 PetscFunctionBegin;
115 PetscValidIntPointer(dim, 2);
116 *dim = ctx->dim;
117 PetscFunctionReturn(0);
118 }
119
120 /*@C
121 DMInterpolationSetDof - Sets the number of fields interpolated at a point for the interpolation context
122
123 Not collective
124
125 Input Parameters:
126 + ctx - the context
127 - dof - the number of fields
128
129 Level: intermediate
130
131 .seealso: DMInterpolationGetDof(), DMInterpolationEvaluate(), DMInterpolationAddPoints()
132 @*/
DMInterpolationSetDof(DMInterpolationInfo ctx,PetscInt dof)133 PetscErrorCode DMInterpolationSetDof(DMInterpolationInfo ctx, PetscInt dof)
134 {
135 PetscFunctionBegin;
136 if (dof < 1) SETERRQ1(ctx->comm, PETSC_ERR_ARG_OUTOFRANGE, "Invalid number of components: %D", dof);
137 ctx->dof = dof;
138 PetscFunctionReturn(0);
139 }
140
141 /*@C
142 DMInterpolationGetDof - Gets the number of fields interpolated at a point for the interpolation context
143
144 Not collective
145
146 Input Parameter:
147 . ctx - the context
148
149 Output Parameter:
150 . dof - the number of fields
151
152 Level: intermediate
153
154 .seealso: DMInterpolationSetDof(), DMInterpolationEvaluate(), DMInterpolationAddPoints()
155 @*/
DMInterpolationGetDof(DMInterpolationInfo ctx,PetscInt * dof)156 PetscErrorCode DMInterpolationGetDof(DMInterpolationInfo ctx, PetscInt *dof)
157 {
158 PetscFunctionBegin;
159 PetscValidIntPointer(dof, 2);
160 *dof = ctx->dof;
161 PetscFunctionReturn(0);
162 }
163
164 /*@C
165 DMInterpolationAddPoints - Add points at which we will interpolate the fields
166
167 Not collective
168
169 Input Parameters:
170 + ctx - the context
171 . n - the number of points
172 - points - the coordinates for each point, an array of size n * dim
173
174 Note: The coordinate information is copied.
175
176 Level: intermediate
177
178 .seealso: DMInterpolationSetDim(), DMInterpolationEvaluate(), DMInterpolationCreate()
179 @*/
DMInterpolationAddPoints(DMInterpolationInfo ctx,PetscInt n,PetscReal points[])180 PetscErrorCode DMInterpolationAddPoints(DMInterpolationInfo ctx, PetscInt n, PetscReal points[])
181 {
182 PetscErrorCode ierr;
183
184 PetscFunctionBegin;
185 if (ctx->dim < 0) SETERRQ(ctx->comm, PETSC_ERR_ARG_WRONGSTATE, "The spatial dimension has not been set");
186 if (ctx->points) SETERRQ(ctx->comm, PETSC_ERR_ARG_WRONGSTATE, "Cannot add points multiple times yet");
187 ctx->nInput = n;
188
189 ierr = PetscMalloc1(n*ctx->dim, &ctx->points);CHKERRQ(ierr);
190 ierr = PetscArraycpy(ctx->points, points, n*ctx->dim);CHKERRQ(ierr);
191 PetscFunctionReturn(0);
192 }
193
194 /*@C
195 DMInterpolationSetUp - Computea spatial indices that add in point location during interpolation
196
197 Collective on ctx
198
199 Input Parameters:
200 + ctx - the context
201 . dm - the DM for the function space used for interpolation
202 - redundantPoints - If PETSC_TRUE, all processes are passing in the same array of points. Otherwise, points need to be communicated among processes.
203
204 Level: intermediate
205
206 .seealso: DMInterpolationEvaluate(), DMInterpolationAddPoints(), DMInterpolationCreate()
207 @*/
DMInterpolationSetUp(DMInterpolationInfo ctx,DM dm,PetscBool redundantPoints)208 PetscErrorCode DMInterpolationSetUp(DMInterpolationInfo ctx, DM dm, PetscBool redundantPoints)
209 {
210 MPI_Comm comm = ctx->comm;
211 PetscScalar *a;
212 PetscInt p, q, i;
213 PetscMPIInt rank, size;
214 PetscErrorCode ierr;
215 Vec pointVec;
216 PetscSF cellSF;
217 PetscLayout layout;
218 PetscReal *globalPoints;
219 PetscScalar *globalPointsScalar;
220 const PetscInt *ranges;
221 PetscMPIInt *counts, *displs;
222 const PetscSFNode *foundCells;
223 const PetscInt *foundPoints;
224 PetscMPIInt *foundProcs, *globalProcs;
225 PetscInt n, N, numFound;
226
227 PetscFunctionBegin;
228 PetscValidHeaderSpecific(dm, DM_CLASSID, 1);
229 ierr = MPI_Comm_size(comm, &size);CHKERRQ(ierr);
230 ierr = MPI_Comm_rank(comm, &rank);CHKERRQ(ierr);
231 if (ctx->dim < 0) SETERRQ(comm, PETSC_ERR_ARG_WRONGSTATE, "The spatial dimension has not been set");
232 /* Locate points */
233 n = ctx->nInput;
234 if (!redundantPoints) {
235 ierr = PetscLayoutCreate(comm, &layout);CHKERRQ(ierr);
236 ierr = PetscLayoutSetBlockSize(layout, 1);CHKERRQ(ierr);
237 ierr = PetscLayoutSetLocalSize(layout, n);CHKERRQ(ierr);
238 ierr = PetscLayoutSetUp(layout);CHKERRQ(ierr);
239 ierr = PetscLayoutGetSize(layout, &N);CHKERRQ(ierr);
240 /* Communicate all points to all processes */
241 ierr = PetscMalloc3(N*ctx->dim,&globalPoints,size,&counts,size,&displs);CHKERRQ(ierr);
242 ierr = PetscLayoutGetRanges(layout, &ranges);CHKERRQ(ierr);
243 for (p = 0; p < size; ++p) {
244 counts[p] = (ranges[p+1] - ranges[p])*ctx->dim;
245 displs[p] = ranges[p]*ctx->dim;
246 }
247 ierr = MPI_Allgatherv(ctx->points, n*ctx->dim, MPIU_REAL, globalPoints, counts, displs, MPIU_REAL, comm);CHKERRQ(ierr);
248 } else {
249 N = n;
250 globalPoints = ctx->points;
251 counts = displs = NULL;
252 layout = NULL;
253 }
254 #if 0
255 ierr = PetscMalloc3(N,&foundCells,N,&foundProcs,N,&globalProcs);CHKERRQ(ierr);
256 /* foundCells[p] = m->locatePoint(&globalPoints[p*ctx->dim]); */
257 #else
258 #if defined(PETSC_USE_COMPLEX)
259 ierr = PetscMalloc1(N*ctx->dim,&globalPointsScalar);CHKERRQ(ierr);
260 for (i=0; i<N*ctx->dim; i++) globalPointsScalar[i] = globalPoints[i];
261 #else
262 globalPointsScalar = globalPoints;
263 #endif
264 ierr = VecCreateSeqWithArray(PETSC_COMM_SELF, ctx->dim, N*ctx->dim, globalPointsScalar, &pointVec);CHKERRQ(ierr);
265 ierr = PetscMalloc2(N,&foundProcs,N,&globalProcs);CHKERRQ(ierr);
266 for (p = 0; p < N; ++p) {foundProcs[p] = size;}
267 cellSF = NULL;
268 ierr = DMLocatePoints(dm, pointVec, DM_POINTLOCATION_REMOVE, &cellSF);CHKERRQ(ierr);
269 ierr = PetscSFGetGraph(cellSF,NULL,&numFound,&foundPoints,&foundCells);CHKERRQ(ierr);
270 #endif
271 for (p = 0; p < numFound; ++p) {
272 if (foundCells[p].index >= 0) foundProcs[foundPoints ? foundPoints[p] : p] = rank;
273 }
274 /* Let the lowest rank process own each point */
275 ierr = MPIU_Allreduce(foundProcs, globalProcs, N, MPI_INT, MPI_MIN, comm);CHKERRQ(ierr);
276 ctx->n = 0;
277 for (p = 0; p < N; ++p) {
278 if (globalProcs[p] == size) SETERRQ4(comm, PETSC_ERR_PLIB, "Point %d: %g %g %g not located in mesh", p, (double)globalPoints[p*ctx->dim+0], (double)(ctx->dim > 1 ? globalPoints[p*ctx->dim+1] : 0.0), (double)(ctx->dim > 2 ? globalPoints[p*ctx->dim+2] : 0.0));
279 else if (globalProcs[p] == rank) ctx->n++;
280 }
281 /* Create coordinates vector and array of owned cells */
282 ierr = PetscMalloc1(ctx->n, &ctx->cells);CHKERRQ(ierr);
283 ierr = VecCreate(comm, &ctx->coords);CHKERRQ(ierr);
284 ierr = VecSetSizes(ctx->coords, ctx->n*ctx->dim, PETSC_DECIDE);CHKERRQ(ierr);
285 ierr = VecSetBlockSize(ctx->coords, ctx->dim);CHKERRQ(ierr);
286 ierr = VecSetType(ctx->coords,VECSTANDARD);CHKERRQ(ierr);
287 ierr = VecGetArray(ctx->coords, &a);CHKERRQ(ierr);
288 for (p = 0, q = 0, i = 0; p < N; ++p) {
289 if (globalProcs[p] == rank) {
290 PetscInt d;
291
292 for (d = 0; d < ctx->dim; ++d, ++i) a[i] = globalPoints[p*ctx->dim+d];
293 ctx->cells[q] = foundCells[q].index;
294 ++q;
295 }
296 }
297 ierr = VecRestoreArray(ctx->coords, &a);CHKERRQ(ierr);
298 #if 0
299 ierr = PetscFree3(foundCells,foundProcs,globalProcs);CHKERRQ(ierr);
300 #else
301 ierr = PetscFree2(foundProcs,globalProcs);CHKERRQ(ierr);
302 ierr = PetscSFDestroy(&cellSF);CHKERRQ(ierr);
303 ierr = VecDestroy(&pointVec);CHKERRQ(ierr);
304 #endif
305 if ((void*)globalPointsScalar != (void*)globalPoints) {ierr = PetscFree(globalPointsScalar);CHKERRQ(ierr);}
306 if (!redundantPoints) {ierr = PetscFree3(globalPoints,counts,displs);CHKERRQ(ierr);}
307 ierr = PetscLayoutDestroy(&layout);CHKERRQ(ierr);
308 PetscFunctionReturn(0);
309 }
310
311 /*@C
312 DMInterpolationGetCoordinates - Gets a Vec with the coordinates of each interpolation point
313
314 Collective on ctx
315
316 Input Parameter:
317 . ctx - the context
318
319 Output Parameter:
320 . coordinates - the coordinates of interpolation points
321
322 Note: The local vector entries correspond to interpolation points lying on this process, according to the associated DM. This is a borrowed vector that the user should not destroy.
323
324 Level: intermediate
325
326 .seealso: DMInterpolationEvaluate(), DMInterpolationAddPoints(), DMInterpolationCreate()
327 @*/
DMInterpolationGetCoordinates(DMInterpolationInfo ctx,Vec * coordinates)328 PetscErrorCode DMInterpolationGetCoordinates(DMInterpolationInfo ctx, Vec *coordinates)
329 {
330 PetscFunctionBegin;
331 PetscValidPointer(coordinates, 2);
332 if (!ctx->coords) SETERRQ(ctx->comm, PETSC_ERR_ARG_WRONGSTATE, "The interpolation context has not been setup.");
333 *coordinates = ctx->coords;
334 PetscFunctionReturn(0);
335 }
336
337 /*@C
338 DMInterpolationGetVector - Gets a Vec which can hold all the interpolated field values
339
340 Collective on ctx
341
342 Input Parameter:
343 . ctx - the context
344
345 Output Parameter:
346 . v - a vector capable of holding the interpolated field values
347
348 Note: This vector should be returned using DMInterpolationRestoreVector().
349
350 Level: intermediate
351
352 .seealso: DMInterpolationRestoreVector(), DMInterpolationEvaluate(), DMInterpolationAddPoints(), DMInterpolationCreate()
353 @*/
DMInterpolationGetVector(DMInterpolationInfo ctx,Vec * v)354 PetscErrorCode DMInterpolationGetVector(DMInterpolationInfo ctx, Vec *v)
355 {
356 PetscErrorCode ierr;
357
358 PetscFunctionBegin;
359 PetscValidPointer(v, 2);
360 if (!ctx->coords) SETERRQ(ctx->comm, PETSC_ERR_ARG_WRONGSTATE, "The interpolation context has not been setup.");
361 ierr = VecCreate(ctx->comm, v);CHKERRQ(ierr);
362 ierr = VecSetSizes(*v, ctx->n*ctx->dof, PETSC_DECIDE);CHKERRQ(ierr);
363 ierr = VecSetBlockSize(*v, ctx->dof);CHKERRQ(ierr);
364 ierr = VecSetType(*v,VECSTANDARD);CHKERRQ(ierr);
365 PetscFunctionReturn(0);
366 }
367
368 /*@C
369 DMInterpolationRestoreVector - Returns a Vec which can hold all the interpolated field values
370
371 Collective on ctx
372
373 Input Parameters:
374 + ctx - the context
375 - v - a vector capable of holding the interpolated field values
376
377 Level: intermediate
378
379 .seealso: DMInterpolationGetVector(), DMInterpolationEvaluate(), DMInterpolationAddPoints(), DMInterpolationCreate()
380 @*/
DMInterpolationRestoreVector(DMInterpolationInfo ctx,Vec * v)381 PetscErrorCode DMInterpolationRestoreVector(DMInterpolationInfo ctx, Vec *v)
382 {
383 PetscErrorCode ierr;
384
385 PetscFunctionBegin;
386 PetscValidPointer(v, 2);
387 if (!ctx->coords) SETERRQ(ctx->comm, PETSC_ERR_ARG_WRONGSTATE, "The interpolation context has not been setup.");
388 ierr = VecDestroy(v);CHKERRQ(ierr);
389 PetscFunctionReturn(0);
390 }
391
DMInterpolate_Triangle_Private(DMInterpolationInfo ctx,DM dm,Vec xLocal,Vec v)392 PETSC_STATIC_INLINE PetscErrorCode DMInterpolate_Triangle_Private(DMInterpolationInfo ctx, DM dm, Vec xLocal, Vec v)
393 {
394 PetscReal *v0, *J, *invJ, detJ;
395 const PetscScalar *coords;
396 PetscScalar *a;
397 PetscInt p;
398 PetscErrorCode ierr;
399
400 PetscFunctionBegin;
401 ierr = PetscMalloc3(ctx->dim,&v0,ctx->dim*ctx->dim,&J,ctx->dim*ctx->dim,&invJ);CHKERRQ(ierr);
402 ierr = VecGetArrayRead(ctx->coords, &coords);CHKERRQ(ierr);
403 ierr = VecGetArray(v, &a);CHKERRQ(ierr);
404 for (p = 0; p < ctx->n; ++p) {
405 PetscInt c = ctx->cells[p];
406 PetscScalar *x = NULL;
407 PetscReal xi[4];
408 PetscInt d, f, comp;
409
410 ierr = DMPlexComputeCellGeometryFEM(dm, c, NULL, v0, J, invJ, &detJ);CHKERRQ(ierr);
411 if (detJ <= 0.0) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid determinant %g for element %D", (double)detJ, c);
412 ierr = DMPlexVecGetClosure(dm, NULL, xLocal, c, NULL, &x);CHKERRQ(ierr);
413 for (comp = 0; comp < ctx->dof; ++comp) a[p*ctx->dof+comp] = x[0*ctx->dof+comp];
414
415 for (d = 0; d < ctx->dim; ++d) {
416 xi[d] = 0.0;
417 for (f = 0; f < ctx->dim; ++f) xi[d] += invJ[d*ctx->dim+f]*0.5*PetscRealPart(coords[p*ctx->dim+f] - v0[f]);
418 for (comp = 0; comp < ctx->dof; ++comp) a[p*ctx->dof+comp] += PetscRealPart(x[(d+1)*ctx->dof+comp] - x[0*ctx->dof+comp])*xi[d];
419 }
420 ierr = DMPlexVecRestoreClosure(dm, NULL, xLocal, c, NULL, &x);CHKERRQ(ierr);
421 }
422 ierr = VecRestoreArray(v, &a);CHKERRQ(ierr);
423 ierr = VecRestoreArrayRead(ctx->coords, &coords);CHKERRQ(ierr);
424 ierr = PetscFree3(v0, J, invJ);CHKERRQ(ierr);
425 PetscFunctionReturn(0);
426 }
427
DMInterpolate_Tetrahedron_Private(DMInterpolationInfo ctx,DM dm,Vec xLocal,Vec v)428 PETSC_STATIC_INLINE PetscErrorCode DMInterpolate_Tetrahedron_Private(DMInterpolationInfo ctx, DM dm, Vec xLocal, Vec v)
429 {
430 PetscReal *v0, *J, *invJ, detJ;
431 const PetscScalar *coords;
432 PetscScalar *a;
433 PetscInt p;
434 PetscErrorCode ierr;
435
436 PetscFunctionBegin;
437 ierr = PetscMalloc3(ctx->dim,&v0,ctx->dim*ctx->dim,&J,ctx->dim*ctx->dim,&invJ);CHKERRQ(ierr);
438 ierr = VecGetArrayRead(ctx->coords, &coords);CHKERRQ(ierr);
439 ierr = VecGetArray(v, &a);CHKERRQ(ierr);
440 for (p = 0; p < ctx->n; ++p) {
441 PetscInt c = ctx->cells[p];
442 const PetscInt order[3] = {2, 1, 3};
443 PetscScalar *x = NULL;
444 PetscReal xi[4];
445 PetscInt d, f, comp;
446
447 ierr = DMPlexComputeCellGeometryFEM(dm, c, NULL, v0, J, invJ, &detJ);CHKERRQ(ierr);
448 if (detJ <= 0.0) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid determinant %g for element %D", (double)detJ, c);
449 ierr = DMPlexVecGetClosure(dm, NULL, xLocal, c, NULL, &x);CHKERRQ(ierr);
450 for (comp = 0; comp < ctx->dof; ++comp) a[p*ctx->dof+comp] = x[0*ctx->dof+comp];
451
452 for (d = 0; d < ctx->dim; ++d) {
453 xi[d] = 0.0;
454 for (f = 0; f < ctx->dim; ++f) xi[d] += invJ[d*ctx->dim+f]*0.5*PetscRealPart(coords[p*ctx->dim+f] - v0[f]);
455 for (comp = 0; comp < ctx->dof; ++comp) a[p*ctx->dof+comp] += PetscRealPart(x[order[d]*ctx->dof+comp] - x[0*ctx->dof+comp])*xi[d];
456 }
457 ierr = DMPlexVecRestoreClosure(dm, NULL, xLocal, c, NULL, &x);CHKERRQ(ierr);
458 }
459 ierr = VecRestoreArray(v, &a);CHKERRQ(ierr);
460 ierr = VecRestoreArrayRead(ctx->coords, &coords);CHKERRQ(ierr);
461 ierr = PetscFree3(v0, J, invJ);CHKERRQ(ierr);
462 PetscFunctionReturn(0);
463 }
464
QuadMap_Private(SNES snes,Vec Xref,Vec Xreal,void * ctx)465 PETSC_STATIC_INLINE PetscErrorCode QuadMap_Private(SNES snes, Vec Xref, Vec Xreal, void *ctx)
466 {
467 const PetscScalar *vertices = (const PetscScalar*) ctx;
468 const PetscScalar x0 = vertices[0];
469 const PetscScalar y0 = vertices[1];
470 const PetscScalar x1 = vertices[2];
471 const PetscScalar y1 = vertices[3];
472 const PetscScalar x2 = vertices[4];
473 const PetscScalar y2 = vertices[5];
474 const PetscScalar x3 = vertices[6];
475 const PetscScalar y3 = vertices[7];
476 const PetscScalar f_1 = x1 - x0;
477 const PetscScalar g_1 = y1 - y0;
478 const PetscScalar f_3 = x3 - x0;
479 const PetscScalar g_3 = y3 - y0;
480 const PetscScalar f_01 = x2 - x1 - x3 + x0;
481 const PetscScalar g_01 = y2 - y1 - y3 + y0;
482 const PetscScalar *ref;
483 PetscScalar *real;
484 PetscErrorCode ierr;
485
486 PetscFunctionBegin;
487 ierr = VecGetArrayRead(Xref, &ref);CHKERRQ(ierr);
488 ierr = VecGetArray(Xreal, &real);CHKERRQ(ierr);
489 {
490 const PetscScalar p0 = ref[0];
491 const PetscScalar p1 = ref[1];
492
493 real[0] = x0 + f_1 * p0 + f_3 * p1 + f_01 * p0 * p1;
494 real[1] = y0 + g_1 * p0 + g_3 * p1 + g_01 * p0 * p1;
495 }
496 ierr = PetscLogFlops(28);CHKERRQ(ierr);
497 ierr = VecRestoreArrayRead(Xref, &ref);CHKERRQ(ierr);
498 ierr = VecRestoreArray(Xreal, &real);CHKERRQ(ierr);
499 PetscFunctionReturn(0);
500 }
501
502 #include <petsc/private/dmimpl.h>
QuadJacobian_Private(SNES snes,Vec Xref,Mat J,Mat M,void * ctx)503 PETSC_STATIC_INLINE PetscErrorCode QuadJacobian_Private(SNES snes, Vec Xref, Mat J, Mat M, void *ctx)
504 {
505 const PetscScalar *vertices = (const PetscScalar*) ctx;
506 const PetscScalar x0 = vertices[0];
507 const PetscScalar y0 = vertices[1];
508 const PetscScalar x1 = vertices[2];
509 const PetscScalar y1 = vertices[3];
510 const PetscScalar x2 = vertices[4];
511 const PetscScalar y2 = vertices[5];
512 const PetscScalar x3 = vertices[6];
513 const PetscScalar y3 = vertices[7];
514 const PetscScalar f_01 = x2 - x1 - x3 + x0;
515 const PetscScalar g_01 = y2 - y1 - y3 + y0;
516 const PetscScalar *ref;
517 PetscErrorCode ierr;
518
519 PetscFunctionBegin;
520 ierr = VecGetArrayRead(Xref, &ref);CHKERRQ(ierr);
521 {
522 const PetscScalar x = ref[0];
523 const PetscScalar y = ref[1];
524 const PetscInt rows[2] = {0, 1};
525 PetscScalar values[4];
526
527 values[0] = (x1 - x0 + f_01*y) * 0.5; values[1] = (x3 - x0 + f_01*x) * 0.5;
528 values[2] = (y1 - y0 + g_01*y) * 0.5; values[3] = (y3 - y0 + g_01*x) * 0.5;
529 ierr = MatSetValues(J, 2, rows, 2, rows, values, INSERT_VALUES);CHKERRQ(ierr);
530 }
531 ierr = PetscLogFlops(30);CHKERRQ(ierr);
532 ierr = VecRestoreArrayRead(Xref, &ref);CHKERRQ(ierr);
533 ierr = MatAssemblyBegin(J, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
534 ierr = MatAssemblyEnd(J, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
535 PetscFunctionReturn(0);
536 }
537
DMInterpolate_Quad_Private(DMInterpolationInfo ctx,DM dm,Vec xLocal,Vec v)538 PETSC_STATIC_INLINE PetscErrorCode DMInterpolate_Quad_Private(DMInterpolationInfo ctx, DM dm, Vec xLocal, Vec v)
539 {
540 DM dmCoord;
541 PetscFE fem = NULL;
542 SNES snes;
543 KSP ksp;
544 PC pc;
545 Vec coordsLocal, r, ref, real;
546 Mat J;
547 PetscTabulation T;
548 const PetscScalar *coords;
549 PetscScalar *a;
550 PetscReal xir[2];
551 PetscInt Nf, p;
552 const PetscInt dof = ctx->dof;
553 PetscErrorCode ierr;
554
555 PetscFunctionBegin;
556 ierr = DMGetNumFields(dm, &Nf);CHKERRQ(ierr);
557 if (Nf) {ierr = DMGetField(dm, 0, NULL, (PetscObject *) &fem);CHKERRQ(ierr);}
558 ierr = DMGetCoordinatesLocal(dm, &coordsLocal);CHKERRQ(ierr);
559 ierr = DMGetCoordinateDM(dm, &dmCoord);CHKERRQ(ierr);
560 ierr = SNESCreate(PETSC_COMM_SELF, &snes);CHKERRQ(ierr);
561 ierr = SNESSetOptionsPrefix(snes, "quad_interp_");CHKERRQ(ierr);
562 ierr = VecCreate(PETSC_COMM_SELF, &r);CHKERRQ(ierr);
563 ierr = VecSetSizes(r, 2, 2);CHKERRQ(ierr);
564 ierr = VecSetType(r,dm->vectype);CHKERRQ(ierr);
565 ierr = VecDuplicate(r, &ref);CHKERRQ(ierr);
566 ierr = VecDuplicate(r, &real);CHKERRQ(ierr);
567 ierr = MatCreate(PETSC_COMM_SELF, &J);CHKERRQ(ierr);
568 ierr = MatSetSizes(J, 2, 2, 2, 2);CHKERRQ(ierr);
569 ierr = MatSetType(J, MATSEQDENSE);CHKERRQ(ierr);
570 ierr = MatSetUp(J);CHKERRQ(ierr);
571 ierr = SNESSetFunction(snes, r, QuadMap_Private, NULL);CHKERRQ(ierr);
572 ierr = SNESSetJacobian(snes, J, J, QuadJacobian_Private, NULL);CHKERRQ(ierr);
573 ierr = SNESGetKSP(snes, &ksp);CHKERRQ(ierr);
574 ierr = KSPGetPC(ksp, &pc);CHKERRQ(ierr);
575 ierr = PCSetType(pc, PCLU);CHKERRQ(ierr);
576 ierr = SNESSetFromOptions(snes);CHKERRQ(ierr);
577
578 ierr = VecGetArrayRead(ctx->coords, &coords);CHKERRQ(ierr);
579 ierr = VecGetArray(v, &a);CHKERRQ(ierr);
580 ierr = PetscFECreateTabulation(fem, 1, 1, xir, 0, &T);CHKERRQ(ierr);
581 for (p = 0; p < ctx->n; ++p) {
582 PetscScalar *x = NULL, *vertices = NULL;
583 PetscScalar *xi;
584 PetscInt c = ctx->cells[p], comp, coordSize, xSize;
585
586 /* Can make this do all points at once */
587 ierr = DMPlexVecGetClosure(dmCoord, NULL, coordsLocal, c, &coordSize, &vertices);CHKERRQ(ierr);
588 if (4*2 != coordSize) SETERRQ2(ctx->comm, PETSC_ERR_ARG_SIZ, "Invalid closure size %D should be %d", coordSize, 4*2);
589 ierr = DMPlexVecGetClosure(dm, NULL, xLocal, c, &xSize, &x);CHKERRQ(ierr);
590 ierr = SNESSetFunction(snes, NULL, NULL, (void*) vertices);CHKERRQ(ierr);
591 ierr = SNESSetJacobian(snes, NULL, NULL, NULL, (void*) vertices);CHKERRQ(ierr);
592 ierr = VecGetArray(real, &xi);CHKERRQ(ierr);
593 xi[0] = coords[p*ctx->dim+0];
594 xi[1] = coords[p*ctx->dim+1];
595 ierr = VecRestoreArray(real, &xi);CHKERRQ(ierr);
596 ierr = SNESSolve(snes, real, ref);CHKERRQ(ierr);
597 ierr = VecGetArray(ref, &xi);CHKERRQ(ierr);
598 xir[0] = PetscRealPart(xi[0]);
599 xir[1] = PetscRealPart(xi[1]);
600 if (4*dof != xSize) {
601 PetscInt d;
602
603 xir[0] = 2.0*xir[0] - 1.0; xir[1] = 2.0*xir[1] - 1.0;
604 ierr = PetscFEComputeTabulation(fem, 1, xir, 0, T);CHKERRQ(ierr);
605 for (comp = 0; comp < dof; ++comp) {
606 a[p*dof+comp] = 0.0;
607 for (d = 0; d < xSize/dof; ++d) {
608 a[p*dof+comp] += x[d*dof+comp]*T->T[0][d*dof+comp];
609 }
610 }
611 } else {
612 for (comp = 0; comp < dof; ++comp)
613 a[p*dof+comp] = x[0*dof+comp]*(1 - xir[0])*(1 - xir[1]) + x[1*dof+comp]*xir[0]*(1 - xir[1]) + x[2*dof+comp]*xir[0]*xir[1] + x[3*dof+comp]*(1 - xir[0])*xir[1];
614 }
615 ierr = VecRestoreArray(ref, &xi);CHKERRQ(ierr);
616 ierr = DMPlexVecRestoreClosure(dmCoord, NULL, coordsLocal, c, &coordSize, &vertices);CHKERRQ(ierr);
617 ierr = DMPlexVecRestoreClosure(dm, NULL, xLocal, c, &xSize, &x);CHKERRQ(ierr);
618 }
619 ierr = PetscTabulationDestroy(&T);CHKERRQ(ierr);
620 ierr = VecRestoreArray(v, &a);CHKERRQ(ierr);
621 ierr = VecRestoreArrayRead(ctx->coords, &coords);CHKERRQ(ierr);
622
623 ierr = SNESDestroy(&snes);CHKERRQ(ierr);
624 ierr = VecDestroy(&r);CHKERRQ(ierr);
625 ierr = VecDestroy(&ref);CHKERRQ(ierr);
626 ierr = VecDestroy(&real);CHKERRQ(ierr);
627 ierr = MatDestroy(&J);CHKERRQ(ierr);
628 PetscFunctionReturn(0);
629 }
630
HexMap_Private(SNES snes,Vec Xref,Vec Xreal,void * ctx)631 PETSC_STATIC_INLINE PetscErrorCode HexMap_Private(SNES snes, Vec Xref, Vec Xreal, void *ctx)
632 {
633 const PetscScalar *vertices = (const PetscScalar*) ctx;
634 const PetscScalar x0 = vertices[0];
635 const PetscScalar y0 = vertices[1];
636 const PetscScalar z0 = vertices[2];
637 const PetscScalar x1 = vertices[9];
638 const PetscScalar y1 = vertices[10];
639 const PetscScalar z1 = vertices[11];
640 const PetscScalar x2 = vertices[6];
641 const PetscScalar y2 = vertices[7];
642 const PetscScalar z2 = vertices[8];
643 const PetscScalar x3 = vertices[3];
644 const PetscScalar y3 = vertices[4];
645 const PetscScalar z3 = vertices[5];
646 const PetscScalar x4 = vertices[12];
647 const PetscScalar y4 = vertices[13];
648 const PetscScalar z4 = vertices[14];
649 const PetscScalar x5 = vertices[15];
650 const PetscScalar y5 = vertices[16];
651 const PetscScalar z5 = vertices[17];
652 const PetscScalar x6 = vertices[18];
653 const PetscScalar y6 = vertices[19];
654 const PetscScalar z6 = vertices[20];
655 const PetscScalar x7 = vertices[21];
656 const PetscScalar y7 = vertices[22];
657 const PetscScalar z7 = vertices[23];
658 const PetscScalar f_1 = x1 - x0;
659 const PetscScalar g_1 = y1 - y0;
660 const PetscScalar h_1 = z1 - z0;
661 const PetscScalar f_3 = x3 - x0;
662 const PetscScalar g_3 = y3 - y0;
663 const PetscScalar h_3 = z3 - z0;
664 const PetscScalar f_4 = x4 - x0;
665 const PetscScalar g_4 = y4 - y0;
666 const PetscScalar h_4 = z4 - z0;
667 const PetscScalar f_01 = x2 - x1 - x3 + x0;
668 const PetscScalar g_01 = y2 - y1 - y3 + y0;
669 const PetscScalar h_01 = z2 - z1 - z3 + z0;
670 const PetscScalar f_12 = x7 - x3 - x4 + x0;
671 const PetscScalar g_12 = y7 - y3 - y4 + y0;
672 const PetscScalar h_12 = z7 - z3 - z4 + z0;
673 const PetscScalar f_02 = x5 - x1 - x4 + x0;
674 const PetscScalar g_02 = y5 - y1 - y4 + y0;
675 const PetscScalar h_02 = z5 - z1 - z4 + z0;
676 const PetscScalar f_012 = x6 - x0 + x1 - x2 + x3 + x4 - x5 - x7;
677 const PetscScalar g_012 = y6 - y0 + y1 - y2 + y3 + y4 - y5 - y7;
678 const PetscScalar h_012 = z6 - z0 + z1 - z2 + z3 + z4 - z5 - z7;
679 const PetscScalar *ref;
680 PetscScalar *real;
681 PetscErrorCode ierr;
682
683 PetscFunctionBegin;
684 ierr = VecGetArrayRead(Xref, &ref);CHKERRQ(ierr);
685 ierr = VecGetArray(Xreal, &real);CHKERRQ(ierr);
686 {
687 const PetscScalar p0 = ref[0];
688 const PetscScalar p1 = ref[1];
689 const PetscScalar p2 = ref[2];
690
691 real[0] = x0 + f_1*p0 + f_3*p1 + f_4*p2 + f_01*p0*p1 + f_12*p1*p2 + f_02*p0*p2 + f_012*p0*p1*p2;
692 real[1] = y0 + g_1*p0 + g_3*p1 + g_4*p2 + g_01*p0*p1 + g_01*p0*p1 + g_12*p1*p2 + g_02*p0*p2 + g_012*p0*p1*p2;
693 real[2] = z0 + h_1*p0 + h_3*p1 + h_4*p2 + h_01*p0*p1 + h_01*p0*p1 + h_12*p1*p2 + h_02*p0*p2 + h_012*p0*p1*p2;
694 }
695 ierr = PetscLogFlops(114);CHKERRQ(ierr);
696 ierr = VecRestoreArrayRead(Xref, &ref);CHKERRQ(ierr);
697 ierr = VecRestoreArray(Xreal, &real);CHKERRQ(ierr);
698 PetscFunctionReturn(0);
699 }
700
HexJacobian_Private(SNES snes,Vec Xref,Mat J,Mat M,void * ctx)701 PETSC_STATIC_INLINE PetscErrorCode HexJacobian_Private(SNES snes, Vec Xref, Mat J, Mat M, void *ctx)
702 {
703 const PetscScalar *vertices = (const PetscScalar*) ctx;
704 const PetscScalar x0 = vertices[0];
705 const PetscScalar y0 = vertices[1];
706 const PetscScalar z0 = vertices[2];
707 const PetscScalar x1 = vertices[9];
708 const PetscScalar y1 = vertices[10];
709 const PetscScalar z1 = vertices[11];
710 const PetscScalar x2 = vertices[6];
711 const PetscScalar y2 = vertices[7];
712 const PetscScalar z2 = vertices[8];
713 const PetscScalar x3 = vertices[3];
714 const PetscScalar y3 = vertices[4];
715 const PetscScalar z3 = vertices[5];
716 const PetscScalar x4 = vertices[12];
717 const PetscScalar y4 = vertices[13];
718 const PetscScalar z4 = vertices[14];
719 const PetscScalar x5 = vertices[15];
720 const PetscScalar y5 = vertices[16];
721 const PetscScalar z5 = vertices[17];
722 const PetscScalar x6 = vertices[18];
723 const PetscScalar y6 = vertices[19];
724 const PetscScalar z6 = vertices[20];
725 const PetscScalar x7 = vertices[21];
726 const PetscScalar y7 = vertices[22];
727 const PetscScalar z7 = vertices[23];
728 const PetscScalar f_xy = x2 - x1 - x3 + x0;
729 const PetscScalar g_xy = y2 - y1 - y3 + y0;
730 const PetscScalar h_xy = z2 - z1 - z3 + z0;
731 const PetscScalar f_yz = x7 - x3 - x4 + x0;
732 const PetscScalar g_yz = y7 - y3 - y4 + y0;
733 const PetscScalar h_yz = z7 - z3 - z4 + z0;
734 const PetscScalar f_xz = x5 - x1 - x4 + x0;
735 const PetscScalar g_xz = y5 - y1 - y4 + y0;
736 const PetscScalar h_xz = z5 - z1 - z4 + z0;
737 const PetscScalar f_xyz = x6 - x0 + x1 - x2 + x3 + x4 - x5 - x7;
738 const PetscScalar g_xyz = y6 - y0 + y1 - y2 + y3 + y4 - y5 - y7;
739 const PetscScalar h_xyz = z6 - z0 + z1 - z2 + z3 + z4 - z5 - z7;
740 const PetscScalar *ref;
741 PetscErrorCode ierr;
742
743 PetscFunctionBegin;
744 ierr = VecGetArrayRead(Xref, &ref);CHKERRQ(ierr);
745 {
746 const PetscScalar x = ref[0];
747 const PetscScalar y = ref[1];
748 const PetscScalar z = ref[2];
749 const PetscInt rows[3] = {0, 1, 2};
750 PetscScalar values[9];
751
752 values[0] = (x1 - x0 + f_xy*y + f_xz*z + f_xyz*y*z) / 2.0;
753 values[1] = (x3 - x0 + f_xy*x + f_yz*z + f_xyz*x*z) / 2.0;
754 values[2] = (x4 - x0 + f_yz*y + f_xz*x + f_xyz*x*y) / 2.0;
755 values[3] = (y1 - y0 + g_xy*y + g_xz*z + g_xyz*y*z) / 2.0;
756 values[4] = (y3 - y0 + g_xy*x + g_yz*z + g_xyz*x*z) / 2.0;
757 values[5] = (y4 - y0 + g_yz*y + g_xz*x + g_xyz*x*y) / 2.0;
758 values[6] = (z1 - z0 + h_xy*y + h_xz*z + h_xyz*y*z) / 2.0;
759 values[7] = (z3 - z0 + h_xy*x + h_yz*z + h_xyz*x*z) / 2.0;
760 values[8] = (z4 - z0 + h_yz*y + h_xz*x + h_xyz*x*y) / 2.0;
761
762 ierr = MatSetValues(J, 3, rows, 3, rows, values, INSERT_VALUES);CHKERRQ(ierr);
763 }
764 ierr = PetscLogFlops(152);CHKERRQ(ierr);
765 ierr = VecRestoreArrayRead(Xref, &ref);CHKERRQ(ierr);
766 ierr = MatAssemblyBegin(J, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
767 ierr = MatAssemblyEnd(J, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
768 PetscFunctionReturn(0);
769 }
770
DMInterpolate_Hex_Private(DMInterpolationInfo ctx,DM dm,Vec xLocal,Vec v)771 PETSC_STATIC_INLINE PetscErrorCode DMInterpolate_Hex_Private(DMInterpolationInfo ctx, DM dm, Vec xLocal, Vec v)
772 {
773 DM dmCoord;
774 SNES snes;
775 KSP ksp;
776 PC pc;
777 Vec coordsLocal, r, ref, real;
778 Mat J;
779 const PetscScalar *coords;
780 PetscScalar *a;
781 PetscInt p;
782 PetscErrorCode ierr;
783
784 PetscFunctionBegin;
785 ierr = DMGetCoordinatesLocal(dm, &coordsLocal);CHKERRQ(ierr);
786 ierr = DMGetCoordinateDM(dm, &dmCoord);CHKERRQ(ierr);
787 ierr = SNESCreate(PETSC_COMM_SELF, &snes);CHKERRQ(ierr);
788 ierr = SNESSetOptionsPrefix(snes, "hex_interp_");CHKERRQ(ierr);
789 ierr = VecCreate(PETSC_COMM_SELF, &r);CHKERRQ(ierr);
790 ierr = VecSetSizes(r, 3, 3);CHKERRQ(ierr);
791 ierr = VecSetType(r,dm->vectype);CHKERRQ(ierr);
792 ierr = VecDuplicate(r, &ref);CHKERRQ(ierr);
793 ierr = VecDuplicate(r, &real);CHKERRQ(ierr);
794 ierr = MatCreate(PETSC_COMM_SELF, &J);CHKERRQ(ierr);
795 ierr = MatSetSizes(J, 3, 3, 3, 3);CHKERRQ(ierr);
796 ierr = MatSetType(J, MATSEQDENSE);CHKERRQ(ierr);
797 ierr = MatSetUp(J);CHKERRQ(ierr);
798 ierr = SNESSetFunction(snes, r, HexMap_Private, NULL);CHKERRQ(ierr);
799 ierr = SNESSetJacobian(snes, J, J, HexJacobian_Private, NULL);CHKERRQ(ierr);
800 ierr = SNESGetKSP(snes, &ksp);CHKERRQ(ierr);
801 ierr = KSPGetPC(ksp, &pc);CHKERRQ(ierr);
802 ierr = PCSetType(pc, PCLU);CHKERRQ(ierr);
803 ierr = SNESSetFromOptions(snes);CHKERRQ(ierr);
804
805 ierr = VecGetArrayRead(ctx->coords, &coords);CHKERRQ(ierr);
806 ierr = VecGetArray(v, &a);CHKERRQ(ierr);
807 for (p = 0; p < ctx->n; ++p) {
808 PetscScalar *x = NULL, *vertices = NULL;
809 PetscScalar *xi;
810 PetscReal xir[3];
811 PetscInt c = ctx->cells[p], comp, coordSize, xSize;
812
813 /* Can make this do all points at once */
814 ierr = DMPlexVecGetClosure(dmCoord, NULL, coordsLocal, c, &coordSize, &vertices);CHKERRQ(ierr);
815 if (8*3 != coordSize) SETERRQ2(ctx->comm, PETSC_ERR_ARG_SIZ, "Invalid closure size %D should be %d", coordSize, 8*3);
816 ierr = DMPlexVecGetClosure(dm, NULL, xLocal, c, &xSize, &x);CHKERRQ(ierr);
817 if (8*ctx->dof != xSize) SETERRQ2(ctx->comm, PETSC_ERR_ARG_SIZ, "Invalid closure size %D should be %D", xSize, 8*ctx->dof);
818 ierr = SNESSetFunction(snes, NULL, NULL, (void*) vertices);CHKERRQ(ierr);
819 ierr = SNESSetJacobian(snes, NULL, NULL, NULL, (void*) vertices);CHKERRQ(ierr);
820 ierr = VecGetArray(real, &xi);CHKERRQ(ierr);
821 xi[0] = coords[p*ctx->dim+0];
822 xi[1] = coords[p*ctx->dim+1];
823 xi[2] = coords[p*ctx->dim+2];
824 ierr = VecRestoreArray(real, &xi);CHKERRQ(ierr);
825 ierr = SNESSolve(snes, real, ref);CHKERRQ(ierr);
826 ierr = VecGetArray(ref, &xi);CHKERRQ(ierr);
827 xir[0] = PetscRealPart(xi[0]);
828 xir[1] = PetscRealPart(xi[1]);
829 xir[2] = PetscRealPart(xi[2]);
830 for (comp = 0; comp < ctx->dof; ++comp) {
831 a[p*ctx->dof+comp] =
832 x[0*ctx->dof+comp]*(1-xir[0])*(1-xir[1])*(1-xir[2]) +
833 x[3*ctx->dof+comp]* xir[0]*(1-xir[1])*(1-xir[2]) +
834 x[2*ctx->dof+comp]* xir[0]* xir[1]*(1-xir[2]) +
835 x[1*ctx->dof+comp]*(1-xir[0])* xir[1]*(1-xir[2]) +
836 x[4*ctx->dof+comp]*(1-xir[0])*(1-xir[1])* xir[2] +
837 x[5*ctx->dof+comp]* xir[0]*(1-xir[1])* xir[2] +
838 x[6*ctx->dof+comp]* xir[0]* xir[1]* xir[2] +
839 x[7*ctx->dof+comp]*(1-xir[0])* xir[1]* xir[2];
840 }
841 ierr = VecRestoreArray(ref, &xi);CHKERRQ(ierr);
842 ierr = DMPlexVecRestoreClosure(dmCoord, NULL, coordsLocal, c, &coordSize, &vertices);CHKERRQ(ierr);
843 ierr = DMPlexVecRestoreClosure(dm, NULL, xLocal, c, &xSize, &x);CHKERRQ(ierr);
844 }
845 ierr = VecRestoreArray(v, &a);CHKERRQ(ierr);
846 ierr = VecRestoreArrayRead(ctx->coords, &coords);CHKERRQ(ierr);
847
848 ierr = SNESDestroy(&snes);CHKERRQ(ierr);
849 ierr = VecDestroy(&r);CHKERRQ(ierr);
850 ierr = VecDestroy(&ref);CHKERRQ(ierr);
851 ierr = VecDestroy(&real);CHKERRQ(ierr);
852 ierr = MatDestroy(&J);CHKERRQ(ierr);
853 PetscFunctionReturn(0);
854 }
855
856 /*@C
857 DMInterpolationEvaluate - Using the input from dm and x, calculates interpolated field values at the interpolation points.
858
859 Input Parameters:
860 + ctx - The DMInterpolationInfo context
861 . dm - The DM
862 - x - The local vector containing the field to be interpolated
863
864 Output Parameters:
865 . v - The vector containing the interpolated values
866
867 Note: A suitable v can be obtained using DMInterpolationGetVector().
868
869 Level: beginner
870
871 .seealso: DMInterpolationGetVector(), DMInterpolationAddPoints(), DMInterpolationCreate()
872 @*/
DMInterpolationEvaluate(DMInterpolationInfo ctx,DM dm,Vec x,Vec v)873 PetscErrorCode DMInterpolationEvaluate(DMInterpolationInfo ctx, DM dm, Vec x, Vec v)
874 {
875 PetscInt dim, coneSize, n;
876 PetscErrorCode ierr;
877
878 PetscFunctionBegin;
879 PetscValidHeaderSpecific(dm, DM_CLASSID, 2);
880 PetscValidHeaderSpecific(x, VEC_CLASSID, 3);
881 PetscValidHeaderSpecific(v, VEC_CLASSID, 4);
882 ierr = VecGetLocalSize(v, &n);CHKERRQ(ierr);
883 if (n != ctx->n*ctx->dof) SETERRQ2(ctx->comm, PETSC_ERR_ARG_SIZ, "Invalid input vector size %D should be %D", n, ctx->n*ctx->dof);
884 if (n) {
885 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr);
886 ierr = DMPlexGetConeSize(dm, ctx->cells[0], &coneSize);CHKERRQ(ierr);
887 if (dim == 2) {
888 if (coneSize == 3) {
889 ierr = DMInterpolate_Triangle_Private(ctx, dm, x, v);CHKERRQ(ierr);
890 } else if (coneSize == 4) {
891 ierr = DMInterpolate_Quad_Private(ctx, dm, x, v);CHKERRQ(ierr);
892 } else SETERRQ1(ctx->comm, PETSC_ERR_ARG_OUTOFRANGE, "Unsupported dimension %D for point interpolation", dim);
893 } else if (dim == 3) {
894 if (coneSize == 4) {
895 ierr = DMInterpolate_Tetrahedron_Private(ctx, dm, x, v);CHKERRQ(ierr);
896 } else {
897 ierr = DMInterpolate_Hex_Private(ctx, dm, x, v);CHKERRQ(ierr);
898 }
899 } else SETERRQ1(ctx->comm, PETSC_ERR_ARG_OUTOFRANGE, "Unsupported dimension %D for point interpolation", dim);
900 }
901 PetscFunctionReturn(0);
902 }
903
904 /*@C
905 DMInterpolationDestroy - Destroys a DMInterpolationInfo context
906
907 Collective on ctx
908
909 Input Parameter:
910 . ctx - the context
911
912 Level: beginner
913
914 .seealso: DMInterpolationEvaluate(), DMInterpolationAddPoints(), DMInterpolationCreate()
915 @*/
DMInterpolationDestroy(DMInterpolationInfo * ctx)916 PetscErrorCode DMInterpolationDestroy(DMInterpolationInfo *ctx)
917 {
918 PetscErrorCode ierr;
919
920 PetscFunctionBegin;
921 PetscValidPointer(ctx, 2);
922 ierr = VecDestroy(&(*ctx)->coords);CHKERRQ(ierr);
923 ierr = PetscFree((*ctx)->points);CHKERRQ(ierr);
924 ierr = PetscFree((*ctx)->cells);CHKERRQ(ierr);
925 ierr = PetscFree(*ctx);CHKERRQ(ierr);
926 *ctx = NULL;
927 PetscFunctionReturn(0);
928 }
929
930 /*@C
931 SNESMonitorFields - Monitors the residual for each field separately
932
933 Collective on SNES
934
935 Input Parameters:
936 + snes - the SNES context
937 . its - iteration number
938 . fgnorm - 2-norm of residual
939 - vf - PetscViewerAndFormat of type ASCII
940
941 Notes:
942 This routine prints the residual norm at each iteration.
943
944 Level: intermediate
945
946 .seealso: SNESMonitorSet(), SNESMonitorDefault()
947 @*/
SNESMonitorFields(SNES snes,PetscInt its,PetscReal fgnorm,PetscViewerAndFormat * vf)948 PetscErrorCode SNESMonitorFields(SNES snes, PetscInt its, PetscReal fgnorm, PetscViewerAndFormat *vf)
949 {
950 PetscViewer viewer = vf->viewer;
951 Vec res;
952 DM dm;
953 PetscSection s;
954 const PetscScalar *r;
955 PetscReal *lnorms, *norms;
956 PetscInt numFields, f, pStart, pEnd, p;
957 PetscErrorCode ierr;
958
959 PetscFunctionBegin;
960 PetscValidHeaderSpecific(viewer,PETSC_VIEWER_CLASSID,4);
961 ierr = SNESGetFunction(snes, &res, NULL, NULL);CHKERRQ(ierr);
962 ierr = SNESGetDM(snes, &dm);CHKERRQ(ierr);
963 ierr = DMGetLocalSection(dm, &s);CHKERRQ(ierr);
964 ierr = PetscSectionGetNumFields(s, &numFields);CHKERRQ(ierr);
965 ierr = PetscSectionGetChart(s, &pStart, &pEnd);CHKERRQ(ierr);
966 ierr = PetscCalloc2(numFields, &lnorms, numFields, &norms);CHKERRQ(ierr);
967 ierr = VecGetArrayRead(res, &r);CHKERRQ(ierr);
968 for (p = pStart; p < pEnd; ++p) {
969 for (f = 0; f < numFields; ++f) {
970 PetscInt fdof, foff, d;
971
972 ierr = PetscSectionGetFieldDof(s, p, f, &fdof);CHKERRQ(ierr);
973 ierr = PetscSectionGetFieldOffset(s, p, f, &foff);CHKERRQ(ierr);
974 for (d = 0; d < fdof; ++d) lnorms[f] += PetscRealPart(PetscSqr(r[foff+d]));
975 }
976 }
977 ierr = VecRestoreArrayRead(res, &r);CHKERRQ(ierr);
978 ierr = MPIU_Allreduce(lnorms, norms, numFields, MPIU_REAL, MPIU_SUM, PetscObjectComm((PetscObject) dm));CHKERRQ(ierr);
979 ierr = PetscViewerPushFormat(viewer,vf->format);CHKERRQ(ierr);
980 ierr = PetscViewerASCIIAddTab(viewer, ((PetscObject) snes)->tablevel);CHKERRQ(ierr);
981 ierr = PetscViewerASCIIPrintf(viewer, "%3D SNES Function norm %14.12e [", its, (double) fgnorm);CHKERRQ(ierr);
982 for (f = 0; f < numFields; ++f) {
983 if (f > 0) {ierr = PetscViewerASCIIPrintf(viewer, ", ");CHKERRQ(ierr);}
984 ierr = PetscViewerASCIIPrintf(viewer, "%14.12e", (double) PetscSqrtReal(norms[f]));CHKERRQ(ierr);
985 }
986 ierr = PetscViewerASCIIPrintf(viewer, "]\n");CHKERRQ(ierr);
987 ierr = PetscViewerASCIISubtractTab(viewer, ((PetscObject) snes)->tablevel);CHKERRQ(ierr);
988 ierr = PetscViewerPopFormat(viewer);CHKERRQ(ierr);
989 ierr = PetscFree2(lnorms, norms);CHKERRQ(ierr);
990 PetscFunctionReturn(0);
991 }
992
993 /********************* Residual Computation **************************/
994
995 /*@
996 DMPlexSNESComputeResidualFEM - Form the local residual F from the local input X using pointwise functions specified by the user
997
998 Input Parameters:
999 + dm - The mesh
1000 . X - Local solution
1001 - user - The user context
1002
1003 Output Parameter:
1004 . F - Local output vector
1005
1006 Level: developer
1007
1008 .seealso: DMPlexComputeJacobianAction()
1009 @*/
DMPlexSNESComputeResidualFEM(DM dm,Vec X,Vec F,void * user)1010 PetscErrorCode DMPlexSNESComputeResidualFEM(DM dm, Vec X, Vec F, void *user)
1011 {
1012 DM plex;
1013 IS allcellIS;
1014 PetscInt Nds, s, depth;
1015 PetscErrorCode ierr;
1016
1017 PetscFunctionBegin;
1018 ierr = DMGetNumDS(dm, &Nds);CHKERRQ(ierr);
1019 ierr = DMSNESConvertPlex(dm, &plex, PETSC_TRUE);CHKERRQ(ierr);
1020 ierr = DMPlexGetDepth(plex, &depth);CHKERRQ(ierr);
1021 ierr = DMGetStratumIS(plex, "dim", depth, &allcellIS);CHKERRQ(ierr);
1022 if (!allcellIS) {ierr = DMGetStratumIS(plex, "depth", depth, &allcellIS);CHKERRQ(ierr);}
1023 for (s = 0; s < Nds; ++s) {
1024 PetscDS ds;
1025 DMLabel label;
1026 IS cellIS;
1027
1028 ierr = DMGetRegionNumDS(dm, s, &label, NULL, &ds);CHKERRQ(ierr);
1029 if (!label) {
1030 ierr = PetscObjectReference((PetscObject) allcellIS);CHKERRQ(ierr);
1031 cellIS = allcellIS;
1032 } else {
1033 IS pointIS;
1034
1035 ierr = DMLabelGetStratumIS(label, 1, &pointIS);CHKERRQ(ierr);
1036 ierr = ISIntersect_Caching_Internal(allcellIS, pointIS, &cellIS);CHKERRQ(ierr);
1037 ierr = ISDestroy(&pointIS);CHKERRQ(ierr);
1038 }
1039 ierr = DMPlexComputeResidual_Internal(plex, cellIS, PETSC_MIN_REAL, X, NULL, 0.0, F, user);CHKERRQ(ierr);
1040 ierr = ISDestroy(&cellIS);CHKERRQ(ierr);
1041 }
1042 ierr = ISDestroy(&allcellIS);CHKERRQ(ierr);
1043 ierr = DMDestroy(&plex);CHKERRQ(ierr);
1044 PetscFunctionReturn(0);
1045 }
1046
1047 /*@
1048 DMPlexSNESComputeBoundaryFEM - Form the boundary values for the local input X
1049
1050 Input Parameters:
1051 + dm - The mesh
1052 - user - The user context
1053
1054 Output Parameter:
1055 . X - Local solution
1056
1057 Level: developer
1058
1059 .seealso: DMPlexComputeJacobianAction()
1060 @*/
DMPlexSNESComputeBoundaryFEM(DM dm,Vec X,void * user)1061 PetscErrorCode DMPlexSNESComputeBoundaryFEM(DM dm, Vec X, void *user)
1062 {
1063 DM plex;
1064 PetscErrorCode ierr;
1065
1066 PetscFunctionBegin;
1067 ierr = DMSNESConvertPlex(dm,&plex,PETSC_TRUE);CHKERRQ(ierr);
1068 ierr = DMPlexInsertBoundaryValues(plex, PETSC_TRUE, X, PETSC_MIN_REAL, NULL, NULL, NULL);CHKERRQ(ierr);
1069 ierr = DMDestroy(&plex);CHKERRQ(ierr);
1070 PetscFunctionReturn(0);
1071 }
1072
1073 /*@
1074 DMPlexComputeJacobianAction - Form the local portion of the Jacobian action Z = J(X) Y at the local solution X using pointwise functions specified by the user.
1075
1076 Input Parameters:
1077 + dm - The mesh
1078 . cellIS -
1079 . t - The time
1080 . X_tShift - The multiplier for the Jacobian with repsect to X_t
1081 . X - Local solution vector
1082 . X_t - Time-derivative of the local solution vector
1083 . Y - Local input vector
1084 - user - The user context
1085
1086 Output Parameter:
1087 . Z - Local output vector
1088
1089 Note:
1090 We form the residual one batch of elements at a time. This allows us to offload work onto an accelerator,
1091 like a GPU, or vectorize on a multicore machine.
1092
1093 Level: developer
1094
1095 .seealso: FormFunctionLocal()
1096 @*/
DMPlexComputeJacobianAction(DM dm,IS cellIS,PetscReal t,PetscReal X_tShift,Vec X,Vec X_t,Vec Y,Vec Z,void * user)1097 PetscErrorCode DMPlexComputeJacobianAction(DM dm, IS cellIS, PetscReal t, PetscReal X_tShift, Vec X, Vec X_t, Vec Y, Vec Z, void *user)
1098 {
1099 DM_Plex *mesh = (DM_Plex *) dm->data;
1100 const char *name = "Jacobian";
1101 DM dmAux, plex, plexAux = NULL;
1102 DMEnclosureType encAux;
1103 Vec A;
1104 PetscDS prob, probAux = NULL;
1105 PetscQuadrature quad;
1106 PetscSection section, globalSection, sectionAux;
1107 PetscScalar *elemMat, *elemMatD, *u, *u_t, *a = NULL, *y, *z;
1108 PetscInt Nf, fieldI, fieldJ;
1109 PetscInt totDim, totDimAux = 0;
1110 const PetscInt *cells;
1111 PetscInt cStart, cEnd, numCells, c;
1112 PetscBool hasDyn;
1113 DMField coordField;
1114 PetscErrorCode ierr;
1115
1116 PetscFunctionBegin;
1117 ierr = PetscLogEventBegin(DMPLEX_JacobianFEM,dm,0,0,0);CHKERRQ(ierr);
1118 ierr = DMSNESConvertPlex(dm, &plex, PETSC_TRUE);CHKERRQ(ierr);
1119 if (!cellIS) {
1120 PetscInt depth;
1121
1122 ierr = DMPlexGetDepth(plex, &depth);CHKERRQ(ierr);
1123 ierr = DMGetStratumIS(plex, "dim", depth, &cellIS);CHKERRQ(ierr);
1124 if (!cellIS) {ierr = DMGetStratumIS(plex, "depth", depth, &cellIS);CHKERRQ(ierr);}
1125 } else {
1126 ierr = PetscObjectReference((PetscObject) cellIS);CHKERRQ(ierr);
1127 }
1128 ierr = DMGetLocalSection(dm, §ion);CHKERRQ(ierr);
1129 ierr = DMGetGlobalSection(dm, &globalSection);CHKERRQ(ierr);
1130 ierr = ISGetLocalSize(cellIS, &numCells);CHKERRQ(ierr);
1131 ierr = ISGetPointRange(cellIS, &cStart, &cEnd, &cells);CHKERRQ(ierr);
1132 ierr = DMGetCellDS(dm, cells ? cells[cStart] : cStart, &prob);CHKERRQ(ierr);
1133 ierr = PetscDSGetNumFields(prob, &Nf);CHKERRQ(ierr);
1134 ierr = PetscDSGetTotalDimension(prob, &totDim);CHKERRQ(ierr);
1135 ierr = PetscDSHasDynamicJacobian(prob, &hasDyn);CHKERRQ(ierr);
1136 hasDyn = hasDyn && (X_tShift != 0.0) ? PETSC_TRUE : PETSC_FALSE;
1137 ierr = PetscObjectQuery((PetscObject) dm, "dmAux", (PetscObject *) &dmAux);CHKERRQ(ierr);
1138 ierr = PetscObjectQuery((PetscObject) dm, "A", (PetscObject *) &A);CHKERRQ(ierr);
1139 if (dmAux) {
1140 ierr = DMGetEnclosureRelation(dmAux, dm, &encAux);CHKERRQ(ierr);
1141 ierr = DMConvert(dmAux, DMPLEX, &plexAux);CHKERRQ(ierr);
1142 ierr = DMGetLocalSection(plexAux, §ionAux);CHKERRQ(ierr);
1143 ierr = DMGetDS(dmAux, &probAux);CHKERRQ(ierr);
1144 ierr = PetscDSGetTotalDimension(probAux, &totDimAux);CHKERRQ(ierr);
1145 }
1146 ierr = VecSet(Z, 0.0);CHKERRQ(ierr);
1147 ierr = PetscMalloc6(numCells*totDim,&u,X_t ? numCells*totDim : 0,&u_t,numCells*totDim*totDim,&elemMat,hasDyn ? numCells*totDim*totDim : 0, &elemMatD,numCells*totDim,&y,totDim,&z);CHKERRQ(ierr);
1148 if (dmAux) {ierr = PetscMalloc1(numCells*totDimAux, &a);CHKERRQ(ierr);}
1149 ierr = DMGetCoordinateField(dm, &coordField);CHKERRQ(ierr);
1150 for (c = cStart; c < cEnd; ++c) {
1151 const PetscInt cell = cells ? cells[c] : c;
1152 const PetscInt cind = c - cStart;
1153 PetscScalar *x = NULL, *x_t = NULL;
1154 PetscInt i;
1155
1156 ierr = DMPlexVecGetClosure(dm, section, X, cell, NULL, &x);CHKERRQ(ierr);
1157 for (i = 0; i < totDim; ++i) u[cind*totDim+i] = x[i];
1158 ierr = DMPlexVecRestoreClosure(dm, section, X, cell, NULL, &x);CHKERRQ(ierr);
1159 if (X_t) {
1160 ierr = DMPlexVecGetClosure(dm, section, X_t, cell, NULL, &x_t);CHKERRQ(ierr);
1161 for (i = 0; i < totDim; ++i) u_t[cind*totDim+i] = x_t[i];
1162 ierr = DMPlexVecRestoreClosure(dm, section, X_t, cell, NULL, &x_t);CHKERRQ(ierr);
1163 }
1164 if (dmAux) {
1165 PetscInt subcell;
1166 ierr = DMGetEnclosurePoint(dmAux, dm, encAux, cell, &subcell);CHKERRQ(ierr);
1167 ierr = DMPlexVecGetClosure(plexAux, sectionAux, A, subcell, NULL, &x);CHKERRQ(ierr);
1168 for (i = 0; i < totDimAux; ++i) a[cind*totDimAux+i] = x[i];
1169 ierr = DMPlexVecRestoreClosure(plexAux, sectionAux, A, subcell, NULL, &x);CHKERRQ(ierr);
1170 }
1171 ierr = DMPlexVecGetClosure(dm, section, Y, cell, NULL, &x);CHKERRQ(ierr);
1172 for (i = 0; i < totDim; ++i) y[cind*totDim+i] = x[i];
1173 ierr = DMPlexVecRestoreClosure(dm, section, Y, cell, NULL, &x);CHKERRQ(ierr);
1174 }
1175 ierr = PetscArrayzero(elemMat, numCells*totDim*totDim);CHKERRQ(ierr);
1176 if (hasDyn) {ierr = PetscArrayzero(elemMatD, numCells*totDim*totDim);CHKERRQ(ierr);}
1177 for (fieldI = 0; fieldI < Nf; ++fieldI) {
1178 PetscFE fe;
1179 PetscInt Nb;
1180 /* Conforming batches */
1181 PetscInt numChunks, numBatches, numBlocks, Ne, blockSize, batchSize;
1182 /* Remainder */
1183 PetscInt Nr, offset, Nq;
1184 PetscQuadrature qGeom = NULL;
1185 PetscInt maxDegree;
1186 PetscFEGeom *cgeomFEM, *chunkGeom = NULL, *remGeom = NULL;
1187
1188 ierr = PetscDSGetDiscretization(prob, fieldI, (PetscObject *) &fe);CHKERRQ(ierr);
1189 ierr = PetscFEGetQuadrature(fe, &quad);CHKERRQ(ierr);
1190 ierr = PetscFEGetDimension(fe, &Nb);CHKERRQ(ierr);
1191 ierr = PetscFEGetTileSizes(fe, NULL, &numBlocks, NULL, &numBatches);CHKERRQ(ierr);
1192 ierr = DMFieldGetDegree(coordField,cellIS,NULL,&maxDegree);CHKERRQ(ierr);
1193 if (maxDegree <= 1) {ierr = DMFieldCreateDefaultQuadrature(coordField,cellIS,&qGeom);CHKERRQ(ierr);}
1194 if (!qGeom) {
1195 ierr = PetscFEGetQuadrature(fe,&qGeom);CHKERRQ(ierr);
1196 ierr = PetscObjectReference((PetscObject)qGeom);CHKERRQ(ierr);
1197 }
1198 ierr = PetscQuadratureGetData(qGeom, NULL, NULL, &Nq, NULL, NULL);CHKERRQ(ierr);
1199 ierr = DMSNESGetFEGeom(coordField,cellIS,qGeom,PETSC_FALSE,&cgeomFEM);CHKERRQ(ierr);
1200 blockSize = Nb;
1201 batchSize = numBlocks * blockSize;
1202 ierr = PetscFESetTileSizes(fe, blockSize, numBlocks, batchSize, numBatches);CHKERRQ(ierr);
1203 numChunks = numCells / (numBatches*batchSize);
1204 Ne = numChunks*numBatches*batchSize;
1205 Nr = numCells % (numBatches*batchSize);
1206 offset = numCells - Nr;
1207 ierr = PetscFEGeomGetChunk(cgeomFEM,0,offset,&chunkGeom);CHKERRQ(ierr);
1208 ierr = PetscFEGeomGetChunk(cgeomFEM,offset,numCells,&remGeom);CHKERRQ(ierr);
1209 for (fieldJ = 0; fieldJ < Nf; ++fieldJ) {
1210 ierr = PetscFEIntegrateJacobian(prob, PETSCFE_JACOBIAN, fieldI, fieldJ, Ne, chunkGeom, u, u_t, probAux, a, t, X_tShift, elemMat);CHKERRQ(ierr);
1211 ierr = PetscFEIntegrateJacobian(prob, PETSCFE_JACOBIAN, fieldI, fieldJ, Nr, remGeom, &u[offset*totDim], u_t ? &u_t[offset*totDim] : NULL, probAux, &a[offset*totDimAux], t, X_tShift, &elemMat[offset*totDim*totDim]);CHKERRQ(ierr);
1212 if (hasDyn) {
1213 ierr = PetscFEIntegrateJacobian(prob, PETSCFE_JACOBIAN_DYN, fieldI, fieldJ, Ne, chunkGeom, u, u_t, probAux, a, t, X_tShift, elemMatD);CHKERRQ(ierr);
1214 ierr = PetscFEIntegrateJacobian(prob, PETSCFE_JACOBIAN_DYN, fieldI, fieldJ, Nr, remGeom, &u[offset*totDim], u_t ? &u_t[offset*totDim] : NULL, probAux, &a[offset*totDimAux], t, X_tShift, &elemMatD[offset*totDim*totDim]);CHKERRQ(ierr);
1215 }
1216 }
1217 ierr = PetscFEGeomRestoreChunk(cgeomFEM,offset,numCells,&remGeom);CHKERRQ(ierr);
1218 ierr = PetscFEGeomRestoreChunk(cgeomFEM,0,offset,&chunkGeom);CHKERRQ(ierr);
1219 ierr = DMSNESRestoreFEGeom(coordField,cellIS,qGeom,PETSC_FALSE,&cgeomFEM);CHKERRQ(ierr);
1220 ierr = PetscQuadratureDestroy(&qGeom);CHKERRQ(ierr);
1221 }
1222 if (hasDyn) {
1223 for (c = 0; c < numCells*totDim*totDim; ++c) elemMat[c] += X_tShift*elemMatD[c];
1224 }
1225 for (c = cStart; c < cEnd; ++c) {
1226 const PetscInt cell = cells ? cells[c] : c;
1227 const PetscInt cind = c - cStart;
1228 const PetscBLASInt M = totDim, one = 1;
1229 const PetscScalar a = 1.0, b = 0.0;
1230
1231 PetscStackCallBLAS("BLASgemv", BLASgemv_("N", &M, &M, &a, &elemMat[cind*totDim*totDim], &M, &y[cind*totDim], &one, &b, z, &one));
1232 if (mesh->printFEM > 1) {
1233 ierr = DMPrintCellMatrix(c, name, totDim, totDim, &elemMat[cind*totDim*totDim]);CHKERRQ(ierr);
1234 ierr = DMPrintCellVector(c, "Y", totDim, &y[cind*totDim]);CHKERRQ(ierr);
1235 ierr = DMPrintCellVector(c, "Z", totDim, z);CHKERRQ(ierr);
1236 }
1237 ierr = DMPlexVecSetClosure(dm, section, Z, cell, z, ADD_VALUES);CHKERRQ(ierr);
1238 }
1239 ierr = PetscFree6(u,u_t,elemMat,elemMatD,y,z);CHKERRQ(ierr);
1240 if (mesh->printFEM) {
1241 ierr = PetscPrintf(PetscObjectComm((PetscObject)Z), "Z:\n");CHKERRQ(ierr);
1242 ierr = VecView(Z, NULL);CHKERRQ(ierr);
1243 }
1244 ierr = PetscFree(a);CHKERRQ(ierr);
1245 ierr = ISDestroy(&cellIS);CHKERRQ(ierr);
1246 ierr = DMDestroy(&plexAux);CHKERRQ(ierr);
1247 ierr = DMDestroy(&plex);CHKERRQ(ierr);
1248 ierr = PetscLogEventEnd(DMPLEX_JacobianFEM,dm,0,0,0);CHKERRQ(ierr);
1249 PetscFunctionReturn(0);
1250 }
1251
1252 /*@
1253 DMPlexSNESComputeJacobianFEM - Form the local portion of the Jacobian matrix J at the local solution X using pointwise functions specified by the user.
1254
1255 Input Parameters:
1256 + dm - The mesh
1257 . X - Local input vector
1258 - user - The user context
1259
1260 Output Parameter:
1261 . Jac - Jacobian matrix
1262
1263 Note:
1264 We form the residual one batch of elements at a time. This allows us to offload work onto an accelerator,
1265 like a GPU, or vectorize on a multicore machine.
1266
1267 Level: developer
1268
1269 .seealso: FormFunctionLocal()
1270 @*/
DMPlexSNESComputeJacobianFEM(DM dm,Vec X,Mat Jac,Mat JacP,void * user)1271 PetscErrorCode DMPlexSNESComputeJacobianFEM(DM dm, Vec X, Mat Jac, Mat JacP,void *user)
1272 {
1273 DM plex;
1274 IS allcellIS;
1275 PetscBool hasJac, hasPrec;
1276 PetscInt Nds, s, depth;
1277 PetscErrorCode ierr;
1278
1279 PetscFunctionBegin;
1280 ierr = DMGetNumDS(dm, &Nds);CHKERRQ(ierr);
1281 ierr = DMSNESConvertPlex(dm, &plex, PETSC_TRUE);CHKERRQ(ierr);
1282 ierr = DMPlexGetDepth(plex, &depth);CHKERRQ(ierr);
1283 ierr = DMGetStratumIS(plex, "dim", depth, &allcellIS);CHKERRQ(ierr);
1284 if (!allcellIS) {ierr = DMGetStratumIS(plex, "depth", depth, &allcellIS);CHKERRQ(ierr);}
1285 for (s = 0; s < Nds; ++s) {
1286 PetscDS ds;
1287 DMLabel label;
1288 IS cellIS;
1289
1290 ierr = DMGetRegionNumDS(dm, s, &label, NULL, &ds);CHKERRQ(ierr);
1291 if (!label) {
1292 ierr = PetscObjectReference((PetscObject) allcellIS);CHKERRQ(ierr);
1293 cellIS = allcellIS;
1294 } else {
1295 IS pointIS;
1296
1297 ierr = DMLabelGetStratumIS(label, 1, &pointIS);CHKERRQ(ierr);
1298 ierr = ISIntersect_Caching_Internal(allcellIS, pointIS, &cellIS);CHKERRQ(ierr);
1299 ierr = ISDestroy(&pointIS);CHKERRQ(ierr);
1300 }
1301 if (!s) {
1302 ierr = PetscDSHasJacobian(ds, &hasJac);CHKERRQ(ierr);
1303 ierr = PetscDSHasJacobianPreconditioner(ds, &hasPrec);CHKERRQ(ierr);
1304 if (hasJac && hasPrec) {ierr = MatZeroEntries(Jac);CHKERRQ(ierr);}
1305 ierr = MatZeroEntries(JacP);CHKERRQ(ierr);
1306 }
1307 ierr = DMPlexComputeJacobian_Internal(plex, cellIS, 0.0, 0.0, X, NULL, Jac, JacP, user);CHKERRQ(ierr);
1308 ierr = ISDestroy(&cellIS);CHKERRQ(ierr);
1309 }
1310 ierr = ISDestroy(&allcellIS);CHKERRQ(ierr);
1311 ierr = DMDestroy(&plex);CHKERRQ(ierr);
1312 PetscFunctionReturn(0);
1313 }
1314
1315 /*
1316 MatComputeNeumannOverlap - Computes an unassembled (Neumann) local overlapping Mat in nonlinear context.
1317
1318 Input Parameters:
1319 + X - SNES linearization point
1320 . ovl - index set of overlapping subdomains
1321
1322 Output Parameter:
1323 . J - unassembled (Neumann) local matrix
1324
1325 Level: intermediate
1326
1327 .seealso: DMCreateNeumannOverlap(), MATIS, PCHPDDMSetAuxiliaryMat()
1328 */
MatComputeNeumannOverlap_Plex(Mat J,PetscReal t,Vec X,Vec X_t,PetscReal s,IS ovl,void * ctx)1329 static PetscErrorCode MatComputeNeumannOverlap_Plex(Mat J, PetscReal t, Vec X, Vec X_t, PetscReal s, IS ovl, void *ctx)
1330 {
1331 SNES snes;
1332 Mat pJ;
1333 DM ovldm,origdm;
1334 DMSNES sdm;
1335 PetscErrorCode (*bfun)(DM,Vec,void*);
1336 PetscErrorCode (*jfun)(DM,Vec,Mat,Mat,void*);
1337 void *bctx,*jctx;
1338 PetscErrorCode ierr;
1339
1340 PetscFunctionBegin;
1341 ierr = PetscObjectQuery((PetscObject)ovl,"_DM_Overlap_HPDDM_MATIS",(PetscObject*)&pJ);CHKERRQ(ierr);
1342 if (!pJ) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Missing overlapping Mat");CHKERRQ(ierr);
1343 ierr = PetscObjectQuery((PetscObject)ovl,"_DM_Original_HPDDM",(PetscObject*)&origdm);CHKERRQ(ierr);
1344 if (!origdm) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Missing original DM");CHKERRQ(ierr);
1345 ierr = MatGetDM(pJ,&ovldm);CHKERRQ(ierr);
1346 ierr = DMSNESGetBoundaryLocal(origdm,&bfun,&bctx);CHKERRQ(ierr);
1347 ierr = DMSNESSetBoundaryLocal(ovldm,bfun,bctx);CHKERRQ(ierr);
1348 ierr = DMSNESGetJacobianLocal(origdm,&jfun,&jctx);CHKERRQ(ierr);
1349 ierr = DMSNESSetJacobianLocal(ovldm,jfun,jctx);CHKERRQ(ierr);
1350 ierr = PetscObjectQuery((PetscObject)ovl,"_DM_Overlap_HPDDM_SNES",(PetscObject*)&snes);CHKERRQ(ierr);
1351 if (!snes) {
1352 ierr = SNESCreate(PetscObjectComm((PetscObject)ovl),&snes);CHKERRQ(ierr);
1353 ierr = SNESSetDM(snes,ovldm);CHKERRQ(ierr);
1354 ierr = PetscObjectCompose((PetscObject)ovl,"_DM_Overlap_HPDDM_SNES",(PetscObject)snes);CHKERRQ(ierr);
1355 ierr = PetscObjectDereference((PetscObject)snes);CHKERRQ(ierr);
1356 }
1357 ierr = DMGetDMSNES(ovldm,&sdm);CHKERRQ(ierr);
1358 ierr = VecLockReadPush(X);CHKERRQ(ierr);
1359 PetscStackPush("SNES user Jacobian function");
1360 ierr = (*sdm->ops->computejacobian)(snes,X,pJ,pJ,sdm->jacobianctx);CHKERRQ(ierr);
1361 PetscStackPop;
1362 ierr = VecLockReadPop(X);CHKERRQ(ierr);
1363 /* this is a no-hop, just in case we decide to change the placeholder for the local Neumann matrix */
1364 {
1365 Mat locpJ;
1366
1367 ierr = MatISGetLocalMat(pJ,&locpJ);CHKERRQ(ierr);
1368 ierr = MatCopy(locpJ,J,SAME_NONZERO_PATTERN);CHKERRQ(ierr);
1369 }
1370 PetscFunctionReturn(0);
1371 }
1372
1373 /*@
1374 DMPlexSetSNESLocalFEM - Use DMPlex's internal FEM routines to compute SNES boundary values, residual, and Jacobian.
1375
1376 Input Parameters:
1377 + dm - The DM object
1378 . boundaryctx - the user context that will be passed to pointwise evaluation of boundary values (see PetscDSAddBoundary())
1379 . residualctx - the user context that will be passed to pointwise evaluation of finite element residual computations (see PetscDSSetResidual())
1380 - jacobianctx - the user context that will be passed to pointwise evaluation of finite element Jacobian construction (see PetscDSSetJacobian())
1381
1382 Level: developer
1383 @*/
DMPlexSetSNESLocalFEM(DM dm,void * boundaryctx,void * residualctx,void * jacobianctx)1384 PetscErrorCode DMPlexSetSNESLocalFEM(DM dm, void *boundaryctx, void *residualctx, void *jacobianctx)
1385 {
1386 PetscErrorCode ierr;
1387
1388 PetscFunctionBegin;
1389 ierr = DMSNESSetBoundaryLocal(dm,DMPlexSNESComputeBoundaryFEM,boundaryctx);CHKERRQ(ierr);
1390 ierr = DMSNESSetFunctionLocal(dm,DMPlexSNESComputeResidualFEM,residualctx);CHKERRQ(ierr);
1391 ierr = DMSNESSetJacobianLocal(dm,DMPlexSNESComputeJacobianFEM,jacobianctx);CHKERRQ(ierr);
1392 ierr = PetscObjectComposeFunction((PetscObject)dm,"MatComputeNeumannOverlap_C",MatComputeNeumannOverlap_Plex);CHKERRQ(ierr);
1393 PetscFunctionReturn(0);
1394 }
1395
1396 /*@C
1397 DMSNESCheckDiscretization - Check the discretization error of the exact solution
1398
1399 Input Parameters:
1400 + snes - the SNES object
1401 . dm - the DM
1402 . t - the time
1403 . u - a DM vector
1404 - tol - A tolerance for the check, or -1 to print the results instead
1405
1406 Output Parameters:
1407 . error - An array which holds the discretization error in each field, or NULL
1408
1409 Note: The user must call PetscDSSetExactSolution() beforehand
1410
1411 Level: developer
1412
1413 .seealso: DNSNESCheckFromOptions(), DMSNESCheckResidual(), DMSNESCheckJacobian(), PetscDSSetExactSolution()
1414 @*/
DMSNESCheckDiscretization(SNES snes,DM dm,PetscReal t,Vec u,PetscReal tol,PetscReal error[])1415 PetscErrorCode DMSNESCheckDiscretization(SNES snes, DM dm, PetscReal t, Vec u, PetscReal tol, PetscReal error[])
1416 {
1417 PetscErrorCode (**exacts)(PetscInt, PetscReal, const PetscReal x[], PetscInt, PetscScalar *u, void *ctx);
1418 void **ectxs;
1419 PetscReal *err;
1420 MPI_Comm comm;
1421 PetscInt Nf, f;
1422 PetscErrorCode ierr;
1423
1424 PetscFunctionBegin;
1425 PetscValidHeaderSpecific(snes, SNES_CLASSID, 1);
1426 PetscValidHeaderSpecific(dm, DM_CLASSID, 2);
1427 PetscValidHeaderSpecific(u, VEC_CLASSID, 3);
1428 if (error) PetscValidRealPointer(error, 6);
1429
1430 ierr = DMComputeExactSolution(dm, t, u, NULL);CHKERRQ(ierr);
1431 ierr = VecViewFromOptions(u, NULL, "-vec_view");CHKERRQ(ierr);
1432
1433 ierr = PetscObjectGetComm((PetscObject) snes, &comm);CHKERRQ(ierr);
1434 ierr = DMGetNumFields(dm, &Nf);CHKERRQ(ierr);
1435 ierr = PetscCalloc3(Nf, &exacts, Nf, &ectxs, PetscMax(1, Nf), &err);CHKERRQ(ierr);
1436 {
1437 PetscInt Nds, s;
1438
1439 ierr = DMGetNumDS(dm, &Nds);CHKERRQ(ierr);
1440 for (s = 0; s < Nds; ++s) {
1441 PetscDS ds;
1442 DMLabel label;
1443 IS fieldIS;
1444 const PetscInt *fields;
1445 PetscInt dsNf, f;
1446
1447 ierr = DMGetRegionNumDS(dm, s, &label, &fieldIS, &ds);CHKERRQ(ierr);
1448 ierr = PetscDSGetNumFields(ds, &dsNf);CHKERRQ(ierr);
1449 ierr = ISGetIndices(fieldIS, &fields);CHKERRQ(ierr);
1450 for (f = 0; f < dsNf; ++f) {
1451 const PetscInt field = fields[f];
1452 ierr = PetscDSGetExactSolution(ds, field, &exacts[field], &ectxs[field]);CHKERRQ(ierr);
1453 }
1454 ierr = ISRestoreIndices(fieldIS, &fields);CHKERRQ(ierr);
1455 }
1456 }
1457 if (Nf > 1) {
1458 ierr = DMComputeL2FieldDiff(dm, t, exacts, ectxs, u, err);CHKERRQ(ierr);
1459 if (tol >= 0.0) {
1460 for (f = 0; f < Nf; ++f) {
1461 if (err[f] > tol) SETERRQ3(comm, PETSC_ERR_ARG_WRONG, "L_2 Error %g for field %D exceeds tolerance %g", (double) err[f], f, (double) tol);
1462 }
1463 } else if (error) {
1464 for (f = 0; f < Nf; ++f) error[f] = err[f];
1465 } else {
1466 ierr = PetscPrintf(comm, "L_2 Error: [");CHKERRQ(ierr);
1467 for (f = 0; f < Nf; ++f) {
1468 if (f) {ierr = PetscPrintf(comm, ", ");CHKERRQ(ierr);}
1469 ierr = PetscPrintf(comm, "%g", (double)err[f]);CHKERRQ(ierr);
1470 }
1471 ierr = PetscPrintf(comm, "]\n");CHKERRQ(ierr);
1472 }
1473 } else {
1474 ierr = DMComputeL2Diff(dm, t, exacts, ectxs, u, &err[0]);CHKERRQ(ierr);
1475 if (tol >= 0.0) {
1476 if (err[0] > tol) SETERRQ2(comm, PETSC_ERR_ARG_WRONG, "L_2 Error %g exceeds tolerance %g", (double) err[0], (double) tol);
1477 } else if (error) {
1478 error[0] = err[0];
1479 } else {
1480 ierr = PetscPrintf(comm, "L_2 Error: %g\n", (double) err[0]);CHKERRQ(ierr);
1481 }
1482 }
1483 ierr = PetscFree3(exacts, ectxs, err);CHKERRQ(ierr);
1484 PetscFunctionReturn(0);
1485 }
1486
1487 /*@C
1488 DMSNESCheckResidual - Check the residual of the exact solution
1489
1490 Input Parameters:
1491 + snes - the SNES object
1492 . dm - the DM
1493 . u - a DM vector
1494 - tol - A tolerance for the check, or -1 to print the results instead
1495
1496 Output Parameters:
1497 . residual - The residual norm of the exact solution, or NULL
1498
1499 Level: developer
1500
1501 .seealso: DNSNESCheckFromOptions(), DMSNESCheckDiscretization(), DMSNESCheckJacobian()
1502 @*/
DMSNESCheckResidual(SNES snes,DM dm,Vec u,PetscReal tol,PetscReal * residual)1503 PetscErrorCode DMSNESCheckResidual(SNES snes, DM dm, Vec u, PetscReal tol, PetscReal *residual)
1504 {
1505 MPI_Comm comm;
1506 Vec r;
1507 PetscReal res;
1508 PetscErrorCode ierr;
1509
1510 PetscFunctionBegin;
1511 PetscValidHeaderSpecific(snes, SNES_CLASSID, 1);
1512 PetscValidHeaderSpecific(dm, DM_CLASSID, 2);
1513 PetscValidHeaderSpecific(u, VEC_CLASSID, 3);
1514 if (residual) PetscValidRealPointer(residual, 5);
1515 ierr = PetscObjectGetComm((PetscObject) snes, &comm);CHKERRQ(ierr);
1516 ierr = DMComputeExactSolution(dm, 0.0, u, NULL);CHKERRQ(ierr);
1517 ierr = VecDuplicate(u, &r);CHKERRQ(ierr);
1518 ierr = SNESComputeFunction(snes, u, r);CHKERRQ(ierr);
1519 ierr = VecNorm(r, NORM_2, &res);CHKERRQ(ierr);
1520 if (tol >= 0.0) {
1521 if (res > tol) SETERRQ2(comm, PETSC_ERR_ARG_WRONG, "L_2 Residual %g exceeds tolerance %g", (double) res, (double) tol);
1522 } else if (residual) {
1523 *residual = res;
1524 } else {
1525 ierr = PetscPrintf(comm, "L_2 Residual: %g\n", (double)res);CHKERRQ(ierr);
1526 ierr = VecChop(r, 1.0e-10);CHKERRQ(ierr);
1527 ierr = PetscObjectSetName((PetscObject) r, "Initial Residual");CHKERRQ(ierr);
1528 ierr = PetscObjectSetOptionsPrefix((PetscObject)r,"res_");CHKERRQ(ierr);
1529 ierr = VecViewFromOptions(r, NULL, "-vec_view");CHKERRQ(ierr);
1530 }
1531 ierr = VecDestroy(&r);CHKERRQ(ierr);
1532 PetscFunctionReturn(0);
1533 }
1534
1535 /*@C
1536 DMSNESCheckJacobian - Check the Jacobian of the exact solution against the residual using the Taylor Test
1537
1538 Input Parameters:
1539 + snes - the SNES object
1540 . dm - the DM
1541 . u - a DM vector
1542 - tol - A tolerance for the check, or -1 to print the results instead
1543
1544 Output Parameters:
1545 + isLinear - Flag indicaing that the function looks linear, or NULL
1546 - convRate - The rate of convergence of the linear model, or NULL
1547
1548 Level: developer
1549
1550 .seealso: DNSNESCheckFromOptions(), DMSNESCheckDiscretization(), DMSNESCheckResidual()
1551 @*/
DMSNESCheckJacobian(SNES snes,DM dm,Vec u,PetscReal tol,PetscBool * isLinear,PetscReal * convRate)1552 PetscErrorCode DMSNESCheckJacobian(SNES snes, DM dm, Vec u, PetscReal tol, PetscBool *isLinear, PetscReal *convRate)
1553 {
1554 MPI_Comm comm;
1555 PetscDS ds;
1556 Mat J, M;
1557 MatNullSpace nullspace;
1558 PetscReal slope, intercept;
1559 PetscBool hasJac, hasPrec, isLin = PETSC_FALSE;
1560 PetscErrorCode ierr;
1561
1562 PetscFunctionBegin;
1563 PetscValidHeaderSpecific(snes, SNES_CLASSID, 1);
1564 PetscValidHeaderSpecific(dm, DM_CLASSID, 2);
1565 PetscValidHeaderSpecific(u, VEC_CLASSID, 3);
1566 if (isLinear) PetscValidBoolPointer(isLinear, 5);
1567 if (convRate) PetscValidRealPointer(convRate, 5);
1568 ierr = PetscObjectGetComm((PetscObject) snes, &comm);CHKERRQ(ierr);
1569 ierr = DMComputeExactSolution(dm, 0.0, u, NULL);CHKERRQ(ierr);
1570 /* Create and view matrices */
1571 ierr = DMCreateMatrix(dm, &J);CHKERRQ(ierr);
1572 ierr = DMGetDS(dm, &ds);CHKERRQ(ierr);
1573 ierr = PetscDSHasJacobian(ds, &hasJac);CHKERRQ(ierr);
1574 ierr = PetscDSHasJacobianPreconditioner(ds, &hasPrec);CHKERRQ(ierr);
1575 if (hasJac && hasPrec) {
1576 ierr = DMCreateMatrix(dm, &M);CHKERRQ(ierr);
1577 ierr = SNESComputeJacobian(snes, u, J, M);CHKERRQ(ierr);
1578 ierr = PetscObjectSetName((PetscObject) M, "Preconditioning Matrix");CHKERRQ(ierr);
1579 ierr = PetscObjectSetOptionsPrefix((PetscObject) M, "jacpre_");CHKERRQ(ierr);
1580 ierr = MatViewFromOptions(M, NULL, "-mat_view");CHKERRQ(ierr);
1581 ierr = MatDestroy(&M);CHKERRQ(ierr);
1582 } else {
1583 ierr = SNESComputeJacobian(snes, u, J, J);CHKERRQ(ierr);
1584 }
1585 ierr = PetscObjectSetName((PetscObject) J, "Jacobian");CHKERRQ(ierr);
1586 ierr = PetscObjectSetOptionsPrefix((PetscObject) J, "jac_");CHKERRQ(ierr);
1587 ierr = MatViewFromOptions(J, NULL, "-mat_view");CHKERRQ(ierr);
1588 /* Check nullspace */
1589 ierr = MatGetNullSpace(J, &nullspace);CHKERRQ(ierr);
1590 if (nullspace) {
1591 PetscBool isNull;
1592 ierr = MatNullSpaceTest(nullspace, J, &isNull);CHKERRQ(ierr);
1593 if (!isNull) SETERRQ(comm, PETSC_ERR_PLIB, "The null space calculated for the system operator is invalid.");
1594 }
1595 ierr = MatNullSpaceDestroy(&nullspace);CHKERRQ(ierr);
1596 /* Taylor test */
1597 {
1598 PetscRandom rand;
1599 Vec du, uhat, r, rhat, df;
1600 PetscReal h;
1601 PetscReal *es, *hs, *errors;
1602 PetscReal hMax = 1.0, hMin = 1e-6, hMult = 0.1;
1603 PetscInt Nv, v;
1604
1605 /* Choose a perturbation direction */
1606 ierr = PetscRandomCreate(comm, &rand);CHKERRQ(ierr);
1607 ierr = VecDuplicate(u, &du);CHKERRQ(ierr);
1608 ierr = VecSetRandom(du, rand);CHKERRQ(ierr);
1609 ierr = PetscRandomDestroy(&rand);CHKERRQ(ierr);
1610 ierr = VecDuplicate(u, &df);CHKERRQ(ierr);
1611 ierr = MatMult(J, du, df);CHKERRQ(ierr);
1612 /* Evaluate residual at u, F(u), save in vector r */
1613 ierr = VecDuplicate(u, &r);CHKERRQ(ierr);
1614 ierr = SNESComputeFunction(snes, u, r);CHKERRQ(ierr);
1615 /* Look at the convergence of our Taylor approximation as we approach u */
1616 for (h = hMax, Nv = 0; h >= hMin; h *= hMult, ++Nv);
1617 ierr = PetscCalloc3(Nv, &es, Nv, &hs, Nv, &errors);CHKERRQ(ierr);
1618 ierr = VecDuplicate(u, &uhat);CHKERRQ(ierr);
1619 ierr = VecDuplicate(u, &rhat);CHKERRQ(ierr);
1620 for (h = hMax, Nv = 0; h >= hMin; h *= hMult, ++Nv) {
1621 ierr = VecWAXPY(uhat, h, du, u);CHKERRQ(ierr);
1622 /* F(\hat u) \approx F(u) + J(u) (uhat - u) = F(u) + h * J(u) du */
1623 ierr = SNESComputeFunction(snes, uhat, rhat);CHKERRQ(ierr);
1624 ierr = VecAXPBYPCZ(rhat, -1.0, -h, 1.0, r, df);CHKERRQ(ierr);
1625 ierr = VecNorm(rhat, NORM_2, &errors[Nv]);CHKERRQ(ierr);
1626
1627 es[Nv] = PetscLog10Real(errors[Nv]);
1628 hs[Nv] = PetscLog10Real(h);
1629 }
1630 ierr = VecDestroy(&uhat);CHKERRQ(ierr);
1631 ierr = VecDestroy(&rhat);CHKERRQ(ierr);
1632 ierr = VecDestroy(&df);CHKERRQ(ierr);
1633 ierr = VecDestroy(&r);CHKERRQ(ierr);
1634 ierr = VecDestroy(&du);CHKERRQ(ierr);
1635 for (v = 0; v < Nv; ++v) {
1636 if ((tol >= 0) && (errors[v] > tol)) break;
1637 else if (errors[v] > PETSC_SMALL) break;
1638 }
1639 if (v == Nv) isLin = PETSC_TRUE;
1640 ierr = PetscLinearRegression(Nv, hs, es, &slope, &intercept);CHKERRQ(ierr);
1641 ierr = PetscFree3(es, hs, errors);CHKERRQ(ierr);
1642 /* Slope should be about 2 */
1643 if (tol >= 0) {
1644 if (!isLin && PetscAbsReal(2 - slope) > tol) SETERRQ1(comm, PETSC_ERR_ARG_WRONG, "Taylor approximation convergence rate should be 2, not %0.2f", (double) slope);
1645 } else if (isLinear || convRate) {
1646 if (isLinear) *isLinear = isLin;
1647 if (convRate) *convRate = slope;
1648 } else {
1649 if (!isLin) {ierr = PetscPrintf(comm, "Taylor approximation converging at order %3.2f\n", (double) slope);CHKERRQ(ierr);}
1650 else {ierr = PetscPrintf(comm, "Function appears to be linear\n");CHKERRQ(ierr);}
1651 }
1652 }
1653 ierr = MatDestroy(&J);CHKERRQ(ierr);
1654 PetscFunctionReturn(0);
1655 }
1656
DMSNESCheck_Internal(SNES snes,DM dm,Vec u)1657 PetscErrorCode DMSNESCheck_Internal(SNES snes, DM dm, Vec u)
1658 {
1659 PetscErrorCode ierr;
1660
1661 PetscFunctionBegin;
1662 ierr = DMSNESCheckDiscretization(snes, dm, 0.0, u, -1.0, NULL);CHKERRQ(ierr);
1663 ierr = DMSNESCheckResidual(snes, dm, u, -1.0, NULL);CHKERRQ(ierr);
1664 ierr = DMSNESCheckJacobian(snes, dm, u, -1.0, NULL, NULL);CHKERRQ(ierr);
1665 PetscFunctionReturn(0);
1666 }
1667
1668 /*@C
1669 DMSNESCheckFromOptions - Check the residual and Jacobian functions using the exact solution by outputting some diagnostic information
1670
1671 Input Parameters:
1672 + snes - the SNES object
1673 - u - representative SNES vector
1674
1675 Note: The user must call PetscDSSetExactSolution() beforehand
1676
1677 Level: developer
1678 @*/
DMSNESCheckFromOptions(SNES snes,Vec u)1679 PetscErrorCode DMSNESCheckFromOptions(SNES snes, Vec u)
1680 {
1681 DM dm;
1682 Vec sol;
1683 PetscBool check;
1684 PetscErrorCode ierr;
1685
1686 PetscFunctionBegin;
1687 ierr = PetscOptionsHasName(((PetscObject)snes)->options,((PetscObject)snes)->prefix, "-dmsnes_check", &check);CHKERRQ(ierr);
1688 if (!check) PetscFunctionReturn(0);
1689 ierr = SNESGetDM(snes, &dm);CHKERRQ(ierr);
1690 ierr = VecDuplicate(u, &sol);CHKERRQ(ierr);
1691 ierr = SNESSetSolution(snes, sol);CHKERRQ(ierr);
1692 ierr = DMSNESCheck_Internal(snes, dm, sol);CHKERRQ(ierr);
1693 ierr = VecDestroy(&sol);CHKERRQ(ierr);
1694 PetscFunctionReturn(0);
1695 }
1696