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 points\n\
6: -psi_axis <axis> \n\
7: -r_min <min> \n\
8: -param <param> \n\n";
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*);
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: MatCreate(PETSC_COMM_WORLD,PETSC_DECIDE,PETSC_DECIDE,grids,grids,&user.J);
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: PetscOptionsGetLogical(PETSC_NULL,"-adicmf_jacobian",&adicmf_jacobian,0);
155: #if defined(PETSC_HAVE_ADIC) && !defined(PETSC_USE_COMPLEX) && !defined(PETSC_USE_SINGLE)
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_LOCAL,&iscoloring);
167: MatFDColoringCreate(user.J,iscoloring,&matfdcoloring);
168: ISColoringDestroy(iscoloring);
169: MatFDColoringSetFunction(matfdcoloring,(int (*)(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,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",PETSC_FILE_CREATE,&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: int FormInitialGuess(AppCtx *user,Vec X)
250: {
251: int i,Mx,ierr,xs,xm;
252: PetscScalar *x;
255: DAGetInfo(user->da,PETSC_IGNORE,&Mx,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,
256: PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE);
258: /*
259: Get a pointer to vector data.
260: - For default PETSc vectors, VecGetArray() returns a pointer to
261: the data array. Otherwise, the routine is implementation dependent.
262: - You MUST call VecRestoreArray() when you no longer need access to
263: the array.
264: */
265: DAVecGetArray(user->da,X,&x);
267: /*
268: Get local grid boundaries (for 2-dimensional DA):
269: xs, ys - starting grid indices (no ghost points)
270: xm, ym - widths of local grid (no ghost points)
272: */
273: DAGetCorners(user->da,&xs,PETSC_NULL,PETSC_NULL,&xm,PETSC_NULL,PETSC_NULL);
275: /*
276: Compute initial guess over the locally owned part of the grid
277: */
278: for (i=xs; i<xs+xm; i++) {
279: x[i] = user->psi_axis + i*(user->psi_bdy - user->psi_axis)/(PetscReal)(Mx-1);
280: }
282: /*
283: Restore vector
284: */
285: DAVecRestoreArray(user->da,X,&x);
287: /*
288: Check to see if we can import an initial guess from disk
289: */
290: {
291: char filename[256];
292: PetscTruth flg;
293: PetscViewer view_in;
294: PetscReal fnorm;
295: Vec Y;
296: PetscOptionsGetString(PETSC_NULL,"-initialGuess",filename,256,&flg);
297: if (flg) {
298: PetscViewerBinaryOpen(PETSC_COMM_WORLD,filename,PETSC_FILE_RDONLY,&view_in);
299: VecLoad(view_in,PETSC_NULL,&Y);
300: PetscViewerDestroy(view_in);
301: VecMax(Y,PETSC_NULL,&user->psi_bdy);
302: SNESDAFormFunction(PETSC_NULL,Y,user->r,(void*)user);
303: VecNorm(user->r,NORM_2,&fnorm);
304: PetscPrintf(PETSC_COMM_WORLD,"In initial guess: psi_bdy = %f, fnorm = %g.\n",user->psi_bdy,fnorm);
305: VecCopy(Y,X);
306: VecDestroy(Y);
307: }
308: }
310: return(0);
311: }
312: /* ------------------------------------------------------------------- */
315: /*
316: FormFunctionLocal - Evaluates nonlinear function, F(x).
317:
318: Process adiC(36): FormFunctionLocal
319:
320: */
321: int FormFunctionLocal(DALocalInfo *info,PetscScalar *x,PetscScalar *f,AppCtx *user)
322: {
323: int ierr,i,xint,xend;
324: PetscReal hx,dhx,r;
325: PetscScalar u,uxx,min = 100000.0,max = -10000.0;
326: PetscScalar psi_0=0.0, psi_a=1.0;
327:
329: for (i=info->xs; i<info->xs + info->xm; i++) {
330: if (x[i] > max) max = x[i];
331: if (x[i] < min) min = x[i];
332: }
333: PetscGlobalMax(&max,&psi_a,PETSC_COMM_WORLD);
334: PetscGlobalMin(&min,&psi_0,PETSC_COMM_WORLD);
336: hx = 1.0/(PetscReal)(info->mx-1); dhx = 1.0/hx;
337:
338: /*
339: Compute function over the locally owned part of the grid
340: */
341: if (info->xs == 0) {
342: xint = info->xs + 1; f[0] = (4.*x[1] - x[2] - 3.*x[0])*dhx; /* f[0] = x[0] - user->psi_axis; */
343: }
344: else {
345: xint = info->xs;
346: }
347: if ((info->xs+info->xm) == info->mx) {
348: xend = info->mx - 1; f[info->mx-1] = -(x[info->mx-1] - user->psi_bdy)*dhx;
349: }
350: else {
351: xend = info->xs + info->xm;
352: }
353:
354: for (i=xint; i<xend; i++) {
355: r = i*hx + user->r_min; /* r coordinate */
356: u = (x[i]-psi_0)/(psi_a - psi_0);
357: uxx = ((r+0.5*hx)*(x[i+1]-x[i]) - (r-0.5*hx)*(x[i]-x[i-1]))*dhx; /* r*nabla^2\psi */
358: f[i] = uxx + r*GdGdPsi(r,u)*hx + r*CurrentWire(r)*hx ;
359: }
361: PetscLogFlops(11*info->ym*info->xm);
362: return(0);
363: }
364: /* ------------------------------------------------------------------- */
367: /*
368: FormJacobian - Evaluates Jacobian matrix.
370: Input Parameters:
371: . snes - the SNES context
372: . x - input vector
373: . ptr - optional user-defined context, as set by SNESSetJacobian()
375: Output Parameters:
376: . A - Jacobian matrix
377: . B - optionally different preconditioning matrix
378: . flag - flag indicating matrix structure
380: */
381: int FormJacobian(SNES snes,Vec X,Mat *J,Mat *B,MatStructure *flag,void *ptr)
382: {
383: AppCtx *user = (AppCtx*)ptr; /* user-defined application context */
384: Mat jac = *B; /* Jacobian matrix */
385: Vec localX;
386: int ierr,i,j;
387: int col[6],row;
388: int xs,xm,Mx,xint,xend;
389: PetscScalar v[6],hx,dhx,*x;
390: PetscReal r, u, psi_0=0.0, psi_a=1.0;
391: int imin, imax;
394: MatZeroEntries(*B);
396: DAGetLocalVector(user->da,&localX);
397: DAGetInfo(user->da,PETSC_IGNORE,&Mx,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,
398: PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE);
400: hx = 1.0/(PetscReal)(Mx-1); dhx = 1.0/hx;
402: imin = 0; imax = Mx-1;
403: VecMin(X,&imin,&psi_0);
404: VecMax(X,&imax,&psi_a);
405: PetscPrintf(PETSC_COMM_WORLD,"psi_0(%d)=%g, psi_a(%d)=%g.\n",imin,psi_0,imax,psi_a);
407: /*
408: Scatter ghost points to local vector, using the 2-step process
409: DAGlobalToLocalBegin(), DAGlobalToLocalEnd().
410: By placing code between these two statements, computations can be
411: done while messages are in transition.
412: */
413: DAGlobalToLocalBegin(user->da,X,INSERT_VALUES,localX);
414: DAGlobalToLocalEnd(user->da,X,INSERT_VALUES,localX);
416: /*
417: Get pointer to vector data
418: */
419: DAVecGetArray(user->da,localX,&x);
421: /*
422: Get local grid boundaries
423: */
424: DAGetCorners(user->da,&xs,PETSC_NULL,PETSC_NULL,&xm,PETSC_NULL,PETSC_NULL);
426: /*
427: Compute entries for the locally owned part of the Jacobian.
428: - Currently, all PETSc parallel matrix formats are partitioned by
429: contiguous chunks of rows across the processors.
430: - Each processor needs to insert only elements that it owns
431: locally (but any non-local elements will be sent to the
432: appropriate processor during matrix assembly).
433: - Here, we set all entries for a particular row at once.
434: - We can set matrix entries either using either
435: MatSetValuesLocal() or MatSetValues(), as discussed above.
436: */
437: if (xs == 0) {
438: xint = xs + 1; /* f[0] = 4.*x[1] - x[2] - 3.*x[0]; */
439: row = 0; /* first row */
440: v[0] = -3.0*dhx; col[0]=row;
441: v[1] = 4.0*dhx; col[1]=row+1;
442: v[2] = -1.0*dhx; col[2]=row+2;
443: MatSetValues(jac,1,&row,3,col,v,ADD_VALUES);
444: }
445: else {
446: xint = xs;
447: }
448: if ((xs+xm) == Mx) {
449: xend = Mx - 1; /* f[Mx-1] = x[Mx-1] - user->psi_bdy; */
450: row = Mx - 1; /* last row */
451: v[0] = -1.0*dhx;
452: MatSetValue(jac,row,row,v[0],ADD_VALUES);
453: }
454: else {
455: xend = xs + xm;
456: }
458: for (i=xint; i<xend; i++) {
459: r = i*hx + user->r_min; /* r coordinate */
460: u = (x[i]-psi_0)/(psi_a - psi_0);
461: /* 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 */
462: row = i;
463: v[0] = (r-0.5*hx)*dhx; col[0] = i-1;
464: v[1] = -2.*r*dhx + hx*r*GdGdPsiPrime(r,u)/(psi_a - psi_0); col[1] = i;
465: v[2] = (r+0.5*hx)*dhx; col[2] = i+1;
466: v[3] = hx*r*GdGdPsiPrime(r,u)*(x[i] - psi_a)/((psi_a - psi_0)*(psi_a - psi_0)); col[3] = imin;
467: v[4] = hx*r*GdGdPsiPrime(r,u)*(psi_0 - x[i])/((psi_a - psi_0)*(psi_a - psi_0)); col[4] = imax;
468: MatSetValues(jac,1,&row,5,col,v,ADD_VALUES);
469: }
470:
471: DAVecRestoreArray(user->da,localX,&x);
472: DARestoreLocalVector(user->da,&localX);
474: /*
475: Assemble matrix, using the 2-step process:
476: MatAssemblyBegin(), MatAssemblyEnd().
477: */
478: MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY);
479: MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY);
481: /*
482: Normally since the matrix has already been assembled above; this
483: would do nothing. But in the matrix free mode -snes_mf_operator
484: this tells the "matrix-free" matrix that a new linear system solve
485: is about to be done.
486: */
487: MatAssemblyBegin(*J,MAT_FINAL_ASSEMBLY);
488: MatAssemblyEnd(*J,MAT_FINAL_ASSEMBLY);
490: /*
491: Set flag to indicate that the Jacobian matrix retains an identical
492: nonzero structure throughout all nonlinear iterations (although the
493: values of the entries change). Thus, we can save some work in setting
494: up the preconditioner (e.g., no need to redo symbolic factorization for
495: ILU/ICC preconditioners).
496: - If the nonzero structure of the matrix is different during
497: successive linear solves, then the flag DIFFERENT_NONZERO_PATTERN
498: must be used instead. If you are unsure whether the matrix
499: structure has changed or not, use the flag DIFFERENT_NONZERO_PATTERN.
500: - Caution: If you specify SAME_NONZERO_PATTERN, PETSc
501: believes your assertion and does not check the structure
502: of the matrix. If you erroneously claim that the structure
503: is the same when it actually is not, the new preconditioner
504: will not function correctly. Thus, use this optimization
505: feature with caution!
506: */
507: *flag = SAME_NONZERO_PATTERN;
509: /*
510: Tell the matrix we will never add a new nonzero location to the
511: matrix. If we do, it will generate an error.
512: */
514: return(0);
515: }