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