Actual source code: ex5s.c

  2: static char help[] = "2d Bratu problem in shared memory parallel with SNES.\n\
  3: We solve the  Bratu (SFI - solid fuel ignition) problem in a 2D rectangular\n\
  4: domain, uses SHARED MEMORY to evaluate the user function.\n\
  5: The command line options include:\n\
  6:   -par <parameter>, where <parameter> indicates the problem's nonlinearity\n\
  7:      problem SFI:  <parameter> = Bratu parameter (0 <= par <= 6.81)\n\
  8:   -mx <xg>, where <xg> = number of grid points in the x-direction\n\
  9:   -my <yg>, where <yg> = number of grid points in the y-direction\n\
 10:   -use_fortran_function: use Fortran coded function, rather than C\n";

 12: /*
 13:              This code compiles ONLY on SGI systems
 14:             ========================================
 15: */
 16: /*T
 17:    Concepts: SNES^parallel Bratu example
 18:    Concepts: shared memory
 19:    Processors: n
 20: T*/

 22: /*

 24:      Programming model: Combination of 
 25:         1) MPI message passing for PETSc routines
 26:         2) automatic loop parallism (using shared memory) for user
 27:            provided function.

 29:        While the user function is being evaluated all MPI processes except process
 30:      0 blocks. Process zero spawns nt threads to evaluate the user function. Once 
 31:      the user function is complete, the worker threads are suspended and all the MPI processes
 32:      continue.

 34:      Other useful options:

 36:        -snes_mf : use matrix free operator and no preconditioner
 37:        -snes_mf_operator : use matrix free operator but compute Jacobian via
 38:                            finite differences to form preconditioner

 40:        Environmental variable:

 42:          setenv MPC_NUM_THREADS nt <- set number of threads processor 0 should 
 43:                                       use to evaluate user provided function

 45:        Note: The number of MPI processes (set with the mpiexec option -np) can 
 46:        be set completely independently from the number of threads process 0 
 47:        uses to evaluate the function (though usually one would make them the same).
 48: */
 49: 
 50: /* ------------------------------------------------------------------------

 52:     Solid Fuel Ignition (SFI) problem.  This problem is modeled by
 53:     the partial differential equation
 54:   
 55:             -Laplacian u - lambda*exp(u) = 0,  0 < x,y < 1,
 56:   
 57:     with boundary conditions
 58:    
 59:              u = 0  for  x = 0, x = 1, y = 0, y = 1.
 60:   
 61:     A finite difference approximation with the usual 5-point stencil
 62:     is used to discretize the boundary value problem to obtain a nonlinear 
 63:     system of equations.

 65:     The uniprocessor version of this code is snes/examples/tutorials/ex4.c
 66:     A parallel distributed memory version is snes/examples/tutorials/ex5.c and ex5f.F

 68:   ------------------------------------------------------------------------- */

 70: /* 
 71:    Include "petscsnes.h" so that we can use SNES solvers.  Note that this
 72:    file automatically includes:
 73:      petscsys.h       - base PETSc routines   petscvec.h - vectors
 74:      petscmat.h - matrices
 75:      petscis.h     - index sets            petscksp.h - Krylov subspace methods
 76:      petscviewer.h - viewers               petscpc.h  - preconditioners
 77:      petscksp.h   - linear solvers
 78: */
 79: #include <petscsnes.h>

 81: /* 
 82:    User-defined application context - contains data needed by the 
 83:    application-provided call-back routines   FormFunction().
 84: */
 85: typedef struct {
 86:    PetscReal   param;          /* test problem parameter */
 87:    int         mx,my;          /* discretization in x, y directions */
 88:    int         rank;           /* processor rank */
 89: } AppCtx;

 91: /* 
 92:    User-defined routines
 93: */

 99: /* 
100:     The main program is written in C while the user provided function
101:  is given in both Fortran and C. The main program could also be written 
102:  in Fortran; the ONE PROBLEM is that VecGetArray() cannot be called from 
103:  Fortran on the SGI machines; thus the routine FormFunctionFortran() must
104:  be written in C.
105: */
106: int main(int argc,char **argv)
107: {
108:   SNES           snes;                /* nonlinear solver */
109:   Vec            x,r;                 /* solution, residual vectors */
110:   AppCtx         user;                /* user-defined work context */
111:   int            its;                 /* iterations for convergence */
112:   int            N,ierr,rstart,rend,*colors,i,ii,ri,rj;
113:   PetscErrorCode (*fnc)(SNES,Vec,Vec,void*);
114:   PetscReal      bratu_lambda_max = 6.81,bratu_lambda_min = 0.;
115:   MatFDColoring  fdcoloring;
116:   ISColoring     iscoloring;
117:   Mat            J;
118:   PetscScalar    zero = 0.0;
119:   PetscBool      flg;

121:   PetscInitialize(&argc,&argv,(char *)0,help);
122:   MPI_Comm_rank(PETSC_COMM_WORLD,&user.rank);

124:   /*
125:      Initialize problem parameters
126:   */
127:   user.mx = 4; user.my = 4; user.param = 6.0;
128:   PetscOptionsGetInt(PETSC_NULL,"-mx",&user.mx,PETSC_NULL);
129:   PetscOptionsGetInt(PETSC_NULL,"-my",&user.my,PETSC_NULL);
130:   PetscOptionsGetReal(PETSC_NULL,"-par",&user.param,PETSC_NULL);
131:   if (user.param >= bratu_lambda_max || user.param <= bratu_lambda_min) SETERRQ(PETSC_COMM_SELF,1,"Lambda is out of range");
132:   N = user.mx*user.my;

134:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
135:      Create nonlinear solver context
136:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

138:   SNESCreate(PETSC_COMM_WORLD,&snes);

140:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
141:      Create vector data structures; set function evaluation routine
142:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

144:   /*
145:       The routine VecCreateShared() creates a parallel vector with each processor
146:     assigned its own segment, BUT, in addition, the first processor has access to the 
147:     entire array. This is to allow the users function to be based on loop level
148:     parallelism rather than MPI.
149:   */
150:   VecCreateShared(PETSC_COMM_WORLD,PETSC_DECIDE,N,&x);
151:   VecDuplicate(x,&r);

153:   PetscOptionsHasName(PETSC_NULL,"-use_fortran_function",&flg);
154:   if (flg) {
155:     fnc = FormFunctionFortran;
156:   } else {
157:     fnc = FormFunction;
158:   }

160:   /* 
161:      Set function evaluation routine and vector
162:   */
163:   SNESSetFunction(snes,r,fnc,&user);

165:   /*
166:        Currently when using VecCreateShared() and using loop level parallelism
167:     to automatically parallelise the user function it makes no sense for the 
168:     Jacobian to be computed via loop level parallelism, because all the threads
169:     would be simultaneously calling MatSetValues() causing a bottle-neck.

171:     Thus this example uses the PETSc Jacobian calculations via finite differencing
172:     to approximate the Jacobian
173:   */

175:   /*

177:   */
178:   VecGetOwnershipRange(r,&rstart,&rend);
179:   PetscMalloc((rend-rstart)*sizeof(PetscInt),&colors);
180:   for (i=rstart; i<rend; i++) {
181:     colors[i - rstart] = 3*((i/user.mx) % 3) + (i % 3);
182:   }
183:   ISColoringCreate(PETSC_COMM_WORLD,3*2+2,rend-rstart,colors,&iscoloring);
184:   PetscFree(colors);

186:   /*
187:      Create and set the nonzero pattern for the Jacobian: This is not done 
188:      particularly efficiently. One should process the boundary nodes separately and 
189:      then use a simple loop for the interior nodes.
190:        Note that for this code we use the "natural" number of the nodes on the 
191:      grid (since that is what is good for the user provided function). In the 
192:      DMDA examples we must use the DMDA numbering where each processor is assigned a
193:      chunk of data.
194:   */
195:   MatCreateMPIAIJ(PETSC_COMM_WORLD,rend-rstart,rend-rstart,N,
196:                          N,5,0,0,0,&J);
197:   for (i=rstart; i<rend; i++) {
198:     rj = i % user.mx;         /* column in grid */
199:     ri = i / user.mx;         /* row in grid */
200:     if (ri != 0) {     /* first row does not have neighbor below */
201:       ii   = i - user.mx;
202:       MatSetValues(J,1,&i,1,&ii,&zero,INSERT_VALUES);
203:     }
204:     if (ri != user.my - 1) { /* last row does not have neighbors above */
205:       ii   = i + user.mx;
206:       MatSetValues(J,1,&i,1,&ii,&zero,INSERT_VALUES);
207:     }
208:     if (rj != 0) {     /* first column does not have neighbor to left */
209:       ii   = i - 1;
210:       MatSetValues(J,1,&i,1,&ii,&zero,INSERT_VALUES);
211:     }
212:     if (rj != user.mx - 1) {     /* last column does not have neighbor to right */
213:       ii   = i + 1;
214:       MatSetValues(J,1,&i,1,&ii,&zero,INSERT_VALUES);
215:     }
216:     MatSetValues(J,1,&i,1,&i,&zero,INSERT_VALUES);
217:   }
218:   MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);
219:   MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);

221:   /*
222:        Create the data structure that SNESDefaultComputeJacobianColor() uses
223:        to compute the actual Jacobians via finite differences.
224:   */
225:   MatFDColoringCreate(J,iscoloring,&fdcoloring);
226:   MatFDColoringSetFunction(fdcoloring,(PetscErrorCode (*)(void))fnc,&user);
227:   MatFDColoringSetFromOptions(fdcoloring);
228:   /*
229:         Tell SNES to use the routine SNESDefaultComputeJacobianColor()
230:       to compute Jacobians.
231:   */
232:   SNESSetJacobian(snes,J,J,SNESDefaultComputeJacobianColor,fdcoloring);
233:   ISColoringDestroy(&iscoloring);


236:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
237:      Customize nonlinear solver; set runtime options
238:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

240:   /*
241:      Set runtime options (e.g., -snes_monitor -snes_rtol <rtol> -ksp_type <type>)
242:   */
243:   SNESSetFromOptions(snes);

245:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
246:      Evaluate initial guess; then solve nonlinear system
247:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
248:   /*
249:      Note: The user should initialize the vector, x, with the initial guess
250:      for the nonlinear solver prior to calling SNESSolve().  In particular,
251:      to employ an initial guess of zero, the user should explicitly set
252:      this vector to zero by calling VecSet().
253:   */
254:   FormInitialGuess(&user,x);
255:   SNESSolve(snes,PETSC_NULL,x);
256:   SNESGetIterationNumber(snes,&its);
257:   PetscPrintf(PETSC_COMM_WORLD,"Number of Newton iterations = %D\n",its);

259:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
260:      Free work space.  All PETSc objects should be destroyed when they
261:      are no longer needed.
262:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
263:   VecDestroy(&x);
264:   VecDestroy(&r);
265:   SNESDestroy(&snes);
266:   PetscFinalize();

268:   return 0;
269: }
270: /* ------------------------------------------------------------------- */

274: /* 
275:    FormInitialGuess - Forms initial approximation.

277:    Input Parameters:
278:    user - user-defined application context
279:    X - vector

281:    Output Parameter:
282:    X - vector
283:  */
284: int FormInitialGuess(AppCtx *user,Vec X)
285: {
286:   int          i,j,row,mx,my,ierr;
287:   PetscReal    one = 1.0,lambda,temp1,temp,hx,hy,hxdhy,hydhx,sc;
288:   PetscScalar  *x;

290:   /*
291:       Process 0 has to wait for all other processes to get here 
292:    before proceeding to write in the shared vector
293:   */
294:   PetscBarrier((PetscObject)X);
295:   if (user->rank) {
296:      /*
297:         All the non-busy processors have to wait here for process 0 to finish
298:         evaluating the function; otherwise they will start using the vector values
299:         before they have been computed
300:      */
301:      PetscBarrier((PetscObject)X);
302:      return 0;
303:   }

305:   mx = user->mx;            my = user->my;            lambda = user->param;
306:   hx = one/(PetscReal)(mx-1);  hy = one/(PetscReal)(my-1);
307:   sc = hx*hy*lambda;        hxdhy = hx/hy;            hydhx = hy/hx;
308:   temp1 = lambda/(lambda + one);

310:   /*
311:      Get a pointer to vector data.
312:        - For default PETSc vectors, VecGetArray() returns a pointer to
313:          the data array.  Otherwise, the routine is implementation dependent.
314:        - You MUST call VecRestoreArray() when you no longer need access to
315:          the array.
316:   */
317:   VecGetArray(X,&x);

319:   /*
320:      Compute initial guess over the locally owned part of the grid
321:   */
322: #pragma arl(4)
323: #pragma distinct (*x,*f)
324: #pragma no side effects (sqrt)
325:   for (j=0; j<my; j++) {
326:     temp = (PetscReal)(PetscMin(j,my-j-1))*hy;
327:     for (i=0; i<mx; i++) {
328:       row = i + j*mx;
329:       if (i == 0 || j == 0 || i == mx-1 || j == my-1) {
330:         x[row] = 0.0;
331:         continue;
332:       }
333:       x[row] = temp1*sqrt(PetscMin((PetscReal)(PetscMin(i,mx-i-1))*hx,temp));
334:     }
335:   }

337:   /*
338:      Restore vector
339:   */
340:   VecRestoreArray(X,&x);

342:   PetscBarrier((PetscObject)X);
343:   return 0;
344: }
345: /* ------------------------------------------------------------------- */
348: /* 
349:    FormFunction - Evaluates nonlinear function, F(x).

351:    Input Parameters:
352: .  snes - the SNES context
353: .  X - input vector
354: .  ptr - optional user-defined context, as set by SNESSetFunction()

356:    Output Parameter:
357: .  F - function vector
358:  */
359: int FormFunction(SNES snes,Vec X,Vec F,void *ptr)
360: {
361:   AppCtx       *user = (AppCtx*)ptr;
362:   int          ierr,i,j,row,mx,my;
363:   PetscReal    two = 2.0,one = 1.0,lambda,hx,hy,hxdhy,hydhx,sc;
364:   PetscScalar  u,uxx,uyy,*x,*f;

366:   /*
367:       Process 0 has to wait for all other processes to get here 
368:    before proceeding to write in the shared vector
369:   */
370:   PetscBarrier((PetscObject)X);

372:   if (user->rank) {
373:      /*
374:         All the non-busy processors have to wait here for process 0 to finish
375:         evaluating the function; otherwise they will start using the vector values
376:         before they have been computed
377:      */
378:      PetscBarrier((PetscObject)X);
379:      return 0;
380:   }

382:   mx = user->mx;            my = user->my;            lambda = user->param;
383:   hx = one/(PetscReal)(mx-1);  hy = one/(PetscReal)(my-1);
384:   sc = hx*hy*lambda;        hxdhy = hx/hy;            hydhx = hy/hx;

386:   /*
387:      Get pointers to vector data
388:   */
389:   VecGetArray(X,&x);
390:   VecGetArray(F,&f);

392:   /*
393:       The next line tells the SGI compiler that x and f contain no overlapping 
394:     regions and thus it can use addition optimizations.
395:   */
396: #pragma arl(4)
397: #pragma distinct (*x,*f)
398: #pragma no side effects (exp)

400:   /*
401:      Compute function over the entire  grid
402:   */
403:   for (j=0; j<my; j++) {
404:     for (i=0; i<mx; i++) {
405:       row = i + j*mx;
406:       if (i == 0 || j == 0 || i == mx-1 || j == my-1) {
407:         f[row] = x[row];
408:         continue;
409:       }
410:       u = x[row];
411:       uxx = (two*u - x[row-1] - x[row+1])*hydhx;
412:       uyy = (two*u - x[row-mx] - x[row+mx])*hxdhy;
413:       f[row] = uxx + uyy - sc*exp(u);
414:     }
415:   }

417:   /*
418:      Restore vectors
419:   */
420:   VecRestoreArray(X,&x);
421:   VecRestoreArray(F,&f);

423:   PetscLogFlops(11.0*(mx-2)*(my-2))
424:   PetscBarrier((PetscObject)X);
425:   return 0;
426: }

428: #if defined(PETSC_HAVE_FORTRAN_CAPS)
429: #define applicationfunctionfortran_ APPLICATIONFUNCTIONFORTRAN
430: #elif !defined(PETSC_HAVE_FORTRAN_UNDERSCORE)
431: #define applicationfunctionfortran_ applicationfunctionfortran
432: #endif

434: /* ------------------------------------------------------------------- */
437: /* 
438:    FormFunctionFortran - Evaluates nonlinear function, F(x) in Fortran.

440: */
441: int FormFunctionFortran(SNES snes,Vec X,Vec F,void *ptr)
442: {
443:   AppCtx  *user = (AppCtx*)ptr;
444:   int     ierr;
445:   PetscScalar  *x,*f;

447:   /*
448:       Process 0 has to wait for all other processes to get here 
449:    before proceeding to write in the shared vector
450:   */
451:   PetscBarrier((PetscObject)snes);
452:   if (!user->rank) {
453:     VecGetArray(X,&x);
454:     VecGetArray(F,&f);
455:     applicationfunctionfortran_(&user->param,&user->mx,&user->my,x,f,&ierr);
456:     VecRestoreArray(X,&x);
457:     VecRestoreArray(F,&f);
458:     PetscLogFlops(11.0*(user->mx-2)*(user->my-2))
459:   }
460:   /*
461:       All the non-busy processors have to wait here for process 0 to finish
462:       evaluating the function; otherwise they will start using the vector values
463:       before they have been computed
464:   */
465:   PetscBarrier((PetscObject)snes);
466:   return 0;
467: }