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