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