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