Actual source code: ex26.c

  1: /*$Id: gs-cylinder.c, xtang Exp $*/

  3: static char help[] = "Grad-Shafranov solver for one dimensional CHI equilibrium.n
  4: The command line options include:n
  5:   -n <n> number of grid pointsn
  6:   -psi_axis <axis> n
  7:   -r_min <min> n
  8:   -param <param> nn";

 10: /*T
 11:    Concepts: SNES^parallel CHI equilibrium
 12:    Concepts: DA^using distributed arrays;
 13:    Processors: n
 14: T*/

 16: /* ------------------------------------------------------------------------

 18:    Grad-Shafranov solver for one dimensional CHI equilibrium
 19:   
 20:     A finite difference approximation with the usual 3-point stencil
 21:     is used to discretize the boundary value problem to obtain a nonlinear 
 22:     system of equations.

 24:     Contributed by Xianzhu Tang, LANL

 26:     An interesting feature of this example is that as you refine the grid
 27:     (with a larger -n <n> you cannot force the residual norm as small. This
 28:     appears to be due to "NOISE" in the function, the FormFunctionLocal() cannot
 29:     be computed as accurately with a finer grid.
 30:   ------------------------------------------------------------------------- */

 32:  #include petscda.h
 33:  #include petscsnes.h

 35: /* 
 36:    User-defined application context - contains data needed by the 
 37:    application-provided call-back routines, FormJacobian() and
 38:    FormFunction().
 39: */
 40: typedef struct {
 41:   DA            da;               /* distributed array data structure */
 42:   Vec           psi,r;            /* solution, residual vectors */
 43:   Mat           A,J;              /* Jacobian matrix */
 44:   Vec           coordinates;      /* grid coordinates */
 45:   PassiveReal   psi_axis,psi_bdy;
 46:   PassiveReal   r_min;
 47:   PassiveReal   param;            /* test problem parameter */
 48: } AppCtx;

 50: #define GdGdPsi(r,u)      (((r) < 0.05) ? 0.0 : (user->param*((r)-0.05)*(1.0-(u)*(u))*(1.0-(u)*(u))))
 51: #define CurrentWire(r)    (((r) < .05) ? -3.E2 : 0.0)
 52: #define GdGdPsiPrime(r,u) (((r) < 0.05) ? 0.0 : -4.*(user->param*((r)-0.05)*(1.0-(u)*(u)))*u)

 54: /* 
 55:    User-defined routines
 56: */
 57: extern int FormInitialGuess(AppCtx*,Vec);
 58: extern int FormJacobian(SNES,Vec,Mat*,Mat*,MatStructure*,void*);
 59: extern int FormFunctionLocal(DALocalInfo*,PetscScalar*,PetscScalar*,AppCtx*);

 61: #undef __FUNCT__
 63: int main(int argc,char **argv)
 64: {
 65:   SNES                   snes;                 /* nonlinear solver */
 66:   AppCtx                 user;                 /* user-defined work context */
 67:   int                    its;                  /* iterations for convergence */
 68:   PetscTruth             fd_jacobian = PETSC_FALSE;
 69:   PetscTruth             adicmf_jacobian = PETSC_FALSE;
 70:   int                    grids = 100, dof = 1, stencil_width = 1;
 71:   int                    ierr;
 72:   PetscReal              fnorm;
 73:   MatFDColoring          matfdcoloring = 0;
 74:   ISColoring             iscoloring;

 76:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 77:      Initialize program
 78:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

 80:   PetscInitialize(&argc,&argv,(char *)0,help);
 81: 
 82:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 83:      Initialize problem parameters
 84:   - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

 86:   user.psi_axis=0.0;
 87:   PetscOptionsGetReal(PETSC_NULL,"-psi_axis",&user.psi_axis,PETSC_NULL);
 88:   user.psi_bdy=1.0;
 89:   PetscOptionsGetReal(PETSC_NULL,"-psi_bdy",&user.psi_bdy,PETSC_NULL);
 90:   user.r_min=0.0;
 91:   PetscOptionsGetReal(PETSC_NULL,"-r_min",&user.r_min,PETSC_NULL);
 92:   user.param=-1.E1;
 93:   PetscOptionsGetReal(PETSC_NULL,"-param",&user.param,PETSC_NULL);

 95:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 96:      Create nonlinear solver context
 97:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
 98:   SNESCreate(PETSC_COMM_WORLD,&snes);
 99: 
100:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
101:      Create distributed array (DA) to manage parallel grid and vectors
102:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
103:   PetscOptionsGetInt(PETSC_NULL,"-n",&grids,PETSC_NULL);
104:   DACreate1d(PETSC_COMM_WORLD,DA_NONPERIODIC,grids,dof,stencil_width,PETSC_NULL,&user.da);
105: 
106:   /*  - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
107:      Extract global vectors from DA; then duplicate for remaining
108:      vectors that are the same types
109:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
110:   DACreateGlobalVector(user.da,&user.psi);
111:   VecDuplicate(user.psi,&user.r);

113:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
114:      Set local function evaluation routine
115:   - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

117:   DASetLocalFunction(user.da,(DALocalFunction1)FormFunctionLocal);

119:   SNESSetFunction(snes,user.r,SNESDAFormFunction,&user);

121:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
122:      Create matrix data structure; set Jacobian evaluation routine

124:      Set Jacobian matrix data structure and default Jacobian evaluation
125:      routine. User can override with:
126:      -snes_mf : matrix-free Newton-Krylov method with no preconditioning
127:                 (unless user explicitly sets preconditioner) 
128:      -snes_mf_operator : form preconditioning matrix as set by the user,
129:                          but use matrix-free approx for Jacobian-vector
130:                          products within Newton-Krylov method

132:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
133:   PetscOptionsGetLogical(PETSC_NULL,"-fd_jacobian",&fd_jacobian,0);
134:   /*
135:        Note that fd_jacobian DOES NOT compute the finite difference approximation to 
136:     the ENTIRE Jacobian. Rather it removes the global coupling from the Jacobian and
137:     computes the finite difference approximation only for the "local" coupling.

139:        Thus running with fd_jacobian and not -snes_mf_operator or -adicmf_jacobian
140:     won't converge.
141:   */
142:   if (!fd_jacobian) {
143:     ierr      = MatCreateMPIAIJ(PETSC_COMM_WORLD,PETSC_DECIDE,PETSC_DECIDE,grids,grids,
144:                                 5,PETSC_NULL,3,PETSC_NULL,&user.J);
145:     ierr      = MatSetFromOptions(user.J);
146:     user.A    = user.J;
147:   }
148:   else {
149:     ierr      = DAGetMatrix(user.da,MATMPIAIJ,&user.J);
150:     user.A    = user.J;
151:   }

153:   PetscOptionsGetLogical(PETSC_NULL,"-adicmf_jacobian",&adicmf_jacobian,0);
154: #if defined(PETSC_HAVE_ADIC) && !defined(PETSC_USE_COMPLEX) && !defined(PETSC_USE_SINGLE)
155:   if (adicmf_jacobian) {
156:     DASetLocalAdicMFFunction(user.da,admf_FormFunctionLocal);
157:     MatRegisterDAAD();
158:     MatCreateDAAD(user.da,&user.A);
159:     MatDAADSetSNES(user.A,snes);
160:     MatDAADSetCtx(user.A,&user);
161:   }
162: #endif

164:   if (fd_jacobian) {
165:     DAGetColoring(user.da,IS_COLORING_LOCAL,&iscoloring);
166:     MatFDColoringCreate(user.J,iscoloring,&matfdcoloring);
167:     ISColoringDestroy(iscoloring);
168:     MatFDColoringSetFunction(matfdcoloring,(int (*)(void))SNESDAFormFunction,&user);
169:     MatFDColoringSetFromOptions(matfdcoloring);
170:     SNESSetJacobian(snes,user.A,user.J,SNESDefaultComputeJacobianColor,matfdcoloring);
171:   } else {
172:     SNESSetJacobian(snes,user.A,user.J,FormJacobian,&user);
173:   }

175:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
176:      Customize nonlinear solver; set runtime options
177:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
178:   SNESSetFromOptions(snes);

180:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
181:      Evaluate initial guess
182:      Note: The user should initialize the vector, x, with the initial guess
183:      for the nonlinear solver prior to calling SNESSolve().  In particular,
184:      to employ an initial guess of zero, the user should explicitly set
185:      this vector to zero by calling VecSet().
186:   */
187:   FormInitialGuess(&user,user.psi);

189:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
190:      Solve nonlinear system
191:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
192:   SNESSolve(snes,user.psi,&its);

194:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
195:      Explicitly check norm of the residual of the solution
196:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
197:   SNESDAFormFunction(snes,user.psi,user.r,(void *)&user);
198:   VecNorm(user.r,NORM_MAX,&fnorm);
199:   PetscPrintf(PETSC_COMM_WORLD,"Number of Newton iterations = %d fnorm %gn",its,fnorm);

201:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 
202:      Output the solution vector
203:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
204:   {
205:     PetscViewer view_out;
206:     PetscViewerBinaryOpen(PETSC_COMM_WORLD,"psi.binary",PETSC_BINARY_CREATE,&view_out);
207:     VecView(user.psi,view_out);
208:     PetscViewerDestroy(view_out);
209:     PetscViewerASCIIOpen(PETSC_COMM_WORLD,"psi.out",&view_out);
210:     VecView(user.psi,view_out);
211:     PetscViewerDestroy(view_out);
212:   }

214:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
215:      Free work space.  All PETSc objects should be destroyed when they
216:      are no longer needed.
217:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

219:   if (user.A != user.J) {
220:     MatDestroy(user.A);
221:   }
222:   MatDestroy(user.J);
223:   if (matfdcoloring) {
224:     MatFDColoringDestroy(matfdcoloring);
225:   }
226:   VecDestroy(user.psi);
227:   VecDestroy(user.r);
228:   SNESDestroy(snes);
229:   DADestroy(user.da);
230:   PetscFinalize();

232:   return(0);
233: }
234: /* ------------------------------------------------------------------- */
235: #undef __FUNCT__
237: /* 
238:    FormInitialGuess - Forms initial approximation.

240:    Input Parameters:
241:    user - user-defined application context
242:    X - vector

244:    Output Parameter:
245:    X - vector
246:  */
247: int FormInitialGuess(AppCtx *user,Vec X)
248: {
249:   int         i,Mx,ierr,xs,xm;
250:   PetscScalar *x;

253:   DAGetInfo(user->da,PETSC_IGNORE,&Mx,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,
254:                    PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE);

256:   /*
257:      Get a pointer to vector data.
258:        - For default PETSc vectors, VecGetArray() returns a pointer to
259:          the data array.  Otherwise, the routine is implementation dependent.
260:        - You MUST call VecRestoreArray() when you no longer need access to
261:          the array.
262:   */
263:   DAVecGetArray(user->da,X,(void**)&x);

265:   /*
266:      Get local grid boundaries (for 2-dimensional DA):
267:        xs, ys   - starting grid indices (no ghost points)
268:        xm, ym   - widths of local grid (no ghost points)

270:   */
271:   DAGetCorners(user->da,&xs,PETSC_NULL,PETSC_NULL,&xm,PETSC_NULL,PETSC_NULL);

273:   /*
274:      Compute initial guess over the locally owned part of the grid
275:   */
276:   for (i=xs; i<xs+xm; i++) {
277:     x[i] = user->psi_axis + i*(user->psi_bdy - user->psi_axis)/(PetscReal)(Mx-1);
278:   }

280:   /*
281:      Restore vector
282:   */
283:   DAVecRestoreArray(user->da,X,(void**)&x);

285:   /* 
286:      Check to see if we can import an initial guess from disk
287:   */
288:   {
289:     char         filename[256];
290:     PetscTruth   flg;
291:     PetscViewer  view_in;
292:     PetscReal    fnorm;
293:     Vec          Y;
294:     PetscOptionsGetString(PETSC_NULL,"-initialGuess",filename,256,&flg);
295:     if (flg) {
296:       PetscViewerBinaryOpen(PETSC_COMM_WORLD,filename,PETSC_BINARY_RDONLY,&view_in);
297:       VecLoad(view_in,&Y);
298:       PetscViewerDestroy(view_in);
299:       VecMax(Y,PETSC_NULL,&user->psi_bdy);
300:       SNESDAFormFunction(PETSC_NULL,Y,user->r,(void*)user);
301:       VecNorm(user->r,NORM_2,&fnorm);
302:       PetscPrintf(PETSC_COMM_WORLD,"In initial guess: psi_bdy = %f, fnorm = %g.n",user->psi_bdy,fnorm);
303:       VecCopy(Y,X);
304:       VecDestroy(Y);
305:     }
306:   }

308:   return(0);
309: }
310: /* ------------------------------------------------------------------- */
311: #undef __FUNCT__
313: /* 
314:    FormFunctionLocal - Evaluates nonlinear function, F(x).
315:    
316:    Process adiC: FormFunctionLocal
317:    
318: */
319: int FormFunctionLocal(DALocalInfo *info,PetscScalar *x,PetscScalar *f,AppCtx *user)
320: {
321:   int         ierr,i,xint,xend;
322:   PetscReal   hx,dhx,r;
323:   PetscScalar u,uxx,min = 100000.0,max = -10000.0;
324:   PetscScalar psi_0=0.0, psi_a=1.0;
325: 
327:   for (i=info->xs; i<info->xs + info->xm; i++) {
328:     if (x[i] > max) max = x[i];
329:     if (x[i] < min) min = x[i];
330:   }
331:   PetscGlobalMax(&max,&psi_a,PETSC_COMM_WORLD);
332:   PetscGlobalMin(&min,&psi_0,PETSC_COMM_WORLD);

334:   hx     = 1.0/(PetscReal)(info->mx-1);  dhx    = 1.0/hx;
335: 
336:   /*
337:     Compute function over the locally owned part of the grid
338:   */
339:   if (info->xs == 0) {
340:     xint = info->xs + 1; f[0] = (4.*x[1] - x[2] - 3.*x[0])*dhx; /* f[0] = x[0] - user->psi_axis; */
341:   }
342:   else {
343:     xint = info->xs;
344:   }
345:   if ((info->xs+info->xm) == info->mx) {
346:     xend = info->mx - 1; f[info->mx-1] = -(x[info->mx-1] - user->psi_bdy)*dhx;
347:   }
348:   else {
349:     xend = info->xs + info->xm;
350:   }
351: 
352:   for (i=xint; i<xend; i++) {
353:     r       = i*hx + user->r_min;   /* r coordinate */
354:     u       = (x[i]-psi_0)/(psi_a - psi_0);
355:     uxx     = ((r+0.5*hx)*(x[i+1]-x[i]) - (r-0.5*hx)*(x[i]-x[i-1]))*dhx; /* r*nabla^2psi */
356:     f[i] = uxx  + r*GdGdPsi(r,u)*hx  + r*CurrentWire(r)*hx ;
357:   }

359:   PetscLogFlops(11*info->ym*info->xm);
360:   return(0);
361: }
362: /* ------------------------------------------------------------------- */
363: #undef __FUNCT__
365: /*
366:    FormJacobian - Evaluates Jacobian matrix.

368:    Input Parameters:
369: .  snes - the SNES context
370: .  x - input vector
371: .  ptr - optional user-defined context, as set by SNESSetJacobian()

373:    Output Parameters:
374: .  A - Jacobian matrix
375: .  B - optionally different preconditioning matrix
376: .  flag - flag indicating matrix structure

378: */
379: int FormJacobian(SNES snes,Vec X,Mat *J,Mat *B,MatStructure *flag,void *ptr)
380: {
381:   AppCtx       *user = (AppCtx*)ptr;  /* user-defined application context */
382:   Mat          jac = *B;                /* Jacobian matrix */
383:   Vec          localX;
384:   int          ierr,i,j;
385:   int          col[6],row;
386:   int          xs,xm,Mx,xint,xend;
387:   PetscScalar  v[6],hx,dhx,*x;
388:   PetscReal    r, u, psi_0=0.0, psi_a=1.0;
389:   int          imin, imax;

392:   MatZeroEntries(*B);

394:   DAGetLocalVector(user->da,&localX);
395:   DAGetInfo(user->da,PETSC_IGNORE,&Mx,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,
396:                    PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE);

398:   hx     = 1.0/(PetscReal)(Mx-1);  dhx    = 1.0/hx;

400:   imin = 0; imax = Mx-1;
401:   VecMin(X,&imin,&psi_0);
402:   VecMax(X,&imax,&psi_a);
403:   PetscPrintf(PETSC_COMM_WORLD,"psi_0(%d)=%g, psi_a(%d)=%g.n",imin,psi_0,imax,psi_a);

405:   /*
406:      Scatter ghost points to local vector, using the 2-step process
407:         DAGlobalToLocalBegin(), DAGlobalToLocalEnd().
408:      By placing code between these two statements, computations can be
409:      done while messages are in transition.
410:   */
411:   DAGlobalToLocalBegin(user->da,X,INSERT_VALUES,localX);
412:   DAGlobalToLocalEnd(user->da,X,INSERT_VALUES,localX);

414:   /*
415:      Get pointer to vector data
416:   */
417:   DAVecGetArray(user->da,localX,(void**)&x);

419:   /*
420:      Get local grid boundaries
421:   */
422:   DAGetCorners(user->da,&xs,PETSC_NULL,PETSC_NULL,&xm,PETSC_NULL,PETSC_NULL);

424:   /* 
425:      Compute entries for the locally owned part of the Jacobian.
426:       - Currently, all PETSc parallel matrix formats are partitioned by
427:         contiguous chunks of rows across the processors. 
428:       - Each processor needs to insert only elements that it owns
429:         locally (but any non-local elements will be sent to the
430:         appropriate processor during matrix assembly). 
431:       - Here, we set all entries for a particular row at once.
432:       - We can set matrix entries either using either
433:         MatSetValuesLocal() or MatSetValues(), as discussed above.
434:   */
435:   if (xs == 0) {
436:     xint = xs + 1; /* f[0] = 4.*x[1] - x[2] - 3.*x[0]; */
437:     row  = 0;     /* first row */
438:     v[0] = -3.0*dhx;                                              col[0]=row;
439:     v[1] = 4.0*dhx;                                               col[1]=row+1;
440:     v[2] = -1.0*dhx;                                              col[2]=row+2;
441:     MatSetValues(jac,1,&row,3,col,v,ADD_VALUES);
442:   }
443:   else {
444:     xint = xs;
445:   }
446:   if ((xs+xm) == Mx) {
447:     xend = Mx - 1;   /* f[Mx-1] = x[Mx-1] - user->psi_bdy; */
448:     row  = Mx - 1;  /* last row */
449:     v[0] = -1.0*dhx;
450:     MatSetValue(jac,row,row,v[0],ADD_VALUES);
451:   }
452:   else {
453:     xend = xs + xm;
454:   }

456:   for (i=xint; i<xend; i++) {
457:     r       = i*hx + user->r_min;   /* r coordinate */
458:     u       = (x[i]-psi_0)/(psi_a - psi_0);
459:     /* uxx     = ((r+0.5*hx)*(x[i+1]-x[i]) - (r-0.5*hx)*(x[i]-x[i-1]))*dhx*dhx; */ /* r*nabla^2psi */
460:     row  = i;
461:     v[0] = (r-0.5*hx)*dhx;                                                              col[0] = i-1;
462:     v[1] = -2.*r*dhx + hx*r*GdGdPsiPrime(r,u)/(psi_a - psi_0);                          col[1] = i;
463:     v[2] = (r+0.5*hx)*dhx;                                                              col[2] = i+1;
464:     v[3] = hx*r*GdGdPsiPrime(r,u)*(x[i] - psi_a)/((psi_a - psi_0)*(psi_a - psi_0));     col[3] = imin;
465:     v[4] = hx*r*GdGdPsiPrime(r,u)*(psi_0 - x[i])/((psi_a - psi_0)*(psi_a - psi_0));     col[4] = imax;
466:     MatSetValues(jac,1,&row,5,col,v,ADD_VALUES);
467:   }
468: 
469:   DAVecRestoreArray(user->da,localX,(void**)&x);
470:   DARestoreLocalVector(user->da,&localX);

472:   /* 
473:      Assemble matrix, using the 2-step process:
474:        MatAssemblyBegin(), MatAssemblyEnd().
475:   */
476:   MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY);
477:   MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY);

479:   /*
480:      Normally since the matrix has already been assembled above; this
481:      would do nothing. But in the matrix free mode -snes_mf_operator
482:      this tells the "matrix-free" matrix that a new linear system solve
483:      is about to be done.
484:   */
485:   MatAssemblyBegin(*J,MAT_FINAL_ASSEMBLY);
486:   MatAssemblyEnd(*J,MAT_FINAL_ASSEMBLY);

488:   /*
489:      Set flag to indicate that the Jacobian matrix retains an identical
490:      nonzero structure throughout all nonlinear iterations (although the
491:      values of the entries change). Thus, we can save some work in setting
492:      up the preconditioner (e.g., no need to redo symbolic factorization for
493:      ILU/ICC preconditioners).
494:       - If the nonzero structure of the matrix is different during
495:         successive linear solves, then the flag DIFFERENT_NONZERO_PATTERN
496:         must be used instead.  If you are unsure whether the matrix
497:         structure has changed or not, use the flag DIFFERENT_NONZERO_PATTERN.
498:       - Caution:  If you specify SAME_NONZERO_PATTERN, PETSc
499:         believes your assertion and does not check the structure
500:         of the matrix.  If you erroneously claim that the structure
501:         is the same when it actually is not, the new preconditioner
502:         will not function correctly.  Thus, use this optimization
503:         feature with caution!
504:   */
505:   *flag = SAME_NONZERO_PATTERN;

507:   /*
508:      Tell the matrix we will never add a new nonzero location to the
509:      matrix. If we do, it will generate an error.
510:   */

512:   return(0);
513: }