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