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