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