1 2 #include <../src/snes/impls/vi/rs/virsimpl.h> /*I "petscsnes.h" I*/ 3 #include <petsc/private/dmimpl.h> 4 #include <petsc/private/vecimpl.h> 5 6 /* 7 SNESVIGetInactiveSet - Gets the global indices for the inactive set variables (these correspond to the degrees of freedom the linear 8 system is solved on) 9 10 Input parameter: 11 . snes - the SNES context 12 13 Output parameter: 14 . inact - inactive set index set 15 16 */ 17 PetscErrorCode SNESVIGetInactiveSet(SNES snes,IS *inact) 18 { 19 SNES_VINEWTONRSLS *vi = (SNES_VINEWTONRSLS*)snes->data; 20 21 PetscFunctionBegin; 22 *inact = vi->IS_inact; 23 PetscFunctionReturn(0); 24 } 25 26 /* 27 Provides a wrapper to a DM to allow it to be used to generated the interpolation/restriction from the DM for the smaller matrices and vectors 28 defined by the reduced space method. 29 30 Simple calls the regular DM interpolation and restricts it to operation on the variables not associated with active constraints. 31 32 */ 33 typedef struct { 34 PetscInt n; /* size of vectors in the reduced DM space */ 35 IS inactive; 36 37 PetscErrorCode (*createinterpolation)(DM,DM,Mat*,Vec*); /* DM's original routines */ 38 PetscErrorCode (*coarsen)(DM, MPI_Comm, DM*); 39 PetscErrorCode (*createglobalvector)(DM,Vec*); 40 PetscErrorCode (*createinjection)(DM,DM,Mat*); 41 PetscErrorCode (*hascreateinjection)(DM,PetscBool*); 42 43 DM dm; /* when destroying this object we need to reset the above function into the base DM */ 44 } DM_SNESVI; 45 46 /* 47 DMCreateGlobalVector_SNESVI - Creates global vector of the size of the reduced space 48 49 */ 50 PetscErrorCode DMCreateGlobalVector_SNESVI(DM dm,Vec *vec) 51 { 52 PetscContainer isnes; 53 DM_SNESVI *dmsnesvi; 54 55 PetscFunctionBegin; 56 CHKERRQ(PetscObjectQuery((PetscObject)dm,"VI",(PetscObject*)&isnes)); 57 PetscCheck(isnes,PetscObjectComm((PetscObject)dm),PETSC_ERR_PLIB,"Composed SNES is missing"); 58 CHKERRQ(PetscContainerGetPointer(isnes,(void**)&dmsnesvi)); 59 CHKERRQ(VecCreateMPI(PetscObjectComm((PetscObject)dm),dmsnesvi->n,PETSC_DETERMINE,vec)); 60 CHKERRQ(VecSetDM(*vec, dm)); 61 PetscFunctionReturn(0); 62 } 63 64 /* 65 DMCreateInterpolation_SNESVI - Modifieds the interpolation obtained from the DM by removing all rows and columns associated with active constraints. 66 67 */ 68 PetscErrorCode DMCreateInterpolation_SNESVI(DM dm1,DM dm2,Mat *mat,Vec *vec) 69 { 70 PetscContainer isnes; 71 DM_SNESVI *dmsnesvi1,*dmsnesvi2; 72 Mat interp; 73 74 PetscFunctionBegin; 75 CHKERRQ(PetscObjectQuery((PetscObject)dm1,"VI",(PetscObject*)&isnes)); 76 PetscCheck(isnes,PetscObjectComm((PetscObject)dm1),PETSC_ERR_PLIB,"Composed VI data structure is missing"); 77 CHKERRQ(PetscContainerGetPointer(isnes,(void**)&dmsnesvi1)); 78 CHKERRQ(PetscObjectQuery((PetscObject)dm2,"VI",(PetscObject*)&isnes)); 79 PetscCheck(isnes,PetscObjectComm((PetscObject)dm2),PETSC_ERR_PLIB,"Composed VI data structure is missing"); 80 CHKERRQ(PetscContainerGetPointer(isnes,(void**)&dmsnesvi2)); 81 82 CHKERRQ((*dmsnesvi1->createinterpolation)(dm1,dm2,&interp,NULL)); 83 CHKERRQ(MatCreateSubMatrix(interp,dmsnesvi2->inactive,dmsnesvi1->inactive,MAT_INITIAL_MATRIX,mat)); 84 CHKERRQ(MatDestroy(&interp)); 85 *vec = NULL; 86 PetscFunctionReturn(0); 87 } 88 89 static PetscErrorCode DMSetVI(DM,IS); 90 static PetscErrorCode DMDestroyVI(DM); 91 92 /* 93 DMCoarsen_SNESVI - Computes the regular coarsened DM then computes additional information about its inactive set 94 95 */ 96 PetscErrorCode DMCoarsen_SNESVI(DM dm1,MPI_Comm comm,DM *dm2) 97 { 98 PetscContainer isnes; 99 DM_SNESVI *dmsnesvi1; 100 Vec finemarked,coarsemarked; 101 IS inactive; 102 Mat inject; 103 const PetscInt *index; 104 PetscInt n,k,cnt = 0,rstart,*coarseindex; 105 PetscScalar *marked; 106 107 PetscFunctionBegin; 108 CHKERRQ(PetscObjectQuery((PetscObject)dm1,"VI",(PetscObject*)&isnes)); 109 PetscCheck(isnes,PetscObjectComm((PetscObject)dm1),PETSC_ERR_PLIB,"Composed VI data structure is missing"); 110 CHKERRQ(PetscContainerGetPointer(isnes,(void**)&dmsnesvi1)); 111 112 /* get the original coarsen */ 113 CHKERRQ((*dmsnesvi1->coarsen)(dm1,comm,dm2)); 114 115 /* not sure why this extra reference is needed, but without the dm2 disappears too early */ 116 /* Updating the KSPCreateVecs() to avoid using DMGetGlobalVector() when matrix is available removes the need for this reference? */ 117 /* CHKERRQ(PetscObjectReference((PetscObject)*dm2));*/ 118 119 /* need to set back global vectors in order to use the original injection */ 120 CHKERRQ(DMClearGlobalVectors(dm1)); 121 122 dm1->ops->createglobalvector = dmsnesvi1->createglobalvector; 123 124 CHKERRQ(DMCreateGlobalVector(dm1,&finemarked)); 125 CHKERRQ(DMCreateGlobalVector(*dm2,&coarsemarked)); 126 127 /* 128 fill finemarked with locations of inactive points 129 */ 130 CHKERRQ(ISGetIndices(dmsnesvi1->inactive,&index)); 131 CHKERRQ(ISGetLocalSize(dmsnesvi1->inactive,&n)); 132 CHKERRQ(VecSet(finemarked,0.0)); 133 for (k=0; k<n; k++) { 134 CHKERRQ(VecSetValue(finemarked,index[k],1.0,INSERT_VALUES)); 135 } 136 CHKERRQ(VecAssemblyBegin(finemarked)); 137 CHKERRQ(VecAssemblyEnd(finemarked)); 138 139 CHKERRQ(DMCreateInjection(*dm2,dm1,&inject)); 140 CHKERRQ(MatRestrict(inject,finemarked,coarsemarked)); 141 CHKERRQ(MatDestroy(&inject)); 142 143 /* 144 create index set list of coarse inactive points from coarsemarked 145 */ 146 CHKERRQ(VecGetLocalSize(coarsemarked,&n)); 147 CHKERRQ(VecGetOwnershipRange(coarsemarked,&rstart,NULL)); 148 CHKERRQ(VecGetArray(coarsemarked,&marked)); 149 for (k=0; k<n; k++) { 150 if (marked[k] != 0.0) cnt++; 151 } 152 CHKERRQ(PetscMalloc1(cnt,&coarseindex)); 153 cnt = 0; 154 for (k=0; k<n; k++) { 155 if (marked[k] != 0.0) coarseindex[cnt++] = k + rstart; 156 } 157 CHKERRQ(VecRestoreArray(coarsemarked,&marked)); 158 CHKERRQ(ISCreateGeneral(PetscObjectComm((PetscObject)coarsemarked),cnt,coarseindex,PETSC_OWN_POINTER,&inactive)); 159 160 CHKERRQ(DMClearGlobalVectors(dm1)); 161 162 dm1->ops->createglobalvector = DMCreateGlobalVector_SNESVI; 163 164 CHKERRQ(DMSetVI(*dm2,inactive)); 165 166 CHKERRQ(VecDestroy(&finemarked)); 167 CHKERRQ(VecDestroy(&coarsemarked)); 168 CHKERRQ(ISDestroy(&inactive)); 169 PetscFunctionReturn(0); 170 } 171 172 PetscErrorCode DMDestroy_SNESVI(DM_SNESVI *dmsnesvi) 173 { 174 PetscFunctionBegin; 175 /* reset the base methods in the DM object that were changed when the DM_SNESVI was reset */ 176 dmsnesvi->dm->ops->createinterpolation = dmsnesvi->createinterpolation; 177 dmsnesvi->dm->ops->coarsen = dmsnesvi->coarsen; 178 dmsnesvi->dm->ops->createglobalvector = dmsnesvi->createglobalvector; 179 dmsnesvi->dm->ops->createinjection = dmsnesvi->createinjection; 180 dmsnesvi->dm->ops->hascreateinjection = dmsnesvi->hascreateinjection; 181 /* need to clear out this vectors because some of them may not have a reference to the DM 182 but they are counted as having references to the DM in DMDestroy() */ 183 CHKERRQ(DMClearGlobalVectors(dmsnesvi->dm)); 184 185 CHKERRQ(ISDestroy(&dmsnesvi->inactive)); 186 CHKERRQ(PetscFree(dmsnesvi)); 187 PetscFunctionReturn(0); 188 } 189 190 /* 191 DMSetVI - Marks a DM as associated with a VI problem. This causes the interpolation/restriction operators to 192 be restricted to only those variables NOT associated with active constraints. 193 194 */ 195 static PetscErrorCode DMSetVI(DM dm,IS inactive) 196 { 197 PetscContainer isnes; 198 DM_SNESVI *dmsnesvi; 199 200 PetscFunctionBegin; 201 if (!dm) PetscFunctionReturn(0); 202 203 CHKERRQ(PetscObjectReference((PetscObject)inactive)); 204 205 CHKERRQ(PetscObjectQuery((PetscObject)dm,"VI",(PetscObject*)&isnes)); 206 if (!isnes) { 207 CHKERRQ(PetscContainerCreate(PetscObjectComm((PetscObject)dm),&isnes)); 208 CHKERRQ(PetscContainerSetUserDestroy(isnes,(PetscErrorCode (*)(void*))DMDestroy_SNESVI)); 209 CHKERRQ(PetscNew(&dmsnesvi)); 210 CHKERRQ(PetscContainerSetPointer(isnes,(void*)dmsnesvi)); 211 CHKERRQ(PetscObjectCompose((PetscObject)dm,"VI",(PetscObject)isnes)); 212 CHKERRQ(PetscContainerDestroy(&isnes)); 213 214 dmsnesvi->createinterpolation = dm->ops->createinterpolation; 215 dm->ops->createinterpolation = DMCreateInterpolation_SNESVI; 216 dmsnesvi->coarsen = dm->ops->coarsen; 217 dm->ops->coarsen = DMCoarsen_SNESVI; 218 dmsnesvi->createglobalvector = dm->ops->createglobalvector; 219 dm->ops->createglobalvector = DMCreateGlobalVector_SNESVI; 220 dmsnesvi->createinjection = dm->ops->createinjection; 221 dm->ops->createinjection = NULL; 222 dmsnesvi->hascreateinjection = dm->ops->hascreateinjection; 223 dm->ops->hascreateinjection = NULL; 224 } else { 225 CHKERRQ(PetscContainerGetPointer(isnes,(void**)&dmsnesvi)); 226 CHKERRQ(ISDestroy(&dmsnesvi->inactive)); 227 } 228 CHKERRQ(DMClearGlobalVectors(dm)); 229 CHKERRQ(ISGetLocalSize(inactive,&dmsnesvi->n)); 230 231 dmsnesvi->inactive = inactive; 232 dmsnesvi->dm = dm; 233 PetscFunctionReturn(0); 234 } 235 236 /* 237 DMDestroyVI - Frees the DM_SNESVI object contained in the DM 238 - also resets the function pointers in the DM for createinterpolation() etc to use the original DM 239 */ 240 static PetscErrorCode DMDestroyVI(DM dm) 241 { 242 PetscFunctionBegin; 243 if (!dm) PetscFunctionReturn(0); 244 CHKERRQ(PetscObjectCompose((PetscObject)dm,"VI",(PetscObject)NULL)); 245 PetscFunctionReturn(0); 246 } 247 248 /* --------------------------------------------------------------------------------------------------------*/ 249 250 PetscErrorCode SNESCreateIndexSets_VINEWTONRSLS(SNES snes,Vec X,Vec F,IS *ISact,IS *ISinact) 251 { 252 PetscFunctionBegin; 253 CHKERRQ(SNESVIGetActiveSetIS(snes,X,F,ISact)); 254 CHKERRQ(ISComplement(*ISact,X->map->rstart,X->map->rend,ISinact)); 255 PetscFunctionReturn(0); 256 } 257 258 /* Create active and inactive set vectors. The local size of this vector is set and petsc computes the global size */ 259 PetscErrorCode SNESCreateSubVectors_VINEWTONRSLS(SNES snes,PetscInt n,Vec *newv) 260 { 261 Vec v; 262 263 PetscFunctionBegin; 264 CHKERRQ(VecCreate(PetscObjectComm((PetscObject)snes),&v)); 265 CHKERRQ(VecSetSizes(v,n,PETSC_DECIDE)); 266 CHKERRQ(VecSetType(v,VECSTANDARD)); 267 *newv = v; 268 PetscFunctionReturn(0); 269 } 270 271 /* Resets the snes PC and KSP when the active set sizes change */ 272 PetscErrorCode SNESVIResetPCandKSP(SNES snes,Mat Amat,Mat Pmat) 273 { 274 KSP snesksp; 275 276 PetscFunctionBegin; 277 CHKERRQ(SNESGetKSP(snes,&snesksp)); 278 CHKERRQ(KSPReset(snesksp)); 279 CHKERRQ(KSPResetFromOptions(snesksp)); 280 281 /* 282 KSP kspnew; 283 PC pcnew; 284 MatSolverType stype; 285 286 CHKERRQ(KSPCreate(PetscObjectComm((PetscObject)snes),&kspnew)); 287 kspnew->pc_side = snesksp->pc_side; 288 kspnew->rtol = snesksp->rtol; 289 kspnew->abstol = snesksp->abstol; 290 kspnew->max_it = snesksp->max_it; 291 CHKERRQ(KSPSetType(kspnew,((PetscObject)snesksp)->type_name)); 292 CHKERRQ(KSPGetPC(kspnew,&pcnew)); 293 CHKERRQ(PCSetType(kspnew->pc,((PetscObject)snesksp->pc)->type_name)); 294 CHKERRQ(PCSetOperators(kspnew->pc,Amat,Pmat)); 295 CHKERRQ(PCFactorGetMatSolverType(snesksp->pc,&stype)); 296 CHKERRQ(PCFactorSetMatSolverType(kspnew->pc,stype)); 297 CHKERRQ(KSPDestroy(&snesksp)); 298 snes->ksp = kspnew; 299 CHKERRQ(PetscLogObjectParent((PetscObject)snes,(PetscObject)kspnew)); 300 CHKERRQ(KSPSetFromOptions(kspnew));*/ 301 PetscFunctionReturn(0); 302 } 303 304 /* Variational Inequality solver using reduce space method. No semismooth algorithm is 305 implemented in this algorithm. It basically identifies the active constraints and does 306 a linear solve on the other variables (those not associated with the active constraints). */ 307 PetscErrorCode SNESSolve_VINEWTONRSLS(SNES snes) 308 { 309 SNES_VINEWTONRSLS *vi = (SNES_VINEWTONRSLS*)snes->data; 310 PetscInt maxits,i,lits; 311 SNESLineSearchReason lssucceed; 312 PetscReal fnorm,gnorm,xnorm=0,ynorm; 313 Vec Y,X,F; 314 KSPConvergedReason kspreason; 315 KSP ksp; 316 PC pc; 317 318 PetscFunctionBegin; 319 /* Multigrid must use Galerkin for coarse grids with active set/reduced space methods; cannot rediscretize on coarser grids*/ 320 CHKERRQ(SNESGetKSP(snes,&ksp)); 321 CHKERRQ(KSPGetPC(ksp,&pc)); 322 CHKERRQ(PCMGSetGalerkin(pc,PC_MG_GALERKIN_BOTH)); 323 324 snes->numFailures = 0; 325 snes->numLinearSolveFailures = 0; 326 snes->reason = SNES_CONVERGED_ITERATING; 327 328 maxits = snes->max_its; /* maximum number of iterations */ 329 X = snes->vec_sol; /* solution vector */ 330 F = snes->vec_func; /* residual vector */ 331 Y = snes->work[0]; /* work vectors */ 332 333 CHKERRQ(SNESLineSearchSetVIFunctions(snes->linesearch, SNESVIProjectOntoBounds, SNESVIComputeInactiveSetFnorm)); 334 CHKERRQ(SNESLineSearchSetVecs(snes->linesearch, X, NULL, NULL, NULL, NULL)); 335 CHKERRQ(SNESLineSearchSetUp(snes->linesearch)); 336 337 CHKERRQ(PetscObjectSAWsTakeAccess((PetscObject)snes)); 338 snes->iter = 0; 339 snes->norm = 0.0; 340 CHKERRQ(PetscObjectSAWsGrantAccess((PetscObject)snes)); 341 342 CHKERRQ(SNESVIProjectOntoBounds(snes,X)); 343 CHKERRQ(SNESComputeFunction(snes,X,F)); 344 CHKERRQ(SNESVIComputeInactiveSetFnorm(snes,F,X,&fnorm)); 345 CHKERRQ(VecNorm(X,NORM_2,&xnorm)); /* xnorm <- ||x|| */ 346 SNESCheckFunctionNorm(snes,fnorm); 347 CHKERRQ(PetscObjectSAWsTakeAccess((PetscObject)snes)); 348 snes->norm = fnorm; 349 CHKERRQ(PetscObjectSAWsGrantAccess((PetscObject)snes)); 350 CHKERRQ(SNESLogConvergenceHistory(snes,fnorm,0)); 351 CHKERRQ(SNESMonitor(snes,0,fnorm)); 352 353 /* test convergence */ 354 CHKERRQ((*snes->ops->converged)(snes,0,0.0,0.0,fnorm,&snes->reason,snes->cnvP)); 355 if (snes->reason) PetscFunctionReturn(0); 356 357 for (i=0; i<maxits; i++) { 358 359 IS IS_act; /* _act -> active set _inact -> inactive set */ 360 IS IS_redact; /* redundant active set */ 361 VecScatter scat_act,scat_inact; 362 PetscInt nis_act,nis_inact; 363 Vec Y_act,Y_inact,F_inact; 364 Mat jac_inact_inact,prejac_inact_inact; 365 PetscBool isequal; 366 367 /* Call general purpose update function */ 368 if (snes->ops->update) { 369 CHKERRQ((*snes->ops->update)(snes, snes->iter)); 370 } 371 CHKERRQ(SNESComputeJacobian(snes,X,snes->jacobian,snes->jacobian_pre)); 372 SNESCheckJacobianDomainerror(snes); 373 374 /* Create active and inactive index sets */ 375 376 /*original 377 CHKERRQ(SNESVICreateIndexSets_RS(snes,X,F,&IS_act,&vi->IS_inact)); 378 */ 379 CHKERRQ(SNESVIGetActiveSetIS(snes,X,F,&IS_act)); 380 381 if (vi->checkredundancy) { 382 CHKERRQ((*vi->checkredundancy)(snes,IS_act,&IS_redact,vi->ctxP)); 383 if (IS_redact) { 384 CHKERRQ(ISSort(IS_redact)); 385 CHKERRQ(ISComplement(IS_redact,X->map->rstart,X->map->rend,&vi->IS_inact)); 386 CHKERRQ(ISDestroy(&IS_redact)); 387 } else { 388 CHKERRQ(ISComplement(IS_act,X->map->rstart,X->map->rend,&vi->IS_inact)); 389 } 390 } else { 391 CHKERRQ(ISComplement(IS_act,X->map->rstart,X->map->rend,&vi->IS_inact)); 392 } 393 394 /* Create inactive set submatrix */ 395 CHKERRQ(MatCreateSubMatrix(snes->jacobian,vi->IS_inact,vi->IS_inact,MAT_INITIAL_MATRIX,&jac_inact_inact)); 396 397 if (0) { /* Dead code (temporary developer hack) */ 398 IS keptrows; 399 CHKERRQ(MatFindNonzeroRows(jac_inact_inact,&keptrows)); 400 if (keptrows) { 401 PetscInt cnt,*nrows,k; 402 const PetscInt *krows,*inact; 403 PetscInt rstart; 404 405 CHKERRQ(MatGetOwnershipRange(jac_inact_inact,&rstart,NULL)); 406 CHKERRQ(MatDestroy(&jac_inact_inact)); 407 CHKERRQ(ISDestroy(&IS_act)); 408 409 CHKERRQ(ISGetLocalSize(keptrows,&cnt)); 410 CHKERRQ(ISGetIndices(keptrows,&krows)); 411 CHKERRQ(ISGetIndices(vi->IS_inact,&inact)); 412 CHKERRQ(PetscMalloc1(cnt,&nrows)); 413 for (k=0; k<cnt; k++) nrows[k] = inact[krows[k]-rstart]; 414 CHKERRQ(ISRestoreIndices(keptrows,&krows)); 415 CHKERRQ(ISRestoreIndices(vi->IS_inact,&inact)); 416 CHKERRQ(ISDestroy(&keptrows)); 417 CHKERRQ(ISDestroy(&vi->IS_inact)); 418 419 CHKERRQ(ISCreateGeneral(PetscObjectComm((PetscObject)snes),cnt,nrows,PETSC_OWN_POINTER,&vi->IS_inact)); 420 CHKERRQ(ISComplement(vi->IS_inact,F->map->rstart,F->map->rend,&IS_act)); 421 CHKERRQ(MatCreateSubMatrix(snes->jacobian,vi->IS_inact,vi->IS_inact,MAT_INITIAL_MATRIX,&jac_inact_inact)); 422 } 423 } 424 CHKERRQ(DMSetVI(snes->dm,vi->IS_inact)); 425 /* remove later */ 426 427 /* 428 CHKERRQ(VecView(vi->xu,PETSC_VIEWER_BINARY_(((PetscObject)(vi->xu))->comm))); 429 CHKERRQ(VecView(vi->xl,PETSC_VIEWER_BINARY_(((PetscObject)(vi->xl))->comm))); 430 CHKERRQ(VecView(X,PETSC_VIEWER_BINARY_(PetscObjectComm((PetscObject)X)))); 431 CHKERRQ(VecView(F,PETSC_VIEWER_BINARY_(PetscObjectComm((PetscObject)F)))); 432 CHKERRQ(ISView(vi->IS_inact,PETSC_VIEWER_BINARY_(PetscObjectComm((PetscObject)vi->IS_inact)))); 433 */ 434 435 /* Get sizes of active and inactive sets */ 436 CHKERRQ(ISGetLocalSize(IS_act,&nis_act)); 437 CHKERRQ(ISGetLocalSize(vi->IS_inact,&nis_inact)); 438 439 /* Create active and inactive set vectors */ 440 CHKERRQ(SNESCreateSubVectors_VINEWTONRSLS(snes,nis_inact,&F_inact)); 441 CHKERRQ(SNESCreateSubVectors_VINEWTONRSLS(snes,nis_act,&Y_act)); 442 CHKERRQ(SNESCreateSubVectors_VINEWTONRSLS(snes,nis_inact,&Y_inact)); 443 444 /* Create scatter contexts */ 445 CHKERRQ(VecScatterCreate(Y,IS_act,Y_act,NULL,&scat_act)); 446 CHKERRQ(VecScatterCreate(Y,vi->IS_inact,Y_inact,NULL,&scat_inact)); 447 448 /* Do a vec scatter to active and inactive set vectors */ 449 CHKERRQ(VecScatterBegin(scat_inact,F,F_inact,INSERT_VALUES,SCATTER_FORWARD)); 450 CHKERRQ(VecScatterEnd(scat_inact,F,F_inact,INSERT_VALUES,SCATTER_FORWARD)); 451 452 CHKERRQ(VecScatterBegin(scat_act,Y,Y_act,INSERT_VALUES,SCATTER_FORWARD)); 453 CHKERRQ(VecScatterEnd(scat_act,Y,Y_act,INSERT_VALUES,SCATTER_FORWARD)); 454 455 CHKERRQ(VecScatterBegin(scat_inact,Y,Y_inact,INSERT_VALUES,SCATTER_FORWARD)); 456 CHKERRQ(VecScatterEnd(scat_inact,Y,Y_inact,INSERT_VALUES,SCATTER_FORWARD)); 457 458 /* Active set direction = 0 */ 459 CHKERRQ(VecSet(Y_act,0)); 460 if (snes->jacobian != snes->jacobian_pre) { 461 CHKERRQ(MatCreateSubMatrix(snes->jacobian_pre,vi->IS_inact,vi->IS_inact,MAT_INITIAL_MATRIX,&prejac_inact_inact)); 462 } else prejac_inact_inact = jac_inact_inact; 463 464 CHKERRQ(ISEqual(vi->IS_inact_prev,vi->IS_inact,&isequal)); 465 if (!isequal) { 466 CHKERRQ(SNESVIResetPCandKSP(snes,jac_inact_inact,prejac_inact_inact)); 467 CHKERRQ(PCFieldSplitRestrictIS(pc,vi->IS_inact)); 468 } 469 470 /* CHKERRQ(ISView(vi->IS_inact,0)); */ 471 /* CHKERRQ(ISView(IS_act,0));*/ 472 /* ierr = MatView(snes->jacobian_pre,0); */ 473 474 CHKERRQ(KSPSetOperators(snes->ksp,jac_inact_inact,prejac_inact_inact)); 475 CHKERRQ(KSPSetUp(snes->ksp)); 476 { 477 PC pc; 478 PetscBool flg; 479 CHKERRQ(KSPGetPC(snes->ksp,&pc)); 480 CHKERRQ(PetscObjectTypeCompare((PetscObject)pc,PCFIELDSPLIT,&flg)); 481 if (flg) { 482 KSP *subksps; 483 CHKERRQ(PCFieldSplitGetSubKSP(pc,NULL,&subksps)); 484 CHKERRQ(KSPGetPC(subksps[0],&pc)); 485 CHKERRQ(PetscFree(subksps)); 486 CHKERRQ(PetscObjectTypeCompare((PetscObject)pc,PCBJACOBI,&flg)); 487 if (flg) { 488 PetscInt n,N = 101*101,j,cnts[3] = {0,0,0}; 489 const PetscInt *ii; 490 491 CHKERRQ(ISGetSize(vi->IS_inact,&n)); 492 CHKERRQ(ISGetIndices(vi->IS_inact,&ii)); 493 for (j=0; j<n; j++) { 494 if (ii[j] < N) cnts[0]++; 495 else if (ii[j] < 2*N) cnts[1]++; 496 else if (ii[j] < 3*N) cnts[2]++; 497 } 498 CHKERRQ(ISRestoreIndices(vi->IS_inact,&ii)); 499 500 CHKERRQ(PCBJacobiSetTotalBlocks(pc,3,cnts)); 501 } 502 } 503 } 504 505 CHKERRQ(KSPSolve(snes->ksp,F_inact,Y_inact)); 506 CHKERRQ(VecScatterBegin(scat_act,Y_act,Y,INSERT_VALUES,SCATTER_REVERSE)); 507 CHKERRQ(VecScatterEnd(scat_act,Y_act,Y,INSERT_VALUES,SCATTER_REVERSE)); 508 CHKERRQ(VecScatterBegin(scat_inact,Y_inact,Y,INSERT_VALUES,SCATTER_REVERSE)); 509 CHKERRQ(VecScatterEnd(scat_inact,Y_inact,Y,INSERT_VALUES,SCATTER_REVERSE)); 510 511 CHKERRQ(VecDestroy(&F_inact)); 512 CHKERRQ(VecDestroy(&Y_act)); 513 CHKERRQ(VecDestroy(&Y_inact)); 514 CHKERRQ(VecScatterDestroy(&scat_act)); 515 CHKERRQ(VecScatterDestroy(&scat_inact)); 516 CHKERRQ(ISDestroy(&IS_act)); 517 if (!isequal) { 518 CHKERRQ(ISDestroy(&vi->IS_inact_prev)); 519 CHKERRQ(ISDuplicate(vi->IS_inact,&vi->IS_inact_prev)); 520 } 521 CHKERRQ(ISDestroy(&vi->IS_inact)); 522 CHKERRQ(MatDestroy(&jac_inact_inact)); 523 if (snes->jacobian != snes->jacobian_pre) { 524 CHKERRQ(MatDestroy(&prejac_inact_inact)); 525 } 526 527 CHKERRQ(KSPGetConvergedReason(snes->ksp,&kspreason)); 528 if (kspreason < 0) { 529 if (++snes->numLinearSolveFailures >= snes->maxLinearSolveFailures) { 530 CHKERRQ(PetscInfo(snes,"iter=%D, number linear solve failures %D greater than current SNES allowed, stopping solve\n",snes->iter,snes->numLinearSolveFailures)); 531 snes->reason = SNES_DIVERGED_LINEAR_SOLVE; 532 break; 533 } 534 } 535 536 CHKERRQ(KSPGetIterationNumber(snes->ksp,&lits)); 537 snes->linear_its += lits; 538 CHKERRQ(PetscInfo(snes,"iter=%D, linear solve iterations=%D\n",snes->iter,lits)); 539 /* 540 if (snes->ops->precheck) { 541 PetscBool changed_y = PETSC_FALSE; 542 CHKERRQ((*snes->ops->precheck)(snes,X,Y,snes->precheck,&changed_y)); 543 } 544 545 if (PetscLogPrintInfo) { 546 CHKERRQ(SNESVICheckResidual_Private(snes,snes->jacobian,F,Y,G,W)); 547 } 548 */ 549 /* Compute a (scaled) negative update in the line search routine: 550 Y <- X - lambda*Y 551 and evaluate G = function(Y) (depends on the line search). 552 */ 553 CHKERRQ(VecCopy(Y,snes->vec_sol_update)); 554 ynorm = 1; gnorm = fnorm; 555 CHKERRQ(SNESLineSearchApply(snes->linesearch, X, F, &gnorm, Y)); 556 CHKERRQ(SNESLineSearchGetReason(snes->linesearch, &lssucceed)); 557 CHKERRQ(SNESLineSearchGetNorms(snes->linesearch, &xnorm, &gnorm, &ynorm)); 558 CHKERRQ(PetscInfo(snes,"fnorm=%18.16e, gnorm=%18.16e, ynorm=%18.16e, lssucceed=%d\n",(double)fnorm,(double)gnorm,(double)ynorm,(int)lssucceed)); 559 if (snes->reason == SNES_DIVERGED_FUNCTION_COUNT) break; 560 if (snes->domainerror) { 561 snes->reason = SNES_DIVERGED_FUNCTION_DOMAIN; 562 CHKERRQ(DMDestroyVI(snes->dm)); 563 PetscFunctionReturn(0); 564 } 565 if (lssucceed) { 566 if (++snes->numFailures >= snes->maxFailures) { 567 PetscBool ismin; 568 snes->reason = SNES_DIVERGED_LINE_SEARCH; 569 CHKERRQ(SNESVICheckLocalMin_Private(snes,snes->jacobian,F,X,gnorm,&ismin)); 570 if (ismin) snes->reason = SNES_DIVERGED_LOCAL_MIN; 571 break; 572 } 573 } 574 CHKERRQ(DMDestroyVI(snes->dm)); 575 /* Update function and solution vectors */ 576 fnorm = gnorm; 577 /* Monitor convergence */ 578 CHKERRQ(PetscObjectSAWsTakeAccess((PetscObject)snes)); 579 snes->iter = i+1; 580 snes->norm = fnorm; 581 snes->xnorm = xnorm; 582 snes->ynorm = ynorm; 583 CHKERRQ(PetscObjectSAWsGrantAccess((PetscObject)snes)); 584 CHKERRQ(SNESLogConvergenceHistory(snes,snes->norm,lits)); 585 CHKERRQ(SNESMonitor(snes,snes->iter,snes->norm)); 586 /* Test for convergence, xnorm = || X || */ 587 if (snes->ops->converged != SNESConvergedSkip) CHKERRQ(VecNorm(X,NORM_2,&xnorm)); 588 CHKERRQ((*snes->ops->converged)(snes,snes->iter,xnorm,ynorm,fnorm,&snes->reason,snes->cnvP)); 589 if (snes->reason) break; 590 } 591 /* make sure that the VI information attached to the DM is removed if the for loop above was broken early due to some exceptional conditional */ 592 CHKERRQ(DMDestroyVI(snes->dm)); 593 if (i == maxits) { 594 CHKERRQ(PetscInfo(snes,"Maximum number of iterations has been reached: %D\n",maxits)); 595 if (!snes->reason) snes->reason = SNES_DIVERGED_MAX_IT; 596 } 597 PetscFunctionReturn(0); 598 } 599 600 PetscErrorCode SNESVISetRedundancyCheck(SNES snes,PetscErrorCode (*func)(SNES,IS,IS*,void*),void *ctx) 601 { 602 SNES_VINEWTONRSLS *vi = (SNES_VINEWTONRSLS*)snes->data; 603 604 PetscFunctionBegin; 605 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 606 vi->checkredundancy = func; 607 vi->ctxP = ctx; 608 PetscFunctionReturn(0); 609 } 610 611 #if defined(PETSC_HAVE_MATLAB_ENGINE) 612 #include <engine.h> 613 #include <mex.h> 614 typedef struct {char *funcname; mxArray *ctx;} SNESMatlabContext; 615 616 PetscErrorCode SNESVIRedundancyCheck_Matlab(SNES snes,IS is_act,IS *is_redact,void *ctx) 617 { 618 SNESMatlabContext *sctx = (SNESMatlabContext*)ctx; 619 int nlhs = 1, nrhs = 5; 620 mxArray *plhs[1], *prhs[5]; 621 long long int l1 = 0, l2 = 0, ls = 0; 622 PetscInt *indices=NULL; 623 624 PetscFunctionBegin; 625 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 626 PetscValidHeaderSpecific(is_act,IS_CLASSID,2); 627 PetscValidPointer(is_redact,3); 628 PetscCheckSameComm(snes,1,is_act,2); 629 630 /* Create IS for reduced active set of size 0, its size and indices will 631 bet set by the Matlab function */ 632 CHKERRQ(ISCreateGeneral(PetscObjectComm((PetscObject)snes),0,indices,PETSC_OWN_POINTER,is_redact)); 633 /* call Matlab function in ctx */ 634 CHKERRQ(PetscArraycpy(&ls,&snes,1)); 635 CHKERRQ(PetscArraycpy(&l1,&is_act,1)); 636 CHKERRQ(PetscArraycpy(&l2,is_redact,1)); 637 prhs[0] = mxCreateDoubleScalar((double)ls); 638 prhs[1] = mxCreateDoubleScalar((double)l1); 639 prhs[2] = mxCreateDoubleScalar((double)l2); 640 prhs[3] = mxCreateString(sctx->funcname); 641 prhs[4] = sctx->ctx; 642 CHKERRQ(mexCallMATLAB(nlhs,plhs,nrhs,prhs,"PetscSNESVIRedundancyCheckInternal")); 643 CHKERRQ(mxGetScalar(plhs[0])); 644 mxDestroyArray(prhs[0]); 645 mxDestroyArray(prhs[1]); 646 mxDestroyArray(prhs[2]); 647 mxDestroyArray(prhs[3]); 648 mxDestroyArray(plhs[0]); 649 PetscFunctionReturn(0); 650 } 651 652 PetscErrorCode SNESVISetRedundancyCheckMatlab(SNES snes,const char *func,mxArray *ctx) 653 { 654 SNESMatlabContext *sctx; 655 656 PetscFunctionBegin; 657 /* currently sctx is memory bleed */ 658 CHKERRQ(PetscNew(&sctx)); 659 CHKERRQ(PetscStrallocpy(func,&sctx->funcname)); 660 sctx->ctx = mxDuplicateArray(ctx); 661 CHKERRQ(SNESVISetRedundancyCheck(snes,SNESVIRedundancyCheck_Matlab,sctx)); 662 PetscFunctionReturn(0); 663 } 664 665 #endif 666 667 /* -------------------------------------------------------------------------- */ 668 /* 669 SNESSetUp_VINEWTONRSLS - Sets up the internal data structures for the later use 670 of the SNESVI nonlinear solver. 671 672 Input Parameter: 673 . snes - the SNES context 674 675 Application Interface Routine: SNESSetUp() 676 677 Notes: 678 For basic use of the SNES solvers, the user need not explicitly call 679 SNESSetUp(), since these actions will automatically occur during 680 the call to SNESSolve(). 681 */ 682 PetscErrorCode SNESSetUp_VINEWTONRSLS(SNES snes) 683 { 684 SNES_VINEWTONRSLS *vi = (SNES_VINEWTONRSLS*) snes->data; 685 PetscInt *indices; 686 PetscInt i,n,rstart,rend; 687 SNESLineSearch linesearch; 688 689 PetscFunctionBegin; 690 CHKERRQ(SNESSetUp_VI(snes)); 691 692 /* Set up previous active index set for the first snes solve 693 vi->IS_inact_prev = 0,1,2,....N */ 694 695 CHKERRQ(VecGetOwnershipRange(snes->vec_sol,&rstart,&rend)); 696 CHKERRQ(VecGetLocalSize(snes->vec_sol,&n)); 697 CHKERRQ(PetscMalloc1(n,&indices)); 698 for (i=0; i < n; i++) indices[i] = rstart + i; 699 CHKERRQ(ISCreateGeneral(PetscObjectComm((PetscObject)snes),n,indices,PETSC_OWN_POINTER,&vi->IS_inact_prev)); 700 701 /* set the line search functions */ 702 if (!snes->linesearch) { 703 CHKERRQ(SNESGetLineSearch(snes, &linesearch)); 704 CHKERRQ(SNESLineSearchSetType(linesearch, SNESLINESEARCHBT)); 705 } 706 PetscFunctionReturn(0); 707 } 708 /* -------------------------------------------------------------------------- */ 709 PetscErrorCode SNESReset_VINEWTONRSLS(SNES snes) 710 { 711 SNES_VINEWTONRSLS *vi = (SNES_VINEWTONRSLS*) snes->data; 712 713 PetscFunctionBegin; 714 CHKERRQ(SNESReset_VI(snes)); 715 CHKERRQ(ISDestroy(&vi->IS_inact_prev)); 716 PetscFunctionReturn(0); 717 } 718 719 /* -------------------------------------------------------------------------- */ 720 /*MC 721 SNESVINEWTONRSLS - Reduced space active set solvers for variational inequalities based on Newton's method 722 723 Options Database: 724 + -snes_type <vinewtonssls,vinewtonrsls> - a semi-smooth solver, a reduced space active set method 725 - -snes_vi_monitor - prints the number of active constraints at each iteration. 726 727 Level: beginner 728 729 References: 730 . * - T. S. Munson, and S. Benson. Flexible Complementarity Solvers for Large Scale 731 Applications, Optimization Methods and Software, 21 (2006). 732 733 .seealso: SNESVISetVariableBounds(), SNESVISetComputeVariableBounds(), SNESCreate(), SNES, SNESSetType(), SNESVINEWTONSSLS, SNESNEWTONTR, SNESLineSearchSetType(),SNESLineSearchSetPostCheck(), SNESLineSearchSetPreCheck() 734 735 M*/ 736 PETSC_EXTERN PetscErrorCode SNESCreate_VINEWTONRSLS(SNES snes) 737 { 738 SNES_VINEWTONRSLS *vi; 739 SNESLineSearch linesearch; 740 741 PetscFunctionBegin; 742 snes->ops->reset = SNESReset_VINEWTONRSLS; 743 snes->ops->setup = SNESSetUp_VINEWTONRSLS; 744 snes->ops->solve = SNESSolve_VINEWTONRSLS; 745 snes->ops->destroy = SNESDestroy_VI; 746 snes->ops->setfromoptions = SNESSetFromOptions_VI; 747 snes->ops->view = NULL; 748 snes->ops->converged = SNESConvergedDefault_VI; 749 750 snes->usesksp = PETSC_TRUE; 751 snes->usesnpc = PETSC_FALSE; 752 753 CHKERRQ(SNESGetLineSearch(snes, &linesearch)); 754 if (!((PetscObject)linesearch)->type_name) { 755 CHKERRQ(SNESLineSearchSetType(linesearch, SNESLINESEARCHBT)); 756 } 757 CHKERRQ(SNESLineSearchBTSetAlpha(linesearch, 0.0)); 758 759 snes->alwayscomputesfinalresidual = PETSC_TRUE; 760 761 CHKERRQ(PetscNewLog(snes,&vi)); 762 snes->data = (void*)vi; 763 vi->checkredundancy = NULL; 764 765 CHKERRQ(PetscObjectComposeFunction((PetscObject)snes,"SNESVISetVariableBounds_C",SNESVISetVariableBounds_VI)); 766 CHKERRQ(PetscObjectComposeFunction((PetscObject)snes,"SNESVISetComputeVariableBounds_C",SNESVISetComputeVariableBounds_VI)); 767 PetscFunctionReturn(0); 768 } 769