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