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