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