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