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