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 ierr = KSPResetFromOptions(snesksp);CHKERRQ(ierr); 301 302 /* 303 KSP kspnew; 304 PC pcnew; 305 MatSolverType stype; 306 307 308 ierr = KSPCreate(PetscObjectComm((PetscObject)snes),&kspnew);CHKERRQ(ierr); 309 kspnew->pc_side = snesksp->pc_side; 310 kspnew->rtol = snesksp->rtol; 311 kspnew->abstol = snesksp->abstol; 312 kspnew->max_it = snesksp->max_it; 313 ierr = KSPSetType(kspnew,((PetscObject)snesksp)->type_name);CHKERRQ(ierr); 314 ierr = KSPGetPC(kspnew,&pcnew);CHKERRQ(ierr); 315 ierr = PCSetType(kspnew->pc,((PetscObject)snesksp->pc)->type_name);CHKERRQ(ierr); 316 ierr = PCSetOperators(kspnew->pc,Amat,Pmat);CHKERRQ(ierr); 317 ierr = PCFactorGetMatSolverType(snesksp->pc,&stype);CHKERRQ(ierr); 318 ierr = PCFactorSetMatSolverType(kspnew->pc,stype);CHKERRQ(ierr); 319 ierr = KSPDestroy(&snesksp);CHKERRQ(ierr); 320 snes->ksp = kspnew; 321 ierr = PetscLogObjectParent((PetscObject)snes,(PetscObject)kspnew);CHKERRQ(ierr); 322 ierr = KSPSetFromOptions(kspnew);CHKERRQ(ierr);*/ 323 PetscFunctionReturn(0); 324 } 325 326 /* Variational Inequality solver using reduce space method. No semismooth algorithm is 327 implemented in this algorithm. It basically identifies the active constraints and does 328 a linear solve on the other variables (those not associated with the active constraints). */ 329 PetscErrorCode SNESSolve_VINEWTONRSLS(SNES snes) 330 { 331 SNES_VINEWTONRSLS *vi = (SNES_VINEWTONRSLS*)snes->data; 332 PetscErrorCode ierr; 333 PetscInt maxits,i,lits; 334 SNESLineSearchReason lssucceed; 335 PetscReal fnorm,gnorm,xnorm=0,ynorm; 336 Vec Y,X,F; 337 KSPConvergedReason kspreason; 338 KSP ksp; 339 PC pc; 340 341 PetscFunctionBegin; 342 /* Multigrid must use Galerkin for coarse grids with active set/reduced space methods; cannot rediscretize on coarser grids*/ 343 ierr = SNESGetKSP(snes,&ksp);CHKERRQ(ierr); 344 ierr = KSPGetPC(ksp,&pc);CHKERRQ(ierr); 345 ierr = PCMGSetGalerkin(pc,PC_MG_GALERKIN_BOTH);CHKERRQ(ierr); 346 347 snes->numFailures = 0; 348 snes->numLinearSolveFailures = 0; 349 snes->reason = SNES_CONVERGED_ITERATING; 350 351 maxits = snes->max_its; /* maximum number of iterations */ 352 X = snes->vec_sol; /* solution vector */ 353 F = snes->vec_func; /* residual vector */ 354 Y = snes->work[0]; /* work vectors */ 355 356 ierr = SNESLineSearchSetVIFunctions(snes->linesearch, SNESVIProjectOntoBounds, SNESVIComputeInactiveSetFnorm);CHKERRQ(ierr); 357 ierr = SNESLineSearchSetVecs(snes->linesearch, X, NULL, NULL, NULL, NULL);CHKERRQ(ierr); 358 ierr = SNESLineSearchSetUp(snes->linesearch);CHKERRQ(ierr); 359 360 ierr = PetscObjectSAWsTakeAccess((PetscObject)snes);CHKERRQ(ierr); 361 snes->iter = 0; 362 snes->norm = 0.0; 363 ierr = PetscObjectSAWsGrantAccess((PetscObject)snes);CHKERRQ(ierr); 364 365 ierr = SNESVIProjectOntoBounds(snes,X);CHKERRQ(ierr); 366 ierr = SNESComputeFunction(snes,X,F);CHKERRQ(ierr); 367 ierr = SNESVIComputeInactiveSetFnorm(snes,F,X,&fnorm);CHKERRQ(ierr); 368 ierr = VecNorm(X,NORM_2,&xnorm);CHKERRQ(ierr); /* xnorm <- ||x|| */ 369 SNESCheckFunctionNorm(snes,fnorm); 370 ierr = PetscObjectSAWsTakeAccess((PetscObject)snes);CHKERRQ(ierr); 371 snes->norm = fnorm; 372 ierr = PetscObjectSAWsGrantAccess((PetscObject)snes);CHKERRQ(ierr); 373 ierr = SNESLogConvergenceHistory(snes,fnorm,0);CHKERRQ(ierr); 374 ierr = SNESMonitor(snes,0,fnorm);CHKERRQ(ierr); 375 376 /* test convergence */ 377 ierr = (*snes->ops->converged)(snes,0,0.0,0.0,fnorm,&snes->reason,snes->cnvP);CHKERRQ(ierr); 378 if (snes->reason) PetscFunctionReturn(0); 379 380 381 for (i=0; i<maxits; i++) { 382 383 IS IS_act; /* _act -> active set _inact -> inactive set */ 384 IS IS_redact; /* redundant active set */ 385 VecScatter scat_act,scat_inact; 386 PetscInt nis_act,nis_inact; 387 Vec Y_act,Y_inact,F_inact; 388 Mat jac_inact_inact,prejac_inact_inact; 389 PetscBool isequal; 390 391 /* Call general purpose update function */ 392 if (snes->ops->update) { 393 ierr = (*snes->ops->update)(snes, snes->iter);CHKERRQ(ierr); 394 } 395 ierr = SNESComputeJacobian(snes,X,snes->jacobian,snes->jacobian_pre);CHKERRQ(ierr); 396 397 398 /* Create active and inactive index sets */ 399 400 /*original 401 ierr = SNESVICreateIndexSets_RS(snes,X,F,&IS_act,&vi->IS_inact);CHKERRQ(ierr); 402 */ 403 ierr = SNESVIGetActiveSetIS(snes,X,F,&IS_act);CHKERRQ(ierr); 404 405 if (vi->checkredundancy) { 406 (*vi->checkredundancy)(snes,IS_act,&IS_redact,vi->ctxP);CHKERRQ(ierr); 407 if (IS_redact) { 408 ierr = ISSort(IS_redact);CHKERRQ(ierr); 409 ierr = ISComplement(IS_redact,X->map->rstart,X->map->rend,&vi->IS_inact);CHKERRQ(ierr); 410 ierr = ISDestroy(&IS_redact);CHKERRQ(ierr); 411 } else { 412 ierr = ISComplement(IS_act,X->map->rstart,X->map->rend,&vi->IS_inact);CHKERRQ(ierr); 413 } 414 } else { 415 ierr = ISComplement(IS_act,X->map->rstart,X->map->rend,&vi->IS_inact);CHKERRQ(ierr); 416 } 417 418 419 /* Create inactive set submatrix */ 420 ierr = MatCreateSubMatrix(snes->jacobian,vi->IS_inact,vi->IS_inact,MAT_INITIAL_MATRIX,&jac_inact_inact);CHKERRQ(ierr); 421 422 if (0) { /* Dead code (temporary developer hack) */ 423 IS keptrows; 424 ierr = MatFindNonzeroRows(jac_inact_inact,&keptrows);CHKERRQ(ierr); 425 if (keptrows) { 426 PetscInt cnt,*nrows,k; 427 const PetscInt *krows,*inact; 428 PetscInt rstart; 429 430 ierr = MatGetOwnershipRange(jac_inact_inact,&rstart,NULL);CHKERRQ(ierr); 431 ierr = MatDestroy(&jac_inact_inact);CHKERRQ(ierr); 432 ierr = ISDestroy(&IS_act);CHKERRQ(ierr); 433 434 ierr = ISGetLocalSize(keptrows,&cnt);CHKERRQ(ierr); 435 ierr = ISGetIndices(keptrows,&krows);CHKERRQ(ierr); 436 ierr = ISGetIndices(vi->IS_inact,&inact);CHKERRQ(ierr); 437 ierr = PetscMalloc1(cnt,&nrows);CHKERRQ(ierr); 438 for (k=0; k<cnt; k++) nrows[k] = inact[krows[k]-rstart]; 439 ierr = ISRestoreIndices(keptrows,&krows);CHKERRQ(ierr); 440 ierr = ISRestoreIndices(vi->IS_inact,&inact);CHKERRQ(ierr); 441 ierr = ISDestroy(&keptrows);CHKERRQ(ierr); 442 ierr = ISDestroy(&vi->IS_inact);CHKERRQ(ierr); 443 444 ierr = ISCreateGeneral(PetscObjectComm((PetscObject)snes),cnt,nrows,PETSC_OWN_POINTER,&vi->IS_inact);CHKERRQ(ierr); 445 ierr = ISComplement(vi->IS_inact,F->map->rstart,F->map->rend,&IS_act);CHKERRQ(ierr); 446 ierr = MatCreateSubMatrix(snes->jacobian,vi->IS_inact,vi->IS_inact,MAT_INITIAL_MATRIX,&jac_inact_inact);CHKERRQ(ierr); 447 } 448 } 449 ierr = DMSetVI(snes->dm,vi->IS_inact);CHKERRQ(ierr); 450 /* remove later */ 451 452 /* 453 ierr = VecView(vi->xu,PETSC_VIEWER_BINARY_(((PetscObject)(vi->xu))->comm));CHKERRQ(ierr); 454 ierr = VecView(vi->xl,PETSC_VIEWER_BINARY_(((PetscObject)(vi->xl))->comm));CHKERRQ(ierr); 455 ierr = VecView(X,PETSC_VIEWER_BINARY_(PetscObjectComm((PetscObject)X)));CHKERRQ(ierr); 456 ierr = VecView(F,PETSC_VIEWER_BINARY_(PetscObjectComm((PetscObject)F)));CHKERRQ(ierr); 457 ierr = ISView(vi->IS_inact,PETSC_VIEWER_BINARY_(PetscObjectComm((PetscObject)vi->IS_inact)));CHKERRQ(ierr); 458 */ 459 460 /* Get sizes of active and inactive sets */ 461 ierr = ISGetLocalSize(IS_act,&nis_act);CHKERRQ(ierr); 462 ierr = ISGetLocalSize(vi->IS_inact,&nis_inact);CHKERRQ(ierr); 463 464 /* Create active and inactive set vectors */ 465 ierr = SNESCreateSubVectors_VINEWTONRSLS(snes,nis_inact,&F_inact);CHKERRQ(ierr); 466 ierr = SNESCreateSubVectors_VINEWTONRSLS(snes,nis_act,&Y_act);CHKERRQ(ierr); 467 ierr = SNESCreateSubVectors_VINEWTONRSLS(snes,nis_inact,&Y_inact);CHKERRQ(ierr); 468 469 /* Create scatter contexts */ 470 ierr = VecScatterCreateWithData(Y,IS_act,Y_act,NULL,&scat_act);CHKERRQ(ierr); 471 ierr = VecScatterCreateWithData(Y,vi->IS_inact,Y_inact,NULL,&scat_inact);CHKERRQ(ierr); 472 473 /* Do a vec scatter to active and inactive set vectors */ 474 ierr = VecScatterBegin(scat_inact,F,F_inact,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 475 ierr = VecScatterEnd(scat_inact,F,F_inact,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 476 477 ierr = VecScatterBegin(scat_act,Y,Y_act,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 478 ierr = VecScatterEnd(scat_act,Y,Y_act,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 479 480 ierr = VecScatterBegin(scat_inact,Y,Y_inact,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 481 ierr = VecScatterEnd(scat_inact,Y,Y_inact,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 482 483 /* Active set direction = 0 */ 484 ierr = VecSet(Y_act,0);CHKERRQ(ierr); 485 if (snes->jacobian != snes->jacobian_pre) { 486 ierr = MatCreateSubMatrix(snes->jacobian_pre,vi->IS_inact,vi->IS_inact,MAT_INITIAL_MATRIX,&prejac_inact_inact);CHKERRQ(ierr); 487 } else prejac_inact_inact = jac_inact_inact; 488 489 ierr = ISEqual(vi->IS_inact_prev,vi->IS_inact,&isequal);CHKERRQ(ierr); 490 if (!isequal) { 491 ierr = SNESVIResetPCandKSP(snes,jac_inact_inact,prejac_inact_inact);CHKERRQ(ierr); 492 ierr = PCFieldSplitRestrictIS(pc,vi->IS_inact);CHKERRQ(ierr); 493 } 494 495 /* ierr = ISView(vi->IS_inact,0);CHKERRQ(ierr); */ 496 /* ierr = ISView(IS_act,0);CHKERRQ(ierr);*/ 497 /* ierr = MatView(snes->jacobian_pre,0); */ 498 499 500 501 ierr = KSPSetOperators(snes->ksp,jac_inact_inact,prejac_inact_inact);CHKERRQ(ierr); 502 ierr = KSPSetUp(snes->ksp);CHKERRQ(ierr); 503 { 504 PC pc; 505 PetscBool flg; 506 ierr = KSPGetPC(snes->ksp,&pc);CHKERRQ(ierr); 507 ierr = PetscObjectTypeCompare((PetscObject)pc,PCFIELDSPLIT,&flg);CHKERRQ(ierr); 508 if (flg) { 509 KSP *subksps; 510 ierr = PCFieldSplitGetSubKSP(pc,NULL,&subksps);CHKERRQ(ierr); 511 ierr = KSPGetPC(subksps[0],&pc);CHKERRQ(ierr); 512 ierr = PetscFree(subksps);CHKERRQ(ierr); 513 ierr = PetscObjectTypeCompare((PetscObject)pc,PCBJACOBI,&flg);CHKERRQ(ierr); 514 if (flg) { 515 PetscInt n,N = 101*101,j,cnts[3] = {0,0,0}; 516 const PetscInt *ii; 517 518 ierr = ISGetSize(vi->IS_inact,&n);CHKERRQ(ierr); 519 ierr = ISGetIndices(vi->IS_inact,&ii);CHKERRQ(ierr); 520 for (j=0; j<n; j++) { 521 if (ii[j] < N) cnts[0]++; 522 else if (ii[j] < 2*N) cnts[1]++; 523 else if (ii[j] < 3*N) cnts[2]++; 524 } 525 ierr = ISRestoreIndices(vi->IS_inact,&ii);CHKERRQ(ierr); 526 527 ierr = PCBJacobiSetTotalBlocks(pc,3,cnts);CHKERRQ(ierr); 528 } 529 } 530 } 531 532 ierr = KSPSolve(snes->ksp,F_inact,Y_inact);CHKERRQ(ierr); 533 ierr = VecScatterBegin(scat_act,Y_act,Y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 534 ierr = VecScatterEnd(scat_act,Y_act,Y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 535 ierr = VecScatterBegin(scat_inact,Y_inact,Y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 536 ierr = VecScatterEnd(scat_inact,Y_inact,Y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 537 538 ierr = VecDestroy(&F_inact);CHKERRQ(ierr); 539 ierr = VecDestroy(&Y_act);CHKERRQ(ierr); 540 ierr = VecDestroy(&Y_inact);CHKERRQ(ierr); 541 ierr = VecScatterDestroy(&scat_act);CHKERRQ(ierr); 542 ierr = VecScatterDestroy(&scat_inact);CHKERRQ(ierr); 543 ierr = ISDestroy(&IS_act);CHKERRQ(ierr); 544 if (!isequal) { 545 ierr = ISDestroy(&vi->IS_inact_prev);CHKERRQ(ierr); 546 ierr = ISDuplicate(vi->IS_inact,&vi->IS_inact_prev);CHKERRQ(ierr); 547 } 548 ierr = ISDestroy(&vi->IS_inact);CHKERRQ(ierr); 549 ierr = MatDestroy(&jac_inact_inact);CHKERRQ(ierr); 550 if (snes->jacobian != snes->jacobian_pre) { 551 ierr = MatDestroy(&prejac_inact_inact);CHKERRQ(ierr); 552 } 553 554 ierr = KSPGetConvergedReason(snes->ksp,&kspreason);CHKERRQ(ierr); 555 if (kspreason < 0) { 556 if (++snes->numLinearSolveFailures >= snes->maxLinearSolveFailures) { 557 ierr = PetscInfo2(snes,"iter=%D, number linear solve failures %D greater than current SNES allowed, stopping solve\n",snes->iter,snes->numLinearSolveFailures);CHKERRQ(ierr); 558 snes->reason = SNES_DIVERGED_LINEAR_SOLVE; 559 break; 560 } 561 } 562 563 ierr = KSPGetIterationNumber(snes->ksp,&lits);CHKERRQ(ierr); 564 snes->linear_its += lits; 565 ierr = PetscInfo2(snes,"iter=%D, linear solve iterations=%D\n",snes->iter,lits);CHKERRQ(ierr); 566 /* 567 if (snes->ops->precheck) { 568 PetscBool changed_y = PETSC_FALSE; 569 ierr = (*snes->ops->precheck)(snes,X,Y,snes->precheck,&changed_y);CHKERRQ(ierr); 570 } 571 572 if (PetscLogPrintInfo) { 573 ierr = SNESVICheckResidual_Private(snes,snes->jacobian,F,Y,G,W);CHKERRQ(ierr); 574 } 575 */ 576 /* Compute a (scaled) negative update in the line search routine: 577 Y <- X - lambda*Y 578 and evaluate G = function(Y) (depends on the line search). 579 */ 580 ierr = VecCopy(Y,snes->vec_sol_update);CHKERRQ(ierr); 581 ynorm = 1; gnorm = fnorm; 582 ierr = SNESLineSearchApply(snes->linesearch, X, F, &gnorm, Y);CHKERRQ(ierr); 583 ierr = SNESLineSearchGetReason(snes->linesearch, &lssucceed);CHKERRQ(ierr); 584 ierr = SNESLineSearchGetNorms(snes->linesearch, &xnorm, &gnorm, &ynorm);CHKERRQ(ierr); 585 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); 586 if (snes->reason == SNES_DIVERGED_FUNCTION_COUNT) break; 587 if (snes->domainerror) { 588 snes->reason = SNES_DIVERGED_FUNCTION_DOMAIN; 589 ierr = DMDestroyVI(snes->dm);CHKERRQ(ierr); 590 PetscFunctionReturn(0); 591 } 592 if (lssucceed) { 593 if (++snes->numFailures >= snes->maxFailures) { 594 PetscBool ismin; 595 snes->reason = SNES_DIVERGED_LINE_SEARCH; 596 ierr = SNESVICheckLocalMin_Private(snes,snes->jacobian,F,X,gnorm,&ismin);CHKERRQ(ierr); 597 if (ismin) snes->reason = SNES_DIVERGED_LOCAL_MIN; 598 break; 599 } 600 } 601 ierr = DMDestroyVI(snes->dm);CHKERRQ(ierr); 602 /* Update function and solution vectors */ 603 fnorm = gnorm; 604 /* Monitor convergence */ 605 ierr = PetscObjectSAWsTakeAccess((PetscObject)snes);CHKERRQ(ierr); 606 snes->iter = i+1; 607 snes->norm = fnorm; 608 ierr = PetscObjectSAWsGrantAccess((PetscObject)snes);CHKERRQ(ierr); 609 ierr = SNESLogConvergenceHistory(snes,snes->norm,lits);CHKERRQ(ierr); 610 ierr = SNESMonitor(snes,snes->iter,snes->norm);CHKERRQ(ierr); 611 /* Test for convergence, xnorm = || X || */ 612 if (snes->ops->converged != SNESConvergedSkip) { ierr = VecNorm(X,NORM_2,&xnorm);CHKERRQ(ierr); } 613 ierr = (*snes->ops->converged)(snes,snes->iter,xnorm,ynorm,fnorm,&snes->reason,snes->cnvP);CHKERRQ(ierr); 614 if (snes->reason) break; 615 } 616 /* 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 */ 617 ierr = DMDestroyVI(snes->dm);CHKERRQ(ierr); 618 if (i == maxits) { 619 ierr = PetscInfo1(snes,"Maximum number of iterations has been reached: %D\n",maxits);CHKERRQ(ierr); 620 if (!snes->reason) snes->reason = SNES_DIVERGED_MAX_IT; 621 } 622 PetscFunctionReturn(0); 623 } 624 625 PetscErrorCode SNESVISetRedundancyCheck(SNES snes,PetscErrorCode (*func)(SNES,IS,IS*,void*),void *ctx) 626 { 627 SNES_VINEWTONRSLS *vi = (SNES_VINEWTONRSLS*)snes->data; 628 629 PetscFunctionBegin; 630 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 631 vi->checkredundancy = func; 632 vi->ctxP = ctx; 633 PetscFunctionReturn(0); 634 } 635 636 #if defined(PETSC_HAVE_MATLAB_ENGINE) 637 #include <engine.h> 638 #include <mex.h> 639 typedef struct {char *funcname; mxArray *ctx;} SNESMatlabContext; 640 641 PetscErrorCode SNESVIRedundancyCheck_Matlab(SNES snes,IS is_act,IS *is_redact,void *ctx) 642 { 643 PetscErrorCode ierr; 644 SNESMatlabContext *sctx = (SNESMatlabContext*)ctx; 645 int nlhs = 1, nrhs = 5; 646 mxArray *plhs[1], *prhs[5]; 647 long long int l1 = 0, l2 = 0, ls = 0; 648 PetscInt *indices=NULL; 649 650 PetscFunctionBegin; 651 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 652 PetscValidHeaderSpecific(is_act,IS_CLASSID,2); 653 PetscValidPointer(is_redact,3); 654 PetscCheckSameComm(snes,1,is_act,2); 655 656 /* Create IS for reduced active set of size 0, its size and indices will 657 bet set by the Matlab function */ 658 ierr = ISCreateGeneral(PetscObjectComm((PetscObject)snes),0,indices,PETSC_OWN_POINTER,is_redact);CHKERRQ(ierr); 659 /* call Matlab function in ctx */ 660 ierr = PetscMemcpy(&ls,&snes,sizeof(snes));CHKERRQ(ierr); 661 ierr = PetscMemcpy(&l1,&is_act,sizeof(is_act));CHKERRQ(ierr); 662 ierr = PetscMemcpy(&l2,is_redact,sizeof(is_act));CHKERRQ(ierr); 663 prhs[0] = mxCreateDoubleScalar((double)ls); 664 prhs[1] = mxCreateDoubleScalar((double)l1); 665 prhs[2] = mxCreateDoubleScalar((double)l2); 666 prhs[3] = mxCreateString(sctx->funcname); 667 prhs[4] = sctx->ctx; 668 ierr = mexCallMATLAB(nlhs,plhs,nrhs,prhs,"PetscSNESVIRedundancyCheckInternal");CHKERRQ(ierr); 669 ierr = mxGetScalar(plhs[0]);CHKERRQ(ierr); 670 mxDestroyArray(prhs[0]); 671 mxDestroyArray(prhs[1]); 672 mxDestroyArray(prhs[2]); 673 mxDestroyArray(prhs[3]); 674 mxDestroyArray(plhs[0]); 675 PetscFunctionReturn(0); 676 } 677 678 PetscErrorCode SNESVISetRedundancyCheckMatlab(SNES snes,const char *func,mxArray *ctx) 679 { 680 PetscErrorCode ierr; 681 SNESMatlabContext *sctx; 682 683 PetscFunctionBegin; 684 /* currently sctx is memory bleed */ 685 ierr = PetscNew(&sctx);CHKERRQ(ierr); 686 ierr = PetscStrallocpy(func,&sctx->funcname);CHKERRQ(ierr); 687 sctx->ctx = mxDuplicateArray(ctx); 688 ierr = SNESVISetRedundancyCheck(snes,SNESVIRedundancyCheck_Matlab,sctx);CHKERRQ(ierr); 689 PetscFunctionReturn(0); 690 } 691 692 #endif 693 694 /* -------------------------------------------------------------------------- */ 695 /* 696 SNESSetUp_VINEWTONRSLS - Sets up the internal data structures for the later use 697 of the SNESVI nonlinear solver. 698 699 Input Parameter: 700 . snes - the SNES context 701 702 Application Interface Routine: SNESSetUp() 703 704 Notes: 705 For basic use of the SNES solvers, the user need not explicitly call 706 SNESSetUp(), since these actions will automatically occur during 707 the call to SNESSolve(). 708 */ 709 PetscErrorCode SNESSetUp_VINEWTONRSLS(SNES snes) 710 { 711 PetscErrorCode ierr; 712 SNES_VINEWTONRSLS *vi = (SNES_VINEWTONRSLS*) snes->data; 713 PetscInt *indices; 714 PetscInt i,n,rstart,rend; 715 SNESLineSearch linesearch; 716 717 PetscFunctionBegin; 718 ierr = SNESSetUp_VI(snes);CHKERRQ(ierr); 719 720 /* Set up previous active index set for the first snes solve 721 vi->IS_inact_prev = 0,1,2,....N */ 722 723 ierr = VecGetOwnershipRange(snes->vec_sol,&rstart,&rend);CHKERRQ(ierr); 724 ierr = VecGetLocalSize(snes->vec_sol,&n);CHKERRQ(ierr); 725 ierr = PetscMalloc1(n,&indices);CHKERRQ(ierr); 726 for (i=0; i < n; i++) indices[i] = rstart + i; 727 ierr = ISCreateGeneral(PetscObjectComm((PetscObject)snes),n,indices,PETSC_OWN_POINTER,&vi->IS_inact_prev);CHKERRQ(ierr); 728 729 /* set the line search functions */ 730 if (!snes->linesearch) { 731 ierr = SNESGetLineSearch(snes, &linesearch);CHKERRQ(ierr); 732 ierr = SNESLineSearchSetType(linesearch, SNESLINESEARCHBT);CHKERRQ(ierr); 733 } 734 PetscFunctionReturn(0); 735 } 736 /* -------------------------------------------------------------------------- */ 737 PetscErrorCode SNESReset_VINEWTONRSLS(SNES snes) 738 { 739 SNES_VINEWTONRSLS *vi = (SNES_VINEWTONRSLS*) snes->data; 740 PetscErrorCode ierr; 741 742 PetscFunctionBegin; 743 ierr = SNESReset_VI(snes);CHKERRQ(ierr); 744 ierr = ISDestroy(&vi->IS_inact_prev);CHKERRQ(ierr); 745 PetscFunctionReturn(0); 746 } 747 748 /* -------------------------------------------------------------------------- */ 749 /*MC 750 SNESVINEWTONRSLS - Reduced space active set solvers for variational inequalities based on Newton's method 751 752 Options Database: 753 + -snes_type <vinewtonssls,vinewtonrsls> - a semi-smooth solver, a reduced space active set method 754 - -snes_vi_monitor - prints the number of active constraints at each iteration. 755 756 Level: beginner 757 758 References: 759 . 1. - T. S. Munson, and S. Benson. Flexible Complementarity Solvers for Large Scale 760 Applications, Optimization Methods and Software, 21 (2006). 761 762 .seealso: SNESVISetVariableBounds(), SNESVISetComputeVariableBounds(), SNESCreate(), SNES, SNESSetType(), SNESVINEWTONSSLS, SNESNEWTONTR, SNESLineSearchSetType(),SNESLineSearchSetPostCheck(), SNESLineSearchSetPreCheck() 763 764 M*/ 765 PETSC_EXTERN PetscErrorCode SNESCreate_VINEWTONRSLS(SNES snes) 766 { 767 PetscErrorCode ierr; 768 SNES_VINEWTONRSLS *vi; 769 770 PetscFunctionBegin; 771 snes->ops->reset = SNESReset_VINEWTONRSLS; 772 snes->ops->setup = SNESSetUp_VINEWTONRSLS; 773 snes->ops->solve = SNESSolve_VINEWTONRSLS; 774 snes->ops->destroy = SNESDestroy_VI; 775 snes->ops->setfromoptions = SNESSetFromOptions_VI; 776 snes->ops->view = NULL; 777 snes->ops->converged = SNESConvergedDefault_VI; 778 779 snes->usesksp = PETSC_TRUE; 780 snes->usesnpc = PETSC_FALSE; 781 782 snes->alwayscomputesfinalresidual = PETSC_TRUE; 783 784 ierr = PetscNewLog(snes,&vi);CHKERRQ(ierr); 785 snes->data = (void*)vi; 786 vi->checkredundancy = NULL; 787 788 ierr = PetscObjectComposeFunction((PetscObject)snes,"SNESVISetVariableBounds_C",SNESVISetVariableBounds_VI);CHKERRQ(ierr); 789 ierr = PetscObjectComposeFunction((PetscObject)snes,"SNESVISetComputeVariableBounds_C",SNESVISetComputeVariableBounds_VI);CHKERRQ(ierr); 790 PetscFunctionReturn(0); 791 } 792 793