Actual source code: ex26.c

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

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

 15: /* ------------------------------------------------------------------------

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

 23:     Contributed by Xianzhu Tang, LANL

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

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

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

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

 53: /* 
 54:    User-defined routines
 55: */

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

196:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
197:      Explicitly check norm of the residual of the solution
198:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
199:   SNESDAFormFunction(snes,user.psi,user.r,(void*)&user);
200:   VecNorm(user.r,NORM_MAX,&fnorm);
201:   PetscPrintf(PETSC_COMM_WORLD,"Number of Newton iterations = %D fnorm %G\n",its,fnorm);

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

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

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

234:   return(0);
235: }
236: /* ------------------------------------------------------------------- */
239: /* 
240:    FormInitialGuess - Forms initial approximation.

242:    Input Parameters:
243:    user - user-defined application context
244:    X - vector

246:    Output Parameter:
247:    X - vector
248:  */
249: PetscErrorCode FormInitialGuess(AppCtx *user,Vec X)
250: {
252:   PetscInt       i,Mx,xs,xm;
253:   PetscScalar    *x;

256:   DAGetInfo(user->da,PETSC_IGNORE,&Mx,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,
257:                    PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE);

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

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

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

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

283:   /*
284:      Restore vector
285:   */
286:   DAVecRestoreArray(user->da,X,&x);

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

311:   return(0);
312: }
313: /* ------------------------------------------------------------------- */
316: /* 
317:    FormFunctionLocal - Evaluates nonlinear function, F(x).
318:    
319:    Process adiC(36): FormFunctionLocal
320:    
321: */
322: PetscErrorCode FormFunctionLocal(DALocalInfo *info,PetscScalar *x,PetscScalar *f,AppCtx *user)
323: {
325:   PetscInt       i,xint,xend;
326:   PetscReal      hx,dhx,r;
327:   PetscScalar    u,uxx,min = 100000.0,max = -10000.0;
328:   PetscScalar    psi_0=0.0, psi_a=1.0;
329: 
331:   for (i=info->xs; i<info->xs + info->xm; i++) {
332:     if (x[i] > max) max = x[i];
333:     if (x[i] < min) min = x[i];
334:   }
335:   PetscGlobalMax(&max,&psi_a,PETSC_COMM_WORLD);
336:   PetscGlobalMin(&min,&psi_0,PETSC_COMM_WORLD);

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

363:   PetscLogFlops(11*info->ym*info->xm);
364:   return(0);
365: }
366: /* ------------------------------------------------------------------- */
369: /*
370:    FormJacobian - Evaluates Jacobian matrix.

372:    Input Parameters:
373: .  snes - the SNES context
374: .  x - input vector
375: .  ptr - optional user-defined context, as set by SNESSetJacobian()

377:    Output Parameters:
378: .  A - Jacobian matrix
379: .  B - optionally different preconditioning matrix
380: .  flag - flag indicating matrix structure

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

395:   MatAssembled(*B,&assembled);
396:   if (assembled) {
397:     MatZeroEntries(*B);
398:   }

400:   DAGetLocalVector(user->da,&localX);
401:   DAGetInfo(user->da,PETSC_IGNORE,&Mx,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,
402:                    PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE);

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

406:   imin = 0; imax = Mx-1;
407:   VecMin(X,&imin,&psi_0);
408:   VecMax(X,&imax,&psi_a);
409:   PetscPrintf(PETSC_COMM_WORLD,"psi_0(%D)=%G, psi_a(%D)=%G.\n",imin,psi_0,imax,psi_a);

411:   /*
412:      Scatter ghost points to local vector, using the 2-step process
413:         DAGlobalToLocalBegin(), DAGlobalToLocalEnd().
414:      By placing code between these two statements, computations can be
415:      done while messages are in transition.
416:   */
417:   DAGlobalToLocalBegin(user->da,X,INSERT_VALUES,localX);
418:   DAGlobalToLocalEnd(user->da,X,INSERT_VALUES,localX);

420:   /*
421:      Get pointer to vector data
422:   */
423:   DAVecGetArray(user->da,localX,&x);

425:   /*
426:      Get local grid boundaries
427:   */
428:   DAGetCorners(user->da,&xs,PETSC_NULL,PETSC_NULL,&xm,PETSC_NULL,PETSC_NULL);

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

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

478:   /* 
479:      Assemble matrix, using the 2-step process:
480:        MatAssemblyBegin(), MatAssemblyEnd().
481:   */
482:   MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY);
483:   MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY);

485:   /*
486:      Normally since the matrix has already been assembled above; this
487:      would do nothing. But in the matrix free mode -snes_mf_operator
488:      this tells the "matrix-free" matrix that a new linear system solve
489:      is about to be done.
490:   */
491:   MatAssemblyBegin(*J,MAT_FINAL_ASSEMBLY);
492:   MatAssemblyEnd(*J,MAT_FINAL_ASSEMBLY);

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

513:   /*
514:      Tell the matrix we will never add a new nonzero location to the
515:      matrix. If we do, it will generate an error.
516:   */

518:   return(0);
519: }