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