Actual source code: ex10.c
1: /*$Id: ex10.c,v 1.29 2001/09/11 16:34:10 bsmith Exp $*/
3: /*
4: Program usage: mpirun -np <procs> usg [-help] [all PETSc options]
5: */
7: #if !defined(PETSC_USE_COMPLEX)
9: static char help[] = "An Unstructured Grid Example.n
10: This example demonstrates how to solve a nonlinear system in paralleln
11: with SNES for an unstructured mesh. The mesh and partitioning informationn
12: is read in an application defined ordering,which is later transformedn
13: into another convenient ordering (called the local ordering). The localn
14: ordering, apart from being efficient on cpu cycles and memory, allowsn
15: the use of the SPMD model of parallel programming. After partitioningn
16: is done, scatters are created between local (sequential)and globaln
17: (distributed) vectors. Finally, we set up the nonlinear solver contextn
18: in the usual way as a structured grid (seen
19: petsc/src/snes/examples/tutorials/ex5.c).n
20: The command line options include:n
21: -vert <Nv>, where Nv is the global number of nodesn
22: -elem <Ne>, where Ne is the global number of elementsn
23: -nl_par <lambda>, where lambda is the multiplier for the non linear term (u*u) termn
24: -lin_par <alpha>, where alpha is the multiplier for the linear term (u) n";
26: /*T
27: Concepts: SNES^unstructured grid
28: Concepts: AO^application to PETSc ordering or vice versa;
29: Concepts: VecScatter^using vector scatter operations;
30: Processors: n
31: T*/
33: /* ------------------------------------------------------------------------
35: PDE Solved : L(u) + lambda*u*u + alpha*u = 0 where L(u) is the Laplacian.
37: The Laplacian is approximated in the following way: each edge is given a weight
38: of one meaning that the diagonal term will have the weight equal to the degree
39: of a node. The off diagonal terms will get a weight of -1.
41: -----------------------------------------------------------------------*/
43: /*
44: Include petscao.h so that we can use AO (Application Ordering) object's services.
45: Include "petscsnes.h" so that we can use SNES solvers. Note that this
46: file automatically includes:
47: petsc.h - base PETSc routines petscvec.h - vectors
48: petscsys.h - system routines petscmat.h - matrices
49: petscis.h - index sets petscksp.h - Krylov subspace methods
50: petscviewer.h - viewers petscpc.h - preconditioners
51: petscsles.h - linear solvers
52: */
53: #include "petscao.h"
54: #include "petscsnes.h"
57: #define MAX_ELEM 500 /* Maximum number of elements */
58: #define MAX_VERT 100 /* Maximum number of vertices */
59: #define MAX_VERT_ELEM 3 /* Vertices per element */
61: /*
62: Application-defined context for problem specific data
63: */
64: typedef struct {
65: int Nvglobal,Nvlocal; /* global and local number of vertices */
66: int Neglobal,Nelocal; /* global and local number of vertices */
67: int AdjM[MAX_VERT][50]; /* adjacency list of a vertex */
68: int itot[MAX_VERT]; /* total number of neighbors for a vertex */
69: int icv[MAX_ELEM][MAX_VERT_ELEM]; /* vertices belonging to an element */
70: int v2p[MAX_VERT]; /* processor number for a vertex */
71: int *locInd,*gloInd; /* local and global orderings for a node */
72: Vec localX,localF; /* local solution (u) and f(u) vectors */
73: PetscReal non_lin_param; /* nonlinear parameter for the PDE */
74: PetscReal lin_param; /* linear parameter for the PDE */
75: VecScatter scatter; /* scatter context for the local and
76: distributed vectors */
77: } AppCtx;
79: /*
80: User-defined routines
81: */
82: int FormJacobian(SNES,Vec,Mat*,Mat*,MatStructure*,void*),
83: FormFunction(SNES,Vec,Vec,void*),
84: FormInitialGuess(AppCtx*,Vec);
86: #undef __FUNCT__
88: int main(int argc,char **argv)
89: {
90: SNES snes; /* SNES context */
91: SNESType type = SNESLS; /* default nonlinear solution method */
92: Vec x,r; /* solution, residual vectors */
93: Mat Jac; /* Jacobian matrix */
94: AppCtx user; /* user-defined application context */
95: AO ao; /* Application Ordering object */
96: IS isglobal,islocal; /* global and local index sets */
97: int rank,size; /* rank of a process, number of processors */
98: int rstart; /* starting index of PETSc ordering for a processor */
99: int nfails; /* number of unsuccessful Newton steps */
100: int bs = 1; /* block size for multicomponent systems */
101: int nvertices; /* number of local plus ghost nodes of a processor */
102: int *pordering; /* PETSc ordering */
103: int *vertices; /* list of all vertices (incl. ghost ones)
104: on a processor */
105: int *verticesmask,*svertices;
106: int *tmp;
107: int i,j,jstart,inode,nb,nbrs,Nvneighborstotal = 0;
108: int ierr,its,N;
109: PetscScalar *xx;
110: char str[256],form[256],part_name[256];
111: FILE *fptr,*fptr1;
112: ISLocalToGlobalMapping isl2g;
113: #if defined (UNUSED_VARIABLES)
114: PetscDraw draw; /* drawing context */
115: PetscScalar *ff,*gg;
116: PetscReal tiny = 1.0e-10,zero = 0.0,one = 1.0,big = 1.0e+10;
117: int *tmp1,*tmp2;
118: #endif
119: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
120: Initialize program
121: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
123: PetscInitialize(&argc,&argv,"options.inf",help);
124: MPI_Comm_rank(MPI_COMM_WORLD,&rank);
125: MPI_Comm_size(MPI_COMM_WORLD,&size);
127: /* The current input file options.inf is for 2 proc run only */
128: if (size != 2) SETERRQ(1,"This Example currently runs on 2 procs only.");
130: /*
131: Initialize problem parameters
132: */
133: user.Nvglobal = 16; /*Global # of vertices */
134: user.Neglobal = 18; /*Global # of elements */
135: PetscOptionsGetInt(PETSC_NULL,"-vert",&user.Nvglobal,PETSC_NULL);
136: PetscOptionsGetInt(PETSC_NULL,"-elem",&user.Neglobal,PETSC_NULL);
137: user.non_lin_param = 0.06;
138: PetscOptionsGetReal(PETSC_NULL,"-nl_par",&user.non_lin_param,PETSC_NULL);
139: user.lin_param = -1.0;
140: PetscOptionsGetReal(PETSC_NULL,"-lin_par",&user.lin_param,PETSC_NULL);
141: user.Nvlocal = 0;
142: user.Nelocal = 0;
144: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
145: Read the mesh and partitioning information
146: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
147:
148: /*
149: Read the mesh and partitioning information from 'adj.in'.
150: The file format is as follows.
151: For each line the first entry is the processor rank where the
152: current node belongs. The second entry is the number of
153: neighbors of a node. The rest of the line is the adjacency
154: list of a node. Currently this file is set up to work on two
155: processors.
157: This is not a very good example of reading input. In the future,
158: we will put an example that shows the style that should be
159: used in a real application, where partitioning will be done
160: dynamically by calling partitioning routines (at present, we have
161: a ready interface to ParMeTiS).
162: */
163: fptr = fopen("adj.in","r");
164: if (!fptr) {
165: SETERRQ(0,"Could not open adj.in")
166: }
167:
168: /*
169: Each processor writes to the file output.<rank> where rank is the
170: processor's rank.
171: */
172: sprintf(part_name,"output.%d",rank);
173: fptr1 = fopen(part_name,"w");
174: if (!fptr1) {
175: SETERRQ(0,"Could no open output file");
176: }
177: PetscMalloc(user.Nvglobal*sizeof(int),&user.gloInd);
178: fprintf(fptr1,"Rank is %dn",rank);
179: for (inode = 0; inode < user.Nvglobal; inode++) {
180: fgets(str,256,fptr);
181: sscanf(str,"%d",&user.v2p[inode]);
182: if (user.v2p[inode] == rank) {
183: fprintf(fptr1,"Node %d belongs to processor %dn",inode,user.v2p[inode]);
184: user.gloInd[user.Nvlocal] = inode;
185: sscanf(str,"%*d %d",&nbrs);
186: fprintf(fptr1,"Number of neighbors for the vertex %d is %dn",inode,nbrs);
187: user.itot[user.Nvlocal] = nbrs;
188: Nvneighborstotal += nbrs;
189: for (i = 0; i < user.itot[user.Nvlocal]; i++){
190: form[0]='0';
191: for (j=0; j < i+2; j++){
192: PetscStrcat(form,"%*d ");
193: }
194: PetscStrcat(form,"%d");
195: sscanf(str,form,&user.AdjM[user.Nvlocal][i]);
196: fprintf(fptr1,"%d ",user.AdjM[user.Nvlocal][i]);
197: }
198: fprintf(fptr1,"n");
199: user.Nvlocal++;
200: }
201: }
202: fprintf(fptr1,"Total # of Local Vertices is %d n",user.Nvlocal);
204: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
205: Create different orderings
206: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
208: /*
209: Create the local ordering list for vertices. First a list using the PETSc global
210: ordering is created. Then we use the AO object to get the PETSc-to-application and
211: application-to-PETSc mappings. Each vertex also gets a local index (stored in the
212: locInd array).
213: */
214: MPI_Scan(&user.Nvlocal,&rstart,1,MPI_INT,MPI_SUM,MPI_COMM_WORLD);
215: rstart -= user.Nvlocal;
216: PetscMalloc(user.Nvlocal*sizeof(int),&pordering);
218: for (i=0; i < user.Nvlocal; i++) {
219: pordering[i] = rstart + i;
220: }
222: /*
223: Create the AO object
224: */
225: AOCreateBasic(MPI_COMM_WORLD,user.Nvlocal,user.gloInd,pordering,&ao);
226: PetscFree(pordering);
227:
228: /*
229: Keep the global indices for later use
230: */
231: PetscMalloc(user.Nvlocal*sizeof(int),&user.locInd);
232: PetscMalloc(Nvneighborstotal*sizeof(int),&tmp);
233:
234: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
235: Demonstrate the use of AO functionality
236: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
238: fprintf(fptr1,"Before AOApplicationToPetsc, local indices are : n");
239: for (i=0; i < user.Nvlocal; i++) {
240: fprintf(fptr1," %d ",user.gloInd[i]);
241: user.locInd[i] = user.gloInd[i];
242: }
243: fprintf(fptr1,"n");
244: jstart = 0;
245: for (i=0; i < user.Nvlocal; i++) {
246: fprintf(fptr1,"Neghbors of local vertex %d are : ",user.gloInd[i]);
247: for (j=0; j < user.itot[i]; j++) {
248: fprintf(fptr1,"%d ",user.AdjM[i][j]);
249: tmp[j + jstart] = user.AdjM[i][j];
250: }
251: jstart += user.itot[i];
252: fprintf(fptr1,"n");
253: }
255: /*
256: Now map the vlocal and neighbor lists to the PETSc ordering
257: */
258: AOApplicationToPetsc(ao,user.Nvlocal,user.locInd);
259: AOApplicationToPetsc(ao,Nvneighborstotal,tmp);
260:
261: fprintf(fptr1,"After AOApplicationToPetsc, local indices are : n");
262: for (i=0; i < user.Nvlocal; i++) {
263: fprintf(fptr1," %d ",user.locInd[i]);
264: }
265: fprintf(fptr1,"n");
267: jstart = 0;
268: for (i=0; i < user.Nvlocal; i++) {
269: fprintf(fptr1,"Neghbors of local vertex %d are : ",user.locInd[i]);
270: for (j=0; j < user.itot[i]; j++) {
271: user.AdjM[i][j] = tmp[j+jstart];
272: fprintf(fptr1,"%d ",user.AdjM[i][j]);
273: }
274: jstart += user.itot[i];
275: fprintf(fptr1,"n");
276: }
278: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
279: Extract the ghost vertex information for each processor
280: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
281: /*
282: Next, we need to generate a list of vertices required for this processor
283: and a local numbering scheme for all vertices required on this processor.
284: vertices - integer array of all vertices needed on this processor in PETSc
285: global numbering; this list consists of first the "locally owned"
286: vertices followed by the ghost vertices.
287: verticesmask - integer array that for each global vertex lists its local
288: vertex number (in vertices) + 1. If the global vertex is not
289: represented on this processor, then the corresponding
290: entry in verticesmask is zero
291:
292: Note: vertices and verticesmask are both Nvglobal in length; this may
293: sound terribly non-scalable, but in fact is not so bad for a reasonable
294: number of processors. Importantly, it allows us to use NO SEARCHING
295: in setting up the data structures.
296: */
297: ierr = PetscMalloc(user.Nvglobal*sizeof(int),&vertices);
298: ierr = PetscMalloc(user.Nvglobal*sizeof(int),&verticesmask);
299: ierr = PetscMemzero(verticesmask,user.Nvglobal*sizeof(int));
300: nvertices = 0;
301:
302: /*
303: First load "owned vertices" into list
304: */
305: for (i=0; i < user.Nvlocal; i++) {
306: vertices[nvertices++] = user.locInd[i];
307: verticesmask[user.locInd[i]] = nvertices;
308: }
309:
310: /*
311: Now load ghost vertices into list
312: */
313: for (i=0; i < user.Nvlocal; i++) {
314: for (j=0; j < user.itot[i]; j++) {
315: nb = user.AdjM[i][j];
316: if (!verticesmask[nb]) {
317: vertices[nvertices++] = nb;
318: verticesmask[nb] = nvertices;
319: }
320: }
321: }
323: fprintf(fptr1,"n");
324: fprintf(fptr1,"The array vertices is :n");
325: for (i=0; i < nvertices; i++) {
326: fprintf(fptr1,"%d ",vertices[i]);
327: }
328: fprintf(fptr1,"n");
329:
330: /*
331: Map the vertices listed in the neighbors to the local numbering from
332: the global ordering that they contained initially.
333: */
334: fprintf(fptr1,"n");
335: fprintf(fptr1,"After mapping neighbors in the local contiguous orderingn");
336: for (i=0; i<user.Nvlocal; i++) {
337: fprintf(fptr1,"Neghbors of local vertex %d are :n",i);
338: for (j = 0; j < user.itot[i]; j++) {
339: nb = user.AdjM[i][j];
340: user.AdjM[i][j] = verticesmask[nb] - 1;
341: fprintf(fptr1,"%d ",user.AdjM[i][j]);
342: }
343: fprintf(fptr1,"n");
344: }
346: N = user.Nvglobal;
347:
348: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
349: Create vector and matrix data structures
350: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
352: /*
353: Create vector data structures
354: */
355: VecCreate(MPI_COMM_WORLD,&x);
356: VecSetSizes(x,user.Nvlocal,N);
357: VecSetFromOptions(x);
358: VecDuplicate(x,&r);
359: VecCreateSeq(MPI_COMM_SELF,bs*nvertices,&user.localX);
360: VecDuplicate(user.localX,&user.localF);
362: /*
363: Create the scatter between the global representation and the
364: local representation
365: */
366: ISCreateStride(MPI_COMM_SELF,bs*nvertices,0,1,&islocal);
367: PetscMalloc(nvertices*sizeof(int),&svertices);
368: for (i=0; i<nvertices; i++) svertices[i] = bs*vertices[i];
369: ISCreateBlock(MPI_COMM_SELF,bs,nvertices,svertices,&isglobal);
370: PetscFree(svertices);
371: VecScatterCreate(x,isglobal,user.localX,islocal,&user.scatter);
373: /*
374: Create matrix data structure; Just to keep the example simple, we have not done any
375: preallocation of memory for the matrix. In real application code with big matrices,
376: preallocation should always be done to expedite the matrix creation.
377: */
378: MatCreate(MPI_COMM_WORLD,PETSC_DECIDE,PETSC_DECIDE,N,N,&Jac);
379: MatSetFromOptions(Jac);
381: /*
382: The following routine allows us to set the matrix values in local ordering
383: */
384: ISLocalToGlobalMappingCreate(MPI_COMM_SELF,bs*nvertices,vertices,&isl2g);
385: MatSetLocalToGlobalMapping(Jac,isl2g);
387: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
388: Create nonlinear solver context
389: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
391: SNESCreate(MPI_COMM_WORLD,&snes);
392: SNESSetType(snes,type);
394: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
395: Set routines for function and Jacobian evaluation
396: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
398: FormInitialGuess(&user,x);
399: SNESSetFunction(snes,r,FormFunction,(void *)&user);
400: SNESSetJacobian(snes,Jac,Jac,FormJacobian,(void *)&user);
402: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
403: Customize nonlinear solver; set runtime options
404: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
406: SNESSetFromOptions(snes);
408: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
409: Evaluate initial guess; then solve nonlinear system
410: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
412: /*
413: Note: The user should initialize the vector, x, with the initial guess
414: for the nonlinear solver prior to calling SNESSolve(). In particular,
415: to employ an initial guess of zero, the user should explicitly set
416: this vector to zero by calling VecSet().
417: */
418: FormInitialGuess(&user,x);
420: /*
421: Print the initial guess
422: */
423: VecGetArray(x,&xx);
424: for (inode = 0; inode < user.Nvlocal; inode++)
425: fprintf(fptr1,"Initial Solution at node %d is %f n",inode,xx[inode]);
426: VecRestoreArray(x,&xx);
428: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
429: Now solve the nonlinear system
430: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
432: SNESSolve(snes,x,&its);
433: SNESGetNumberUnsuccessfulSteps(snes,&nfails);
434:
435: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
436: Print the output : solution vector and other information
437: Each processor writes to the file output.<rank> where rank is the
438: processor's rank.
439: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
441: VecGetArray(x,&xx);
442: for (inode = 0; inode < user.Nvlocal; inode++)
443: fprintf(fptr1,"Solution at node %d is %f n",inode,xx[inode]);
444: VecRestoreArray(x,&xx);
445: fclose(fptr1);
446: PetscPrintf(MPI_COMM_WORLD,"number of Newton iterations = %d, ",its);
447: PetscPrintf(MPI_COMM_WORLD,"number of unsuccessful steps = %dn",nfails);
449: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
450: Free work space. All PETSc objects should be destroyed when they
451: are no longer needed.
452: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
454: VecDestroy(x);
455: VecDestroy(r);
456: VecDestroy(user.localX);
457: VecDestroy(user.localF);
458: MatDestroy(Jac); SNESDestroy(snes);
459: /*PetscDrawDestroy(draw);*/
460: PetscFinalize();
462: return 0;
463: }
464: #undef __FUNCT__
466: /* -------------------- Form initial approximation ----------------- */
468: /*
469: FormInitialGuess - Forms initial approximation.
471: Input Parameters:
472: user - user-defined application context
473: X - vector
475: Output Parameter:
476: X - vector
477: */
478: int FormInitialGuess(AppCtx *user,Vec X)
479: {
480: int i,Nvlocal,ierr;
481: int *gloInd;
482: PetscScalar *x;
483: #if defined (UNUSED_VARIABLES)
484: PetscReal temp1,temp,hx,hy,hxdhy,hydhx,sc;
485: int Neglobal,Nvglobal,j,row;
486: PetscReal alpha,lambda;
488: Nvglobal = user->Nvglobal;
489: Neglobal = user->Neglobal;
490: lambda = user->non_lin_param;
491: alpha = user->lin_param;
492: #endif
494: Nvlocal = user->Nvlocal;
495: gloInd = user->gloInd;
497: /*
498: Get a pointer to vector data.
499: - For default PETSc vectors, VecGetArray() returns a pointer to
500: the data array. Otherwise, the routine is implementation dependent.
501: - You MUST call VecRestoreArray() when you no longer need access to
502: the array.
503: */
504: VecGetArray(X,&x);
506: /*
507: Compute initial guess over the locally owned part of the grid
508: */
509: for (i=0; i < Nvlocal; i++) {
510: x[i] = (PetscReal)gloInd[i];
511: }
513: /*
514: Restore vector
515: */
516: VecRestoreArray(X,&x);
517: return 0;
518: }
519: #undef __FUNCT__
521: /* -------------------- Evaluate Function F(x) --------------------- */
522: /*
523: FormFunction - Evaluates nonlinear function, F(x).
525: Input Parameters:
526: . snes - the SNES context
527: . X - input vector
528: . ptr - optional user-defined context, as set by SNESSetFunction()
530: Output Parameter:
531: . F - function vector
532: */
533: int FormFunction(SNES snes,Vec X,Vec F,void *ptr)
534: {
535: AppCtx *user = (AppCtx*)ptr;
536: int ierr,i,j,Nvlocal;
537: PetscReal alpha,lambda;
538: PetscScalar *x,*f;
539: VecScatter scatter;
540: Vec localX = user->localX;
541: #if defined (UNUSED_VARIABLES)
542: PetscScalar ut,ub,ul,ur,u,*g,sc,uyy,uxx;
543: PetscReal hx,hy,hxdhy,hydhx;
544: PetscReal two = 2.0,one = 1.0;
545: int Nvglobal,Neglobal,row;
546: int *gloInd;
548: Nvglobal = user->Nvglobal;
549: Neglobal = user->Neglobal;
550: gloInd = user->gloInd;
551: #endif
553: Nvlocal = user->Nvlocal;
554: lambda = user->non_lin_param;
555: alpha = user->lin_param;
556: scatter = user->scatter;
558: /*
559: PDE : L(u) + lambda*u*u +alpha*u = 0 where L(u) is the approximate Laplacian as
560: described in the beginning of this code
561:
562: First scatter the distributed vector X into local vector localX (that includes
563: values for ghost nodes. If we wish,we can put some other work between
564: VecScatterBegin() and VecScatterEnd() to overlap the communication with
565: computation.
566: */
567: VecScatterBegin(X,localX,INSERT_VALUES,SCATTER_FORWARD,scatter);
568: VecScatterEnd(X,localX,INSERT_VALUES,SCATTER_FORWARD,scatter);
570: /*
571: Get pointers to vector data
572: */
573: VecGetArray(localX,&x);
574: VecGetArray(F,&f);
576: /*
577: Now compute the f(x). As mentioned earlier, the computed Laplacian is just an
578: approximate one chosen for illustrative purpose only. Another point to notice
579: is that this is a local (completly parallel) calculation. In practical application
580: codes, function calculation time is a dominat portion of the overall execution time.
581: */
582: for (i=0; i < Nvlocal; i++) {
583: f[i] = (user->itot[i] - alpha)*x[i] - lambda*x[i]*x[i];
584: for (j = 0; j < user->itot[i]; j++) {
585: f[i] -= x[user->AdjM[i][j]];
586: }
587: }
589: /*
590: Restore vectors
591: */
592: VecRestoreArray(localX,&x);
593: VecRestoreArray(F,&f);
594: /*VecView(F,PETSC_VIEWER_STDOUT_WORLD);*/
596: return 0;
597: }
599: #undef __FUNCT__
601: /* -------------------- Evaluate Jacobian F'(x) -------------------- */
602: /*
603: FormJacobian - Evaluates Jacobian matrix.
605: Input Parameters:
606: . snes - the SNES context
607: . X - input vector
608: . ptr - optional user-defined context, as set by SNESSetJacobian()
610: Output Parameters:
611: . A - Jacobian matrix
612: . B - optionally different preconditioning matrix
613: . flag - flag indicating matrix structure
615: */
616: int FormJacobian(SNES snes,Vec X,Mat *J,Mat *B,MatStructure *flag,void *ptr)
617: {
618: AppCtx *user = (AppCtx*)ptr;
619: Mat jac = *B;
620: int i,j,Nvlocal,col[50],ierr;
621: PetscScalar alpha,lambda,value[50];
622: Vec localX = user->localX;
623: VecScatter scatter;
624: PetscScalar *x;
625: #if defined (UNUSED_VARIABLES)
626: PetscScalar two = 2.0,one = 1.0;
627: int row,Nvglobal,Neglobal;
628: int *gloInd;
630: Nvglobal = user->Nvglobal;
631: Neglobal = user->Neglobal;
632: gloInd = user->gloInd;
633: #endif
634:
635: /*printf("Entering into FormJacobian n");*/
636: Nvlocal = user->Nvlocal;
637: lambda = user->non_lin_param;
638: alpha = user->lin_param;
639: scatter = user->scatter;
641: /*
642: PDE : L(u) + lambda*u*u +alpha*u = 0 where L(u) is the approximate Laplacian as
643: described in the beginning of this code
644:
645: First scatter the distributed vector X into local vector localX (that includes
646: values for ghost nodes. If we wish, we can put some other work between
647: VecScatterBegin() and VecScatterEnd() to overlap the communication with
648: computation.
649: */
650: VecScatterBegin(X,localX,INSERT_VALUES,SCATTER_FORWARD,scatter);
651: VecScatterEnd(X,localX,INSERT_VALUES,SCATTER_FORWARD,scatter);
652:
653: /*
654: Get pointer to vector data
655: */
656: VecGetArray(localX,&x);
658: for (i=0; i < Nvlocal; i++) {
659: col[0] = i;
660: value[0] = user->itot[i] - 2.0*lambda*x[i] - alpha;
661: for (j = 0; j < user->itot[i]; j++) {
662: col[j+1] = user->AdjM[i][j];
663: value[j+1] = -1.0;
664: }
666: /*
667: Set the matrix values in the local ordering. Note that in order to use this
668: feature we must call the routine MatSetLocalToGlobalMapping() after the
669: matrix has been created.
670: */
671: MatSetValuesLocal(jac,1,&i,1+user->itot[i],col,value,INSERT_VALUES);
672: }
674: /*
675: Assemble matrix, using the 2-step process:
676: MatAssemblyBegin(), MatAssemblyEnd().
677: Between these two calls, the pointer to vector data has been restored to
678: demonstrate the use of overlapping communicationn with computation.
679: */
680: MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY);
681: VecRestoreArray(localX,&x);
682: MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY);
684: /*
685: Set flag to indicate that the Jacobian matrix retains an identical
686: nonzero structure throughout all nonlinear iterations (although the
687: values of the entries change). Thus, we can save some work in setting
688: up the preconditioner (e.g., no need to redo symbolic factorization for
689: ILU/ICC preconditioners).
690: - If the nonzero structure of the matrix is different during
691: successive linear solves, then the flag DIFFERENT_NONZERO_PATTERN
692: must be used instead. If you are unsure whether the matrix
693: structure has changed or not, use the flag DIFFERENT_NONZERO_PATTERN.
694: - Caution: If you specify SAME_NONZERO_PATTERN, PETSc
695: believes your assertion and does not check the structure
696: of the matrix. If you erroneously claim that the structure
697: is the same when it actually is not, the new preconditioner
698: will not function correctly. Thus, use this optimization
699: feature with caution!
700: */
701: *flag = SAME_NONZERO_PATTERN;
703: /*
704: Tell the matrix we will never add a new nonzero location to the
705: matrix. If we do, it will generate an error.
706: */
707: MatSetOption(jac,MAT_NEW_NONZERO_LOCATION_ERR);
708: /* MatView(jac,PETSC_VIEWER_STDOUT_SELF); */
709: return 0;
710: }
711: #else
712: #include <stdio.h>
713: int main(int argc,char **args)
714: {
715: fprintf(stdout,"This example does not work for complex numbers.n");
716: return 0;
717: }
718: #endif