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: }