1 #include <petsc/private/dmpleximpl.h> /*I "petscdmplex.h" I*/
2 #include <petsc/private/tsimpl.h> /*I "petscts.h" I*/
3 #include <petsc/private/snesimpl.h>
4 #include <petscds.h>
5 #include <petscfv.h>
6
DMTSConvertPlex(DM dm,DM * plex,PetscBool copy)7 static PetscErrorCode DMTSConvertPlex(DM dm, DM *plex, PetscBool copy)
8 {
9 PetscBool isPlex;
10 PetscErrorCode ierr;
11
12 PetscFunctionBegin;
13 ierr = PetscObjectTypeCompare((PetscObject) dm, DMPLEX, &isPlex);CHKERRQ(ierr);
14 if (isPlex) {
15 *plex = dm;
16 ierr = PetscObjectReference((PetscObject) dm);CHKERRQ(ierr);
17 } else {
18 ierr = PetscObjectQuery((PetscObject) dm, "dm_plex", (PetscObject *) plex);CHKERRQ(ierr);
19 if (!*plex) {
20 ierr = DMConvert(dm,DMPLEX,plex);CHKERRQ(ierr);
21 ierr = PetscObjectCompose((PetscObject) dm, "dm_plex", (PetscObject) *plex);CHKERRQ(ierr);
22 if (copy) {
23 PetscInt i;
24 PetscObject obj;
25 const char *comps[3] = {"A","dmAux","dmCh"};
26
27 ierr = DMCopyDMTS(dm, *plex);CHKERRQ(ierr);
28 ierr = DMCopyDMSNES(dm, *plex);CHKERRQ(ierr);
29 for (i = 0; i < 3; i++) {
30 ierr = PetscObjectQuery((PetscObject) dm, comps[i], &obj);CHKERRQ(ierr);
31 ierr = PetscObjectCompose((PetscObject) *plex, comps[i], obj);CHKERRQ(ierr);
32 }
33 }
34 } else {
35 ierr = PetscObjectReference((PetscObject) *plex);CHKERRQ(ierr);
36 }
37 }
38 PetscFunctionReturn(0);
39 }
40
41 /*@
42 DMPlexTSComputeRHSFunctionFVM - Form the local forcing F from the local input X using pointwise functions specified by the user
43
44 Input Parameters:
45 + dm - The mesh
46 . t - The time
47 . locX - Local solution
48 - user - The user context
49
50 Output Parameter:
51 . F - Global output vector
52
53 Level: developer
54
55 .seealso: DMPlexComputeJacobianActionFEM()
56 @*/
DMPlexTSComputeRHSFunctionFVM(DM dm,PetscReal time,Vec locX,Vec F,void * user)57 PetscErrorCode DMPlexTSComputeRHSFunctionFVM(DM dm, PetscReal time, Vec locX, Vec F, void *user)
58 {
59 Vec locF;
60 IS cellIS;
61 DM plex;
62 PetscInt depth;
63 PetscErrorCode ierr;
64
65 PetscFunctionBegin;
66 ierr = DMTSConvertPlex(dm,&plex,PETSC_TRUE);CHKERRQ(ierr);
67 ierr = DMPlexGetDepth(plex, &depth);CHKERRQ(ierr);
68 ierr = DMGetStratumIS(plex, "dim", depth, &cellIS);CHKERRQ(ierr);
69 if (!cellIS) {
70 ierr = DMGetStratumIS(plex, "depth", depth, &cellIS);CHKERRQ(ierr);
71 }
72 ierr = DMGetLocalVector(plex, &locF);CHKERRQ(ierr);
73 ierr = VecZeroEntries(locF);CHKERRQ(ierr);
74 ierr = DMPlexComputeResidual_Internal(plex, cellIS, time, locX, NULL, time, locF, user);CHKERRQ(ierr);
75 ierr = DMLocalToGlobalBegin(plex, locF, ADD_VALUES, F);CHKERRQ(ierr);
76 ierr = DMLocalToGlobalEnd(plex, locF, ADD_VALUES, F);CHKERRQ(ierr);
77 ierr = DMRestoreLocalVector(plex, &locF);CHKERRQ(ierr);
78 ierr = ISDestroy(&cellIS);CHKERRQ(ierr);
79 ierr = DMDestroy(&plex);CHKERRQ(ierr);
80 PetscFunctionReturn(0);
81 }
82
83 /*@
84 DMPlexTSComputeBoundary - Insert the essential boundary values for the local input X and/or its time derivative X_t using pointwise functions specified by the user
85
86 Input Parameters:
87 + dm - The mesh
88 . t - The time
89 . locX - Local solution
90 . locX_t - Local solution time derivative, or NULL
91 - user - The user context
92
93 Level: developer
94
95 .seealso: DMPlexComputeJacobianActionFEM()
96 @*/
DMPlexTSComputeBoundary(DM dm,PetscReal time,Vec locX,Vec locX_t,void * user)97 PetscErrorCode DMPlexTSComputeBoundary(DM dm, PetscReal time, Vec locX, Vec locX_t, void *user)
98 {
99 DM plex;
100 Vec faceGeometryFVM = NULL;
101 PetscInt Nf, f;
102 PetscErrorCode ierr;
103
104 PetscFunctionBegin;
105 ierr = DMTSConvertPlex(dm, &plex, PETSC_TRUE);CHKERRQ(ierr);
106 ierr = DMGetNumFields(plex, &Nf);CHKERRQ(ierr);
107 if (!locX_t) {
108 /* This is the RHS part */
109 for (f = 0; f < Nf; f++) {
110 PetscObject obj;
111 PetscClassId id;
112
113 ierr = DMGetField(plex, f, NULL, &obj);CHKERRQ(ierr);
114 ierr = PetscObjectGetClassId(obj, &id);CHKERRQ(ierr);
115 if (id == PETSCFV_CLASSID) {
116 ierr = DMPlexGetGeometryFVM(plex, &faceGeometryFVM, NULL, NULL);CHKERRQ(ierr);
117 break;
118 }
119 }
120 }
121 ierr = DMPlexInsertBoundaryValues(plex, PETSC_TRUE, locX, time, faceGeometryFVM, NULL, NULL);CHKERRQ(ierr);
122 ierr = DMPlexInsertTimeDerivativeBoundaryValues(plex, PETSC_TRUE, locX_t, time, faceGeometryFVM, NULL, NULL);CHKERRQ(ierr);
123 ierr = DMDestroy(&plex);CHKERRQ(ierr);
124 PetscFunctionReturn(0);
125 }
126
127 /*@
128 DMPlexTSComputeIFunctionFEM - Form the local residual F from the local input X using pointwise functions specified by the user
129
130 Input Parameters:
131 + dm - The mesh
132 . t - The time
133 . locX - Local solution
134 . locX_t - Local solution time derivative, or NULL
135 - user - The user context
136
137 Output Parameter:
138 . locF - Local output vector
139
140 Level: developer
141
142 .seealso: DMPlexComputeJacobianActionFEM()
143 @*/
DMPlexTSComputeIFunctionFEM(DM dm,PetscReal time,Vec locX,Vec locX_t,Vec locF,void * user)144 PetscErrorCode DMPlexTSComputeIFunctionFEM(DM dm, PetscReal time, Vec locX, Vec locX_t, Vec locF, void *user)
145 {
146 DM plex;
147 IS cellIS;
148 PetscInt depth;
149 PetscErrorCode ierr;
150
151 PetscFunctionBegin;
152 ierr = DMTSConvertPlex(dm,&plex,PETSC_TRUE);CHKERRQ(ierr);
153 ierr = DMPlexGetDepth(plex, &depth);CHKERRQ(ierr);
154 ierr = DMGetStratumIS(plex, "dim", depth, &cellIS);CHKERRQ(ierr);
155 if (!cellIS) {
156 ierr = DMGetStratumIS(plex, "depth", depth, &cellIS);CHKERRQ(ierr);
157 }
158 ierr = DMPlexComputeResidual_Internal(plex, cellIS, time, locX, locX_t, time, locF, user);CHKERRQ(ierr);
159 ierr = ISDestroy(&cellIS);CHKERRQ(ierr);
160 ierr = DMDestroy(&plex);CHKERRQ(ierr);
161 PetscFunctionReturn(0);
162 }
163
164 /*@
165 DMPlexTSComputeIJacobianFEM - Form the local Jacobian J from the local input X using pointwise functions specified by the user
166
167 Input Parameters:
168 + dm - The mesh
169 . t - The time
170 . locX - Local solution
171 . locX_t - Local solution time derivative, or NULL
172 . X_tshift - The multiplicative parameter for dF/du_t
173 - user - The user context
174
175 Output Parameter:
176 . locF - Local output vector
177
178 Level: developer
179
180 .seealso: DMPlexComputeJacobianActionFEM()
181 @*/
DMPlexTSComputeIJacobianFEM(DM dm,PetscReal time,Vec locX,Vec locX_t,PetscReal X_tShift,Mat Jac,Mat JacP,void * user)182 PetscErrorCode DMPlexTSComputeIJacobianFEM(DM dm, PetscReal time, Vec locX, Vec locX_t, PetscReal X_tShift, Mat Jac, Mat JacP, void *user)
183 {
184 DM plex;
185 PetscDS prob;
186 PetscBool hasJac, hasPrec;
187 IS cellIS;
188 PetscInt depth;
189 PetscErrorCode ierr;
190
191 PetscFunctionBegin;
192 ierr = DMTSConvertPlex(dm,&plex,PETSC_TRUE);CHKERRQ(ierr);
193 ierr = DMGetDS(dm, &prob);CHKERRQ(ierr);
194 ierr = PetscDSHasJacobian(prob, &hasJac);CHKERRQ(ierr);
195 ierr = PetscDSHasJacobianPreconditioner(prob, &hasPrec);CHKERRQ(ierr);
196 if (hasJac && hasPrec) {ierr = MatZeroEntries(Jac);CHKERRQ(ierr);}
197 ierr = MatZeroEntries(JacP);CHKERRQ(ierr);
198 ierr = DMPlexGetDepth(plex,&depth);CHKERRQ(ierr);
199 ierr = DMGetStratumIS(plex, "dim", depth, &cellIS);CHKERRQ(ierr);
200 if (!cellIS) {ierr = DMGetStratumIS(plex, "depth", depth, &cellIS);CHKERRQ(ierr);}
201 ierr = DMPlexComputeJacobian_Internal(plex, cellIS, time, X_tShift, locX, locX_t, Jac, JacP, user);CHKERRQ(ierr);
202 ierr = ISDestroy(&cellIS);CHKERRQ(ierr);
203 ierr = DMDestroy(&plex);CHKERRQ(ierr);
204 PetscFunctionReturn(0);
205 }
206
207 /*@C
208 DMTSCheckResidual - Check the residual of the exact solution
209
210 Input Parameters:
211 + ts - the TS object
212 . dm - the DM
213 . t - the time
214 . u - a DM vector
215 . u_t - a DM vector
216 - tol - A tolerance for the check, or -1 to print the results instead
217
218 Output Parameters:
219 . residual - The residual norm of the exact solution, or NULL
220
221 Level: developer
222
223 .seealso: DNTSCheckFromOptions(), DMTSCheckJacobian(), DNSNESCheckFromOptions(), DMSNESCheckDiscretization(), DMSNESCheckJacobian()
224 @*/
DMTSCheckResidual(TS ts,DM dm,PetscReal t,Vec u,Vec u_t,PetscReal tol,PetscReal * residual)225 PetscErrorCode DMTSCheckResidual(TS ts, DM dm, PetscReal t, Vec u, Vec u_t, PetscReal tol, PetscReal *residual)
226 {
227 MPI_Comm comm;
228 Vec r;
229 PetscReal res;
230 PetscErrorCode ierr;
231
232 PetscFunctionBegin;
233 PetscValidHeaderSpecific(ts, TS_CLASSID, 1);
234 PetscValidHeaderSpecific(dm, DM_CLASSID, 2);
235 PetscValidHeaderSpecific(u, VEC_CLASSID, 3);
236 if (residual) PetscValidRealPointer(residual, 5);
237 ierr = PetscObjectGetComm((PetscObject) ts, &comm);CHKERRQ(ierr);
238 ierr = DMComputeExactSolution(dm, t, u, u_t);CHKERRQ(ierr);
239 ierr = VecDuplicate(u, &r);CHKERRQ(ierr);
240 ierr = TSComputeIFunction(ts, t, u, u_t, r, PETSC_FALSE);CHKERRQ(ierr);
241 ierr = VecNorm(r, NORM_2, &res);CHKERRQ(ierr);
242 if (tol >= 0.0) {
243 if (res > tol) SETERRQ2(comm, PETSC_ERR_ARG_WRONG, "L_2 Residual %g exceeds tolerance %g", (double) res, (double) tol);
244 } else if (residual) {
245 *residual = res;
246 } else {
247 ierr = PetscPrintf(comm, "L_2 Residual: %g\n", (double)res);CHKERRQ(ierr);
248 ierr = VecChop(r, 1.0e-10);CHKERRQ(ierr);
249 ierr = PetscObjectCompose((PetscObject) r, "__Vec_bc_zero__", (PetscObject) dm);CHKERRQ(ierr);
250 ierr = PetscObjectSetName((PetscObject) r, "Initial Residual");CHKERRQ(ierr);
251 ierr = PetscObjectSetOptionsPrefix((PetscObject)r,"res_");CHKERRQ(ierr);
252 ierr = VecViewFromOptions(r, NULL, "-vec_view");CHKERRQ(ierr);
253 ierr = PetscObjectCompose((PetscObject) r, "__Vec_bc_zero__", NULL);CHKERRQ(ierr);
254 }
255 ierr = VecDestroy(&r);CHKERRQ(ierr);
256 PetscFunctionReturn(0);
257 }
258
259 /*@C
260 DMTSCheckJacobian - Check the Jacobian of the exact solution against the residual using the Taylor Test
261
262 Input Parameters:
263 + ts - the TS object
264 . dm - the DM
265 . t - the time
266 . u - a DM vector
267 . u_t - a DM vector
268 - tol - A tolerance for the check, or -1 to print the results instead
269
270 Output Parameters:
271 + isLinear - Flag indicaing that the function looks linear, or NULL
272 - convRate - The rate of convergence of the linear model, or NULL
273
274 Level: developer
275
276 .seealso: DNTSCheckFromOptions(), DMTSCheckResidual(), DNSNESCheckFromOptions(), DMSNESCheckDiscretization(), DMSNESCheckResidual()
277 @*/
DMTSCheckJacobian(TS ts,DM dm,PetscReal t,Vec u,Vec u_t,PetscReal tol,PetscBool * isLinear,PetscReal * convRate)278 PetscErrorCode DMTSCheckJacobian(TS ts, DM dm, PetscReal t, Vec u, Vec u_t, PetscReal tol, PetscBool *isLinear, PetscReal *convRate)
279 {
280 MPI_Comm comm;
281 PetscDS ds;
282 Mat J, M;
283 MatNullSpace nullspace;
284 PetscReal dt, shift, slope, intercept;
285 PetscBool hasJac, hasPrec, isLin = PETSC_FALSE;
286 PetscErrorCode ierr;
287
288 PetscFunctionBegin;
289 PetscValidHeaderSpecific(ts, TS_CLASSID, 1);
290 PetscValidHeaderSpecific(dm, DM_CLASSID, 2);
291 PetscValidHeaderSpecific(u, VEC_CLASSID, 3);
292 if (isLinear) PetscValidBoolPointer(isLinear, 5);
293 if (convRate) PetscValidRealPointer(convRate, 5);
294 ierr = PetscObjectGetComm((PetscObject) ts, &comm);CHKERRQ(ierr);
295 ierr = DMComputeExactSolution(dm, t, u, u_t);CHKERRQ(ierr);
296 /* Create and view matrices */
297 ierr = TSGetTimeStep(ts, &dt);CHKERRQ(ierr);
298 shift = 1.0/dt;
299 ierr = DMCreateMatrix(dm, &J);CHKERRQ(ierr);
300 ierr = DMGetDS(dm, &ds);CHKERRQ(ierr);
301 ierr = PetscDSHasJacobian(ds, &hasJac);CHKERRQ(ierr);
302 ierr = PetscDSHasJacobianPreconditioner(ds, &hasPrec);CHKERRQ(ierr);
303 if (hasJac && hasPrec) {
304 ierr = DMCreateMatrix(dm, &M);CHKERRQ(ierr);
305 ierr = TSComputeIJacobian(ts, t, u, u_t, shift, J, M, PETSC_FALSE);CHKERRQ(ierr);
306 ierr = PetscObjectSetName((PetscObject) M, "Preconditioning Matrix");CHKERRQ(ierr);
307 ierr = PetscObjectSetOptionsPrefix((PetscObject) M, "jacpre_");CHKERRQ(ierr);
308 ierr = MatViewFromOptions(M, NULL, "-mat_view");CHKERRQ(ierr);
309 ierr = MatDestroy(&M);CHKERRQ(ierr);
310 } else {
311 ierr = TSComputeIJacobian(ts, t, u, u_t, shift, J, J, PETSC_FALSE);CHKERRQ(ierr);
312 }
313 ierr = PetscObjectSetName((PetscObject) J, "Jacobian");CHKERRQ(ierr);
314 ierr = PetscObjectSetOptionsPrefix((PetscObject) J, "jac_");CHKERRQ(ierr);
315 ierr = MatViewFromOptions(J, NULL, "-mat_view");CHKERRQ(ierr);
316 /* Check nullspace */
317 ierr = MatGetNullSpace(J, &nullspace);CHKERRQ(ierr);
318 if (nullspace) {
319 PetscBool isNull;
320 ierr = MatNullSpaceTest(nullspace, J, &isNull);CHKERRQ(ierr);
321 if (!isNull) SETERRQ(comm, PETSC_ERR_PLIB, "The null space calculated for the system operator is invalid.");
322 }
323 ierr = MatNullSpaceDestroy(&nullspace);CHKERRQ(ierr);
324 /* Taylor test */
325 {
326 PetscRandom rand;
327 Vec du, uhat, uhat_t, r, rhat, df;
328 PetscReal h;
329 PetscReal *es, *hs, *errors;
330 PetscReal hMax = 1.0, hMin = 1e-6, hMult = 0.1;
331 PetscInt Nv, v;
332
333 /* Choose a perturbation direction */
334 ierr = PetscRandomCreate(comm, &rand);CHKERRQ(ierr);
335 ierr = VecDuplicate(u, &du);CHKERRQ(ierr);
336 ierr = VecSetRandom(du, rand);CHKERRQ(ierr);
337 ierr = PetscRandomDestroy(&rand);CHKERRQ(ierr);
338 ierr = VecDuplicate(u, &df);CHKERRQ(ierr);
339 ierr = MatMult(J, du, df);CHKERRQ(ierr);
340 /* Evaluate residual at u, F(u), save in vector r */
341 ierr = VecDuplicate(u, &r);CHKERRQ(ierr);
342 ierr = TSComputeIFunction(ts, t, u, u_t, r, PETSC_FALSE);CHKERRQ(ierr);
343 /* Look at the convergence of our Taylor approximation as we approach u */
344 for (h = hMax, Nv = 0; h >= hMin; h *= hMult, ++Nv);
345 ierr = PetscCalloc3(Nv, &es, Nv, &hs, Nv, &errors);CHKERRQ(ierr);
346 ierr = VecDuplicate(u, &uhat);CHKERRQ(ierr);
347 ierr = VecDuplicate(u, &uhat_t);CHKERRQ(ierr);
348 ierr = VecDuplicate(u, &rhat);CHKERRQ(ierr);
349 for (h = hMax, Nv = 0; h >= hMin; h *= hMult, ++Nv) {
350 ierr = VecWAXPY(uhat, h, du, u);CHKERRQ(ierr);
351 ierr = VecWAXPY(uhat_t, h*shift, du, u_t);CHKERRQ(ierr);
352 /* F(\hat u, \hat u_t) \approx F(u, u_t) + J(u, u_t) (uhat - u) + J_t(u, u_t) (uhat_t - u_t) = F(u) + h * J(u) du + h * shift * J_t(u) du = F(u) + h F' du */
353 ierr = TSComputeIFunction(ts, t, uhat, uhat_t, rhat, PETSC_FALSE);CHKERRQ(ierr);
354 ierr = VecAXPBYPCZ(rhat, -1.0, -h, 1.0, r, df);CHKERRQ(ierr);
355 ierr = VecNorm(rhat, NORM_2, &errors[Nv]);CHKERRQ(ierr);
356
357 es[Nv] = PetscLog10Real(errors[Nv]);
358 hs[Nv] = PetscLog10Real(h);
359 }
360 ierr = VecDestroy(&uhat);CHKERRQ(ierr);
361 ierr = VecDestroy(&uhat_t);CHKERRQ(ierr);
362 ierr = VecDestroy(&rhat);CHKERRQ(ierr);
363 ierr = VecDestroy(&df);CHKERRQ(ierr);
364 ierr = VecDestroy(&r);CHKERRQ(ierr);
365 ierr = VecDestroy(&du);CHKERRQ(ierr);
366 for (v = 0; v < Nv; ++v) {
367 if ((tol >= 0) && (errors[v] > tol)) break;
368 else if (errors[v] > PETSC_SMALL) break;
369 }
370 if (v == Nv) isLin = PETSC_TRUE;
371 ierr = PetscLinearRegression(Nv, hs, es, &slope, &intercept);CHKERRQ(ierr);
372 ierr = PetscFree3(es, hs, errors);CHKERRQ(ierr);
373 /* Slope should be about 2 */
374 if (tol >= 0) {
375 if (!isLin && PetscAbsReal(2 - slope) > tol) SETERRQ1(comm, PETSC_ERR_ARG_WRONG, "Taylor approximation convergence rate should be 2, not %0.2f", (double) slope);
376 } else if (isLinear || convRate) {
377 if (isLinear) *isLinear = isLin;
378 if (convRate) *convRate = slope;
379 } else {
380 if (!isLin) {ierr = PetscPrintf(comm, "Taylor approximation converging at order %3.2f\n", (double) slope);CHKERRQ(ierr);}
381 else {ierr = PetscPrintf(comm, "Function appears to be linear\n");CHKERRQ(ierr);}
382 }
383 }
384 ierr = MatDestroy(&J);CHKERRQ(ierr);
385 PetscFunctionReturn(0);
386 }
387
388 /*@C
389 DMTSCheckFromOptions - Check the residual and Jacobian functions using the exact solution by outputting some diagnostic information
390
391 Input Parameters:
392 + ts - the TS object
393 - u - representative TS vector
394
395 Note: The user must call PetscDSSetExactSolution() beforehand
396
397 Level: developer
398 @*/
DMTSCheckFromOptions(TS ts,Vec u)399 PetscErrorCode DMTSCheckFromOptions(TS ts, Vec u)
400 {
401 DM dm;
402 SNES snes;
403 Vec sol, u_t;
404 PetscReal t;
405 PetscBool check;
406 PetscErrorCode ierr;
407
408 PetscFunctionBegin;
409 ierr = PetscOptionsHasName(((PetscObject)ts)->options,((PetscObject)ts)->prefix, "-dmts_check", &check);CHKERRQ(ierr);
410 if (!check) PetscFunctionReturn(0);
411 ierr = VecDuplicate(u, &sol);CHKERRQ(ierr);
412 ierr = VecCopy(u, sol);CHKERRQ(ierr);
413 ierr = TSSetSolution(ts, u);CHKERRQ(ierr);
414 ierr = TSGetDM(ts, &dm);CHKERRQ(ierr);
415 ierr = TSSetUp(ts);CHKERRQ(ierr);
416 ierr = TSGetSNES(ts, &snes);CHKERRQ(ierr);
417 ierr = SNESSetSolution(snes, u);CHKERRQ(ierr);
418
419 ierr = TSGetTime(ts, &t);CHKERRQ(ierr);
420 ierr = DMSNESCheckDiscretization(snes, dm, t, sol, -1.0, NULL);CHKERRQ(ierr);
421 ierr = DMGetGlobalVector(dm, &u_t);CHKERRQ(ierr);
422 ierr = DMTSCheckResidual(ts, dm, t, sol, u_t, -1.0, NULL);CHKERRQ(ierr);
423 ierr = DMTSCheckJacobian(ts, dm, t, sol, u_t, -1.0, NULL, NULL);CHKERRQ(ierr);
424 ierr = DMRestoreGlobalVector(dm, &u_t);CHKERRQ(ierr);
425
426 ierr = VecDestroy(&sol);CHKERRQ(ierr);
427 PetscFunctionReturn(0);
428 }
429