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