1 /* 2 Defines a SNES that can consist of a collection of SNESes on patches of the domain 3 */ 4 #include <petsc/private/snesimpl.h> /*I "petscsnes.h" I*/ 5 #include <petsc/private/pcpatchimpl.h> /* We need internal access to PCPatch right now, until that part is moved to Plex */ 6 #include <petscsf.h> 7 #include <petscsection.h> 8 9 typedef struct { 10 PC pc; /* The linear patch preconditioner */ 11 } SNES_Patch; 12 13 static PetscErrorCode SNESPatchComputeResidual_Private(SNES snes, Vec x, Vec F, void *ctx) 14 { 15 PC pc = (PC) ctx; 16 PC_PATCH *pcpatch = (PC_PATCH *) pc->data; 17 PetscInt pt, size, i; 18 const PetscInt *indices; 19 const PetscScalar *X; 20 PetscScalar *XWithAll; 21 PetscErrorCode ierr; 22 23 PetscFunctionBegin; 24 25 /* scatter from x to patch->patchStateWithAll[pt] */ 26 pt = pcpatch->currentPatch; 27 ierr = ISGetSize(pcpatch->dofMappingWithoutToWithAll[pt], &size);CHKERRQ(ierr); 28 29 ierr = ISGetIndices(pcpatch->dofMappingWithoutToWithAll[pt], &indices);CHKERRQ(ierr); 30 ierr = VecGetArrayRead(x, &X);CHKERRQ(ierr); 31 ierr = VecGetArray(pcpatch->patchStateWithAll, &XWithAll);CHKERRQ(ierr); 32 33 for (i = 0; i < size; ++i) { 34 XWithAll[indices[i]] = X[i]; 35 } 36 37 ierr = VecRestoreArray(pcpatch->patchStateWithAll, &XWithAll);CHKERRQ(ierr); 38 ierr = VecRestoreArrayRead(x, &X);CHKERRQ(ierr); 39 ierr = ISRestoreIndices(pcpatch->dofMappingWithoutToWithAll[pt], &indices);CHKERRQ(ierr); 40 41 ierr = PCPatchComputeFunction_Internal(pc, pcpatch->patchStateWithAll, F, pt);CHKERRQ(ierr); 42 PetscFunctionReturn(0); 43 } 44 45 static PetscErrorCode SNESPatchComputeJacobian_Private(SNES snes, Vec x, Mat J, Mat M, void *ctx) 46 { 47 PC pc = (PC) ctx; 48 PC_PATCH *pcpatch = (PC_PATCH *) pc->data; 49 PetscInt pt, size, i; 50 const PetscInt *indices; 51 const PetscScalar *X; 52 PetscScalar *XWithAll; 53 PetscErrorCode ierr; 54 55 PetscFunctionBegin; 56 /* scatter from x to patch->patchStateWithAll[pt] */ 57 pt = pcpatch->currentPatch; 58 ierr = ISGetSize(pcpatch->dofMappingWithoutToWithAll[pt], &size);CHKERRQ(ierr); 59 60 ierr = ISGetIndices(pcpatch->dofMappingWithoutToWithAll[pt], &indices);CHKERRQ(ierr); 61 ierr = VecGetArrayRead(x, &X);CHKERRQ(ierr); 62 ierr = VecGetArray(pcpatch->patchStateWithAll, &XWithAll);CHKERRQ(ierr); 63 64 for (i = 0; i < size; ++i) { 65 XWithAll[indices[i]] = X[i]; 66 } 67 68 ierr = VecRestoreArray(pcpatch->patchStateWithAll, &XWithAll);CHKERRQ(ierr); 69 ierr = VecRestoreArrayRead(x, &X);CHKERRQ(ierr); 70 ierr = ISRestoreIndices(pcpatch->dofMappingWithoutToWithAll[pt], &indices);CHKERRQ(ierr); 71 72 ierr = PCPatchComputeOperator_Internal(pc, pcpatch->patchStateWithAll, M, pcpatch->currentPatch, PETSC_FALSE);CHKERRQ(ierr); 73 PetscFunctionReturn(0); 74 } 75 76 static PetscErrorCode PCSetUp_PATCH_Nonlinear(PC pc) 77 { 78 PC_PATCH *patch = (PC_PATCH *) pc->data; 79 const char *prefix; 80 PetscInt i, pStart, dof, maxDof = -1; 81 PetscErrorCode ierr; 82 83 PetscFunctionBegin; 84 if (!pc->setupcalled) { 85 ierr = PetscMalloc1(patch->npatch, &patch->solver);CHKERRQ(ierr); 86 ierr = PCGetOptionsPrefix(pc, &prefix);CHKERRQ(ierr); 87 ierr = PetscSectionGetChart(patch->gtolCounts, &pStart, NULL);CHKERRQ(ierr); 88 for (i = 0; i < patch->npatch; ++i) { 89 SNES snes; 90 91 ierr = SNESCreate(PETSC_COMM_SELF, &snes);CHKERRQ(ierr); 92 ierr = SNESSetOptionsPrefix(snes, prefix);CHKERRQ(ierr); 93 ierr = SNESAppendOptionsPrefix(snes, "sub_");CHKERRQ(ierr); 94 ierr = PetscObjectIncrementTabLevel((PetscObject) snes, (PetscObject) pc, 2);CHKERRQ(ierr); 95 ierr = PetscLogObjectParent((PetscObject) pc, (PetscObject) snes);CHKERRQ(ierr); 96 patch->solver[i] = (PetscObject) snes; 97 98 ierr = PetscSectionGetDof(patch->gtolCountsWithAll, i+pStart, &dof);CHKERRQ(ierr); 99 maxDof = PetscMax(maxDof, dof); 100 } 101 ierr = VecDuplicate(patch->localUpdate, &patch->localState);CHKERRQ(ierr); 102 ierr = VecDuplicate(patch->patchRHS, &patch->patchResidual);CHKERRQ(ierr); 103 ierr = VecDuplicate(patch->patchUpdate, &patch->patchState);CHKERRQ(ierr); 104 105 ierr = VecCreateSeq(PETSC_COMM_SELF, maxDof, &patch->patchStateWithAll);CHKERRQ(ierr); 106 ierr = VecSetUp(patch->patchStateWithAll);CHKERRQ(ierr); 107 } 108 for (i = 0; i < patch->npatch; ++i) { 109 SNES snes = (SNES) patch->solver[i]; 110 111 ierr = SNESSetFunction(snes, patch->patchResidual, SNESPatchComputeResidual_Private, pc);CHKERRQ(ierr); 112 ierr = SNESSetJacobian(snes, patch->mat[i], patch->mat[i], SNESPatchComputeJacobian_Private, pc);CHKERRQ(ierr); 113 } 114 if (!pc->setupcalled && patch->optionsSet) for (i = 0; i < patch->npatch; ++i) {ierr = SNESSetFromOptions((SNES) patch->solver[i]);CHKERRQ(ierr);} 115 PetscFunctionReturn(0); 116 } 117 118 static PetscErrorCode PCApply_PATCH_Nonlinear(PC pc, PetscInt i, Vec patchRHS, Vec patchUpdate) 119 { 120 PC_PATCH *patch = (PC_PATCH *) pc->data; 121 PetscInt pStart; 122 PetscErrorCode ierr; 123 124 PetscFunctionBegin; 125 patch->currentPatch = i; 126 ierr = PetscLogEventBegin(PC_Patch_Solve, pc, 0, 0, 0);CHKERRQ(ierr); 127 128 /* Scatter the overlapped global state to our patch state vector */ 129 ierr = PetscSectionGetChart(patch->gtolCounts, &pStart, NULL);CHKERRQ(ierr); 130 ierr = PCPatch_ScatterLocal_Private(pc, i+pStart, patch->localState, patch->patchState, INSERT_VALUES, SCATTER_FORWARD, SCATTER_INTERIOR);CHKERRQ(ierr); 131 ierr = PCPatch_ScatterLocal_Private(pc, i+pStart, patch->localState, patch->patchStateWithAll, INSERT_VALUES, SCATTER_FORWARD, SCATTER_WITHALL);CHKERRQ(ierr); 132 133 /* Set initial guess to be current state*/ 134 ierr = VecCopy(patch->patchState, patchUpdate);CHKERRQ(ierr); 135 /* Solve for new state */ 136 ierr = SNESSolve((SNES) patch->solver[i], patchRHS, patchUpdate);CHKERRQ(ierr); 137 /* To compute update, subtract off previous state */ 138 ierr = VecAXPY(patchUpdate, -1.0, patch->patchState);CHKERRQ(ierr); 139 140 ierr = PetscLogEventEnd(PC_Patch_Solve, pc, 0, 0, 0);CHKERRQ(ierr); 141 PetscFunctionReturn(0); 142 } 143 144 static PetscErrorCode PCReset_PATCH_Nonlinear(PC pc) 145 { 146 PC_PATCH *patch = (PC_PATCH *) pc->data; 147 PetscInt i; 148 PetscErrorCode ierr; 149 150 PetscFunctionBegin; 151 if (patch->solver) { 152 for (i = 0; i < patch->npatch; ++i) {ierr = SNESReset((SNES) patch->solver[i]);CHKERRQ(ierr);} 153 } 154 155 ierr = VecDestroy(&patch->patchResidual);CHKERRQ(ierr); 156 ierr = VecDestroy(&patch->patchState);CHKERRQ(ierr); 157 ierr = VecDestroy(&patch->patchStateWithAll);CHKERRQ(ierr); 158 159 ierr = VecDestroy(&patch->localState);CHKERRQ(ierr); 160 PetscFunctionReturn(0); 161 } 162 163 static PetscErrorCode PCDestroy_PATCH_Nonlinear(PC pc) 164 { 165 PC_PATCH *patch = (PC_PATCH *) pc->data; 166 PetscInt i; 167 PetscErrorCode ierr; 168 169 PetscFunctionBegin; 170 if (patch->solver) { 171 for (i = 0; i < patch->npatch; ++i) {ierr = SNESDestroy((SNES *) &patch->solver[i]);CHKERRQ(ierr);} 172 ierr = PetscFree(patch->solver);CHKERRQ(ierr); 173 } 174 PetscFunctionReturn(0); 175 } 176 177 static PetscErrorCode PCUpdateMultiplicative_PATCH_Nonlinear(PC pc, PetscInt i, PetscInt pStart) 178 { 179 PC_PATCH *patch = (PC_PATCH *) pc->data; 180 PetscErrorCode ierr; 181 182 PetscFunctionBegin; 183 ierr = PCPatch_ScatterLocal_Private(pc, i + pStart, patch->patchUpdate, patch->localState, ADD_VALUES, SCATTER_REVERSE, SCATTER_INTERIOR);CHKERRQ(ierr); 184 PetscFunctionReturn(0); 185 } 186 187 static PetscErrorCode SNESSetUp_Patch(SNES snes) 188 { 189 SNES_Patch *patch = (SNES_Patch *) snes->data; 190 DM dm; 191 Mat dummy; 192 Vec F; 193 PetscInt n, N; 194 PetscErrorCode ierr; 195 196 PetscFunctionBegin; 197 ierr = SNESGetDM(snes, &dm);CHKERRQ(ierr); 198 ierr = PCSetDM(patch->pc, dm);CHKERRQ(ierr); 199 ierr = SNESGetFunction(snes, &F, NULL, NULL);CHKERRQ(ierr); 200 ierr = VecGetLocalSize(F, &n);CHKERRQ(ierr); 201 ierr = VecGetSize(F, &N);CHKERRQ(ierr); 202 ierr = MatCreateShell(PetscObjectComm((PetscObject) snes), n, n, N, N, (void *) snes, &dummy);CHKERRQ(ierr); 203 ierr = PCSetOperators(patch->pc, dummy, dummy);CHKERRQ(ierr); 204 ierr = MatDestroy(&dummy);CHKERRQ(ierr); 205 ierr = PCSetUp(patch->pc);CHKERRQ(ierr); 206 /* allocate workspace */ 207 PetscFunctionReturn(0); 208 } 209 210 static PetscErrorCode SNESReset_Patch(SNES snes) 211 { 212 SNES_Patch *patch = (SNES_Patch *) snes->data; 213 PetscErrorCode ierr; 214 215 PetscFunctionBegin; 216 ierr = PCReset(patch->pc);CHKERRQ(ierr); 217 PetscFunctionReturn(0); 218 } 219 220 static PetscErrorCode SNESDestroy_Patch(SNES snes) 221 { 222 SNES_Patch *patch = (SNES_Patch *) snes->data; 223 PetscErrorCode ierr; 224 225 PetscFunctionBegin; 226 ierr = SNESReset_Patch(snes);CHKERRQ(ierr); 227 ierr = PCDestroy(&patch->pc);CHKERRQ(ierr); 228 ierr = PetscFree(snes->data);CHKERRQ(ierr); 229 PetscFunctionReturn(0); 230 } 231 232 static PetscErrorCode SNESSetFromOptions_Patch(PetscOptionItems *PetscOptionsObject, SNES snes) 233 { 234 SNES_Patch *patch = (SNES_Patch *) snes->data; 235 const char *prefix; 236 PetscErrorCode ierr; 237 238 PetscFunctionBegin; 239 ierr = PetscObjectGetOptionsPrefix((PetscObject)snes, &prefix);CHKERRQ(ierr); 240 ierr = PetscObjectSetOptionsPrefix((PetscObject)patch->pc, prefix);CHKERRQ(ierr); 241 ierr = PCSetFromOptions(patch->pc);CHKERRQ(ierr); 242 PetscFunctionReturn(0); 243 } 244 245 static PetscErrorCode SNESView_Patch(SNES snes,PetscViewer viewer) 246 { 247 SNES_Patch *patch = (SNES_Patch *) snes->data; 248 PetscBool iascii; 249 PetscErrorCode ierr; 250 251 PetscFunctionBegin; 252 ierr = PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &iascii);CHKERRQ(ierr); 253 if (iascii) { 254 ierr = PetscViewerASCIIPrintf(viewer,"SNESPATCH\n");CHKERRQ(ierr); 255 } 256 ierr = PetscViewerASCIIPushTab(viewer);CHKERRQ(ierr); 257 ierr = PCView(patch->pc, viewer);CHKERRQ(ierr); 258 ierr = PetscViewerASCIIPopTab(viewer);CHKERRQ(ierr); 259 PetscFunctionReturn(0); 260 } 261 262 static PetscErrorCode SNESSolve_Patch(SNES snes) 263 { 264 SNES_Patch *patch = (SNES_Patch *) snes->data; 265 PC_PATCH *pcpatch = (PC_PATCH *) patch->pc->data; 266 SNESLineSearch ls; 267 Vec rhs, update, state, residual; 268 const PetscScalar *globalState = NULL; 269 PetscScalar *localState = NULL; 270 PetscInt its = 0; 271 PetscReal xnorm = 0.0, ynorm = 0.0, fnorm = 0.0; 272 PetscErrorCode ierr; 273 274 PetscFunctionBegin; 275 ierr = SNESGetSolution(snes, &state);CHKERRQ(ierr); 276 ierr = SNESGetSolutionUpdate(snes, &update);CHKERRQ(ierr); 277 ierr = SNESGetRhs(snes, &rhs);CHKERRQ(ierr); 278 279 ierr = SNESGetFunction(snes, &residual, NULL, NULL);CHKERRQ(ierr); 280 ierr = SNESGetLineSearch(snes, &ls);CHKERRQ(ierr); 281 282 ierr = SNESSetConvergedReason(snes, SNES_CONVERGED_ITERATING);CHKERRQ(ierr); 283 ierr = VecSet(update, 0.0);CHKERRQ(ierr); 284 ierr = SNESComputeFunction(snes, state, residual);CHKERRQ(ierr); 285 286 ierr = VecNorm(state, NORM_2, &xnorm);CHKERRQ(ierr); 287 ierr = VecNorm(residual, NORM_2, &fnorm);CHKERRQ(ierr); 288 snes->ttol = fnorm*snes->rtol; 289 290 if (snes->ops->converged) { 291 ierr = (*snes->ops->converged)(snes,its,xnorm,ynorm,fnorm,&snes->reason,snes->cnvP);CHKERRQ(ierr); 292 } else { 293 ierr = SNESConvergedSkip(snes,its,xnorm,ynorm,fnorm,&snes->reason,0);CHKERRQ(ierr); 294 } 295 ierr = SNESLogConvergenceHistory(snes, fnorm, 0);CHKERRQ(ierr); /* should we count lits from the patches? */ 296 ierr = SNESMonitor(snes, its, fnorm);CHKERRQ(ierr); 297 298 /* The main solver loop */ 299 for (its = 0; its < snes->max_its; its++) { 300 301 ierr = SNESSetIterationNumber(snes, its);CHKERRQ(ierr); 302 303 /* Scatter state vector to overlapped vector on all patches. 304 The vector pcpatch->localState is scattered to each patch 305 in PCApply_PATCH_Nonlinear. */ 306 ierr = VecGetArrayRead(state, &globalState);CHKERRQ(ierr); 307 ierr = VecGetArray(pcpatch->localState, &localState);CHKERRQ(ierr); 308 ierr = PetscSFBcastBegin(pcpatch->sectionSF, MPIU_SCALAR, globalState, localState);CHKERRQ(ierr); 309 ierr = PetscSFBcastEnd(pcpatch->sectionSF, MPIU_SCALAR, globalState, localState);CHKERRQ(ierr); 310 ierr = VecRestoreArray(pcpatch->localState, &localState);CHKERRQ(ierr); 311 ierr = VecRestoreArrayRead(state, &globalState);CHKERRQ(ierr); 312 313 /* The looping over patches happens here */ 314 ierr = PCApply(patch->pc, rhs, update);CHKERRQ(ierr); 315 316 /* Apply a line search. This will often be basic with 317 damping = 1/(max number of patches a dof can be in), 318 but not always */ 319 ierr = VecScale(update, -1.0);CHKERRQ(ierr); 320 ierr = SNESLineSearchApply(ls, state, residual, &fnorm, update);CHKERRQ(ierr); 321 322 ierr = VecNorm(state, NORM_2, &xnorm);CHKERRQ(ierr); 323 ierr = VecNorm(update, NORM_2, &ynorm);CHKERRQ(ierr); 324 325 if (snes->ops->converged) { 326 ierr = (*snes->ops->converged)(snes,its,xnorm,ynorm,fnorm,&snes->reason,snes->cnvP);CHKERRQ(ierr); 327 } else { 328 ierr = SNESConvergedSkip(snes,its,xnorm,ynorm,fnorm,&snes->reason,0);CHKERRQ(ierr); 329 } 330 ierr = SNESLogConvergenceHistory(snes, fnorm, 0);CHKERRQ(ierr); /* FIXME: should we count lits? */ 331 ierr = SNESMonitor(snes, its, fnorm);CHKERRQ(ierr); 332 } 333 334 if (its == snes->max_its) { ierr = SNESSetConvergedReason(snes, SNES_DIVERGED_MAX_IT);CHKERRQ(ierr); } 335 PetscFunctionReturn(0); 336 } 337 338 /*MC 339 SNESPATCH - Solve a nonlinear problem by composing together many nonlinear solvers on patches 340 341 Level: intermediate 342 343 .seealso: SNESCreate(), SNESSetType(), SNESType (for list of available types), SNES, 344 PCPATCH 345 346 References: 347 . 1. - Peter R. Brune, Matthew G. Knepley, Barry F. Smith, and Xuemin Tu, "Composing Scalable Nonlinear Algebraic Solvers", SIAM Review, 57(4), 2015 348 349 M*/ 350 PETSC_EXTERN PetscErrorCode SNESCreate_Patch(SNES snes) 351 { 352 PetscErrorCode ierr; 353 SNES_Patch *patch; 354 PC_PATCH *patchpc; 355 SNESLineSearch linesearch; 356 357 PetscFunctionBegin; 358 ierr = PetscNewLog(snes, &patch);CHKERRQ(ierr); 359 360 snes->ops->solve = SNESSolve_Patch; 361 snes->ops->setup = SNESSetUp_Patch; 362 snes->ops->reset = SNESReset_Patch; 363 snes->ops->destroy = SNESDestroy_Patch; 364 snes->ops->setfromoptions = SNESSetFromOptions_Patch; 365 snes->ops->view = SNESView_Patch; 366 367 ierr = SNESGetLineSearch(snes,&linesearch);CHKERRQ(ierr); 368 ierr = SNESLineSearchSetType(linesearch,SNESLINESEARCHBASIC);CHKERRQ(ierr); 369 snes->usesksp = PETSC_FALSE; 370 371 snes->alwayscomputesfinalresidual = PETSC_FALSE; 372 373 snes->data = (void *) patch; 374 ierr = PCCreate(PetscObjectComm((PetscObject) snes), &patch->pc);CHKERRQ(ierr); 375 ierr = PCSetType(patch->pc, PCPATCH);CHKERRQ(ierr); 376 377 patchpc = (PC_PATCH*) patch->pc->data; 378 patchpc->classname = "snes"; 379 patchpc->isNonlinear = PETSC_TRUE; 380 381 patchpc->setupsolver = PCSetUp_PATCH_Nonlinear; 382 patchpc->applysolver = PCApply_PATCH_Nonlinear; 383 patchpc->resetsolver = PCReset_PATCH_Nonlinear; 384 patchpc->destroysolver = PCDestroy_PATCH_Nonlinear; 385 patchpc->updatemultiplicative = PCUpdateMultiplicative_PATCH_Nonlinear; 386 387 PetscFunctionReturn(0); 388 } 389 390 PetscErrorCode SNESPatchSetDiscretisationInfo(SNES snes, PetscInt nsubspaces, DM *dms, PetscInt *bs, PetscInt *nodesPerCell, const PetscInt **cellNodeMap, 391 const PetscInt *subspaceOffsets, PetscInt numGhostBcs, const PetscInt *ghostBcNodes, PetscInt numGlobalBcs, const PetscInt *globalBcNodes) 392 { 393 SNES_Patch *patch = (SNES_Patch *) snes->data; 394 PetscErrorCode ierr; 395 DM dm; 396 397 PetscFunctionBegin; 398 ierr = SNESGetDM(snes, &dm);CHKERRQ(ierr); 399 if (!dm) SETERRQ(PetscObjectComm((PetscObject)snes), PETSC_ERR_ARG_WRONGSTATE, "DM not yet set on patch SNES\n"); 400 ierr = PCSetDM(patch->pc, dm);CHKERRQ(ierr); 401 ierr = PCPatchSetDiscretisationInfo(patch->pc, nsubspaces, dms, bs, nodesPerCell, cellNodeMap, subspaceOffsets, numGhostBcs, ghostBcNodes, numGlobalBcs, globalBcNodes);CHKERRQ(ierr); 402 PetscFunctionReturn(0); 403 } 404 405 PetscErrorCode SNESPatchSetComputeOperator(SNES snes, PetscErrorCode (*func)(PC, PetscInt, Vec, Mat, IS, PetscInt, const PetscInt *, const PetscInt *, void *), void *ctx) 406 { 407 SNES_Patch *patch = (SNES_Patch *) snes->data; 408 PetscErrorCode ierr; 409 410 PetscFunctionBegin; 411 ierr = PCPatchSetComputeOperator(patch->pc, func, ctx);CHKERRQ(ierr); 412 PetscFunctionReturn(0); 413 } 414 415 PetscErrorCode SNESPatchSetComputeFunction(SNES snes, PetscErrorCode (*func)(PC, PetscInt, Vec, Vec, IS, PetscInt, const PetscInt *, const PetscInt *, void *), void *ctx) 416 { 417 SNES_Patch *patch = (SNES_Patch *) snes->data; 418 PetscErrorCode ierr; 419 420 PetscFunctionBegin; 421 ierr = PCPatchSetComputeFunction(patch->pc, func, ctx);CHKERRQ(ierr); 422 PetscFunctionReturn(0); 423 } 424 425 PetscErrorCode SNESPatchSetConstructType(SNES snes, PCPatchConstructType ctype, PetscErrorCode (*func)(PC, PetscInt *, IS **, IS *, void *), void *ctx) 426 { 427 SNES_Patch *patch = (SNES_Patch *) snes->data; 428 PetscErrorCode ierr; 429 430 PetscFunctionBegin; 431 ierr = PCPatchSetConstructType(patch->pc, ctype, func, ctx);CHKERRQ(ierr); 432 PetscFunctionReturn(0); 433 } 434 435 PetscErrorCode SNESPatchSetCellNumbering(SNES snes, PetscSection cellNumbering) 436 { 437 SNES_Patch *patch = (SNES_Patch *) snes->data; 438 PetscErrorCode ierr; 439 440 PetscFunctionBegin; 441 ierr = PCPatchSetCellNumbering(patch->pc, cellNumbering);CHKERRQ(ierr); 442 PetscFunctionReturn(0); 443 } 444