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