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