Actual source code: ex46.c

petsc-master 2016-09-28
Report Typos and Errors
  1: static char help[] = "Surface processes in geophysics.\n\n";

  3: /*T
  4:    Concepts: SNES^parallel Surface process example
  5:    Concepts: DMDA^using distributed arrays;
  6:    Concepts: IS coloirng types;
  7:    Processors: n
  8: T*/


 11:  #include <petscsnes.h>
 12:  #include <petscdm.h>
 13:  #include <petscdmda.h>

 15: /*
 16:    User-defined application context - contains data needed by the
 17:    application-provided call-back routines, FormJacobianLocal() and
 18:    FormFunctionLocal().
 19: */
 20: typedef struct {
 21:   PetscReal   D;  /* The diffusion coefficient */
 22:   PetscReal   K;  /* The advection coefficient */
 23:   PetscInt    m;  /* Exponent for A */
 24: } AppCtx;

 26: /*
 27:    User-defined routines
 28: */
 29: extern PetscErrorCode FormFunctionLocal(DMDALocalInfo*,PetscScalar**,PetscScalar**,AppCtx*);
 30: extern PetscErrorCode FormJacobianLocal(DMDALocalInfo*,PetscScalar**,Mat,AppCtx*);

 34: int main(int argc,char **argv)
 35: {
 36:   SNES           snes;                         /* nonlinear solver */
 37:   AppCtx         user;                         /* user-defined work context */
 38:   PetscInt       its;                          /* iterations for convergence */
 40:   DM             da;

 42:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 43:      Initialize program
 44:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

 46:   PetscInitialize(&argc,&argv,(char*)0,help);if (ierr) return ierr;

 48:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 49:      Initialize problem parameters
 50:   - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
 51:   PetscOptionsBegin(PETSC_COMM_WORLD, "", "Surface Process Problem Options", "SNES");
 52:   user.D = 1.0;
 53:   PetscOptionsReal("-D", "The diffusion coefficient D", __FILE__, user.D, &user.D, NULL);
 54:   user.K = 1.0;
 55:   PetscOptionsReal("-K", "The advection coefficient K", __FILE__, user.K, &user.K, NULL);
 56:   user.m = 1;
 57:   PetscOptionsInt("-m", "The exponent for A", __FILE__, user.m, &user.m, NULL);
 58:   PetscOptionsEnd();

 60:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 61:      Create distributed array (DMDA) to manage parallel grid and vectors
 62:   - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
 63:   DMDACreate2d(PETSC_COMM_WORLD, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE,DMDA_STENCIL_STAR,4,4,PETSC_DECIDE,PETSC_DECIDE,1,1,NULL,NULL,&da);
 64:   DMSetFromOptions(da);
 65:   DMSetUp(da);
 66:   DMDASetUniformCoordinates(da, 0.0, 1.0, 0.0, 1.0, 0.0, 1.0);
 67:   DMSetApplicationContext(da,&user);
 68:   SNESCreate(PETSC_COMM_WORLD, &snes);
 69:   SNESSetDM(snes, da);

 71:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 72:      Set local function evaluation routine
 73:   - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
 74:   DMDASNESSetFunctionLocal(da,INSERT_VALUES,(PetscErrorCode (*)(DMDALocalInfo*,void*,void*,void*))FormFunctionLocal,&user);

 76:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 77:      Customize solver; set runtime options
 78:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
 79:   SNESSetFromOptions(snes);


 82:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 83:      Solve nonlinear system
 84:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
 85:   SNESSolve(snes,0,0);
 86:   SNESGetIterationNumber(snes,&its);

 88:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 89:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
 90:   PetscPrintf(PETSC_COMM_WORLD,"Number of SNES iterations = %D\n",its);

 92:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 93:      Free work space.  All PETSc objects should be destroyed when they
 94:      are no longer needed.
 95:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

 97:   SNESDestroy(&snes);
 98:   DMDestroy(&da);

100:   PetscFinalize();
101:   return ierr;
102: }

106: PetscScalar funcU(DMDACoor2d *coords)
107: {
108:   return coords->x + coords->y;
109: }

113: PetscScalar funcA(PetscScalar z, AppCtx *user)
114: {
115:   PetscScalar v = 1.0;
116:   PetscInt    i;

118:   for (i = 0; i < user->m; ++i) v *= z;
119:   return v;
120: }

124: PetscScalar funcADer(PetscScalar z, AppCtx *user)
125: {
126:   PetscScalar v = 1.0;
127:   PetscInt    i;

129:   for (i = 0; i < user->m-1; ++i) v *= z;
130:   return (PetscScalar)user->m*v;
131: }

135: /*
136:    FormFunctionLocal - Evaluates nonlinear function, F(x).
137: */
138: PetscErrorCode FormFunctionLocal(DMDALocalInfo *info,PetscScalar **x,PetscScalar **f,AppCtx *user)
139: {
140:   DM             coordDA;
141:   Vec            coordinates;
142:   DMDACoor2d     **coords;
143:   PetscScalar    u, ux, uy, uxx, uyy;
144:   PetscReal      D, K, hx, hy, hxdhy, hydhx;
145:   PetscInt       i,j;

149:   D     = user->D;
150:   K     = user->K;
151:   hx    = 1.0/(PetscReal)(info->mx-1);
152:   hy    = 1.0/(PetscReal)(info->my-1);
153:   hxdhy = hx/hy;
154:   hydhx = hy/hx;
155:   /*
156:      Compute function over the locally owned part of the grid
157:   */
158:   DMGetCoordinateDM(info->da, &coordDA);
159:   DMGetCoordinates(info->da, &coordinates);
160:   DMDAVecGetArray(coordDA, coordinates, &coords);
161:   for (j=info->ys; j<info->ys+info->ym; j++) {
162:     for (i=info->xs; i<info->xs+info->xm; i++) {
163:       if (i == 0 || j == 0 || i == info->mx-1 || j == info->my-1) f[j][i] = x[j][i];
164:       else {
165:         u       = x[j][i];
166:         ux      = (x[j][i+1] - x[j][i])/hx;
167:         uy      = (x[j+1][i] - x[j][i])/hy;
168:         uxx     = (2.0*u - x[j][i-1] - x[j][i+1])*hydhx;
169:         uyy     = (2.0*u - x[j-1][i] - x[j+1][i])*hxdhy;
170:         f[j][i] = D*(uxx + uyy) - (K*funcA(x[j][i], user)*PetscSqrtScalar(ux*ux + uy*uy) + funcU(&coords[j][i]))*hx*hy;
171:         if (PetscIsInfOrNanScalar(f[j][i])) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_FP, "Invalid residual: %g", PetscRealPart(f[j][i]));
172:       }
173:     }
174:   }
175:   DMDAVecRestoreArray(coordDA, coordinates, &coords);
176:   PetscLogFlops(11*info->ym*info->xm);
177:   return(0);
178: }

182: /*
183:    FormJacobianLocal - Evaluates Jacobian matrix.
184: */
185: PetscErrorCode FormJacobianLocal(DMDALocalInfo *info,PetscScalar **x,Mat jac,AppCtx *user)
186: {
187:   MatStencil     col[5], row;
188:   PetscScalar    D, K, A, v[5], hx, hy, hxdhy, hydhx, ux, uy;
189:   PetscReal      normGradZ;
190:   PetscInt       i, j,k;

194:   D     = user->D;
195:   K     = user->K;
196:   hx    = 1.0/(PetscReal)(info->mx-1);
197:   hy    = 1.0/(PetscReal)(info->my-1);
198:   hxdhy = hx/hy;
199:   hydhx = hy/hx;

201:   /*
202:      Compute entries for the locally owned part of the Jacobian.
203:       - Currently, all PETSc parallel matrix formats are partitioned by
204:         contiguous chunks of rows across the processors.
205:       - Each processor needs to insert only elements that it owns
206:         locally (but any non-local elements will be sent to the
207:         appropriate processor during matrix assembly).
208:       - Here, we set all entries for a particular row at once.
209:       - We can set matrix entries either using either
210:         MatSetValuesLocal() or MatSetValues(), as discussed above.
211:   */
212:   for (j=info->ys; j<info->ys+info->ym; j++) {
213:     for (i=info->xs; i<info->xs+info->xm; i++) {
214:       row.j = j; row.i = i;
215:       if (i == 0 || j == 0 || i == info->mx-1 || j == info->my-1) {
216:         /* boundary points */
217:         v[0] = 1.0;
218:         MatSetValuesStencil(jac,1,&row,1,&row,v,INSERT_VALUES);
219:       } else {
220:         /* interior grid points */
221:         ux        = (x[j][i+1] - x[j][i])/hx;
222:         uy        = (x[j+1][i] - x[j][i])/hy;
223:         normGradZ = PetscRealPart(PetscSqrtScalar(ux*ux + uy*uy));
224:         /* PetscPrintf(PETSC_COMM_SELF, "i: %d j: %d normGradZ: %g\n", i, j, normGradZ); */
225:         if (normGradZ < 1.0e-8) normGradZ = 1.0e-8;
226:         A = funcA(x[j][i], user);

228:         v[0] = -D*hxdhy;                                                                          col[0].j = j - 1; col[0].i = i;
229:         v[1] = -D*hydhx;                                                                          col[1].j = j;     col[1].i = i-1;
230:         v[2] = D*2.0*(hydhx + hxdhy) + K*(funcADer(x[j][i], user)*normGradZ - A/normGradZ)*hx*hy; col[2].j = row.j; col[2].i = row.i;
231:         v[3] = -D*hydhx + K*A*hx*hy/(2.0*normGradZ);                                              col[3].j = j;     col[3].i = i+1;
232:         v[4] = -D*hxdhy + K*A*hx*hy/(2.0*normGradZ);                                              col[4].j = j + 1; col[4].i = i;
233:         for (k = 0; k < 5; ++k) {
234:           if (PetscIsInfOrNanScalar(v[k])) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_FP, "Invalid residual: %g", PetscRealPart(v[k]));
235:         }
236:         MatSetValuesStencil(jac,1,&row,5,col,v,INSERT_VALUES);
237:       }
238:     }
239:   }

241:   /*
242:      Assemble matrix, using the 2-step process:
243:        MatAssemblyBegin(), MatAssemblyEnd().
244:   */
245:   MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY);
246:   MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY);
247:   /*
248:      Tell the matrix we will never add a new nonzero location to the
249:      matrix. If we do, it will generate an error.
250:   */
251:   MatSetOption(jac,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE);
252:   return(0);
253: }