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