Actual source code: ex43-44.h

  1: /*
  2:   Accelerates Newton's method by solving a small problem defined by those elements with large residual plus one level of overlap

  4:   This is a toy code for playing around

  6:   Counts residual entries as small if they are less then .2 times the maximum
  7:   Decides to solve a reduced problem if the number of large entries is less than 20 percent of all entries (and this has been true for criteria_reduce iterations)
  8: */
 9:  #include petscsnes.h



 15: typedef struct {
 16:   Vec        xwork,fwork;
 17:   VecScatter scatter;
 18:   SNES       snes;
 19:   IS         is;
 20: } SubCtx;

 24: PetscErrorCode FormFunctionSub(SNES snes,Vec x,Vec f,void *ictx)
 25: {
 27:   SubCtx         *ctx = (SubCtx*) ictx;

 30:   VecScatterBegin(ctx->scatter,x,ctx->xwork,INSERT_VALUES,SCATTER_REVERSE);
 31:   VecScatterEnd(ctx->scatter,x,ctx->xwork,INSERT_VALUES,SCATTER_REVERSE);
 32:   SNESComputeFunction(ctx->snes,ctx->xwork,ctx->fwork);
 33:   VecScatterBegin(ctx->scatter,ctx->fwork,f,INSERT_VALUES,SCATTER_FORWARD);
 34:   VecScatterEnd(ctx->scatter,ctx->fwork,f,INSERT_VALUES,SCATTER_FORWARD);
 35:   return(0);
 36: }

 40: PetscErrorCode FormJacobianSub(SNES snes,Vec x,Mat *A, Mat *B, MatStructure *str,void *ictx)
 41: {
 43:   SubCtx         *ctx = (SubCtx*) ictx;
 44:   Mat            As,Bs;

 47:   VecScatterBegin(ctx->scatter,x,ctx->xwork,INSERT_VALUES,SCATTER_REVERSE);
 48:   VecScatterEnd(ctx->scatter,x,ctx->xwork,INSERT_VALUES,SCATTER_REVERSE);
 49:   SNESGetJacobian(ctx->snes,&As,&Bs,PETSC_NULL,PETSC_NULL);
 50:   SNESComputeJacobian(ctx->snes,ctx->xwork,&As,&Bs,str);
 51:   if (*B) {
 52:     MatGetSubMatrix(Bs,ctx->is,ctx->is,MAT_REUSE_MATRIX,B);
 53:   } else {
 54:     MatGetSubMatrix(Bs,ctx->is,ctx->is,MAT_INITIAL_MATRIX,B);
 55:   }
 56:   if (!*A) {
 57:     *A = *B;
 58:     PetscObjectReference((PetscObject)*A);
 59:   }
 60:   if (*A != *B) {
 61:     MatAssemblyBegin(*A,MAT_FINAL_ASSEMBLY);
 62:     MatAssemblyEnd(*A,MAT_FINAL_ASSEMBLY);
 63:   }
 64:   return(0);
 65: }

 69: PetscErrorCode SolveSubproblem(SNES snes)
 70: {
 72:   Vec            residual,solution;
 73:   PetscReal      rmax;
 74:   PetscInt       i,n,cnt,*indices;
 75:   PetscScalar    *r;
 76:   SNES           snessub;
 77:   Vec            x,f;
 78:   SubCtx         ctx;
 79:   Mat            mat;

 82:   ctx.snes = snes;
 83:   SNESGetSolution(snes,&solution);
 84:   SNESGetFunction(snes,&residual,0,0);
 85:   VecNorm(residual,NORM_INFINITY,&rmax);
 86:   VecGetLocalSize(residual,&n);
 87:   VecGetArray(residual,&r);
 88:   cnt  = 0;
 89:   for (i=0; i<n; i++) {
 90:     if (PetscAbsScalar(r[i]) > .20*rmax ) cnt++;
 91:   }
 92:   PetscMalloc(cnt*sizeof(PetscInt),&indices);
 93:   cnt  = 0;
 94:   for (i=0; i<n; i++) {
 95:     if (PetscAbsScalar(r[i]) > .20*rmax ) indices[cnt++] = i;
 96:   }
 97:   if (cnt > .2*n) return(0);

 99:   printf("number in subproblem %d\n",cnt);
100:   VecRestoreArray(residual,&r);
101:   ISCreateGeneral(PETSC_COMM_WORLD,cnt,indices,&ctx.is);
102:   PetscFree(indices);

104:   SNESGetJacobian(snes,0,&mat,0,0);
105:   MatIncreaseOverlap(mat,1,&ctx.is,2);
106:   ISSort(ctx.is);
107:   ISGetLocalSize(ctx.is,&cnt);
108:   printf("number in subproblem %d\n",cnt);
109:   VecCreate(PETSC_COMM_WORLD,&x);
110:   VecSetSizes(x,cnt,cnt);
111:   VecSetType(x,VECSEQ);
112:   VecDuplicate(x,&f);
113:   VecScatterCreate(solution,ctx.is,x,PETSC_NULL,&ctx.scatter);

115:   VecDuplicate(solution,&ctx.xwork);
116:   VecCopy(solution,ctx.xwork);
117:   VecDuplicate(residual,&ctx.fwork);
118:   VecScatterBegin(ctx.scatter,solution,x,INSERT_VALUES,SCATTER_FORWARD);
119:   VecScatterEnd(ctx.scatter,solution,x,INSERT_VALUES,SCATTER_FORWARD);
120:   SNESCreate(PETSC_COMM_WORLD,&snessub);
121:   SNESSetOptionsPrefix(snessub,"sub_");
122:   SNESSetFunction(snessub,f,FormFunctionSub,(void*)&ctx);
123:   SNESSetJacobian(snessub,0,0,FormJacobianSub,(void*)&ctx);
124:   SNESSetFromOptions(snessub);
125:   SNESSolve(snessub,PETSC_NULL,x);
126:   VecScatterBegin(ctx.scatter,x,solution,INSERT_VALUES,SCATTER_REVERSE);
127:   VecScatterEnd(ctx.scatter,x,solution,INSERT_VALUES,SCATTER_REVERSE);

129:   ISDestroy(ctx.is);
130:   VecDestroy(ctx.xwork);
131:   VecDestroy(ctx.fwork);
132:   VecScatterDestroy(ctx.scatter);
133:   VecDestroy(x);
134:   VecDestroy(f);
135:   SNESDestroy(snessub);
136:   return(0);
137: }


141: static PetscInt CountGood = 0;

145: PetscErrorCode  MonitorRange(SNES snes,PetscInt it,PetscReal rnorm,void *dummy)
146: {
147:   PetscErrorCode          ierr;
148:   PetscReal               perc;

150:   SNESMonitorRange_Private(snes,it,&perc);
151:   if (perc < .20) CountGood++;
152:   else CountGood = 0;
153:   return(0);
154: }