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