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