Actual source code: ex13.c

petsc-master 2017-01-20
Report Typos and Errors

  2: static char help[] = "Solves a variable Poisson problem with KSP.\n\n";

  4: /*T
  5:    Concepts: KSP^basic sequential example
  6:    Concepts: KSP^Laplacian, 2d
  7:    Concepts: Laplacian, 2d
  8:    Processors: 1
  9: T*/

 11: /*
 12:   Include "petscksp.h" so that we can use KSP solvers.  Note that this file
 13:   automatically includes:
 14:      petscsys.h       - base PETSc routines   petscvec.h - vectors
 15:      petscmat.h - matrices
 16:      petscis.h     - index sets            petscksp.h - Krylov subspace methods
 17:      petscviewer.h - viewers               petscpc.h  - preconditioners
 18: */
 19:  #include <petscksp.h>

 21: /*
 22:     User-defined context that contains all the data structures used
 23:     in the linear solution process.
 24: */
 25: typedef struct {
 26:   Vec         x,b;       /* solution vector, right-hand-side vector */
 27:   Mat         A;          /* sparse matrix */
 28:   KSP         ksp;       /* linear solver context */
 29:   PetscInt    m,n;       /* grid dimensions */
 30:   PetscScalar hx2,hy2;   /* 1/(m+1)*(m+1) and 1/(n+1)*(n+1) */
 31: } UserCtx;

 33: extern PetscErrorCode UserInitializeLinearSolver(PetscInt,PetscInt,UserCtx*);
 34: extern PetscErrorCode UserFinalizeLinearSolver(UserCtx*);
 35: extern PetscErrorCode UserDoLinearSolver(PetscScalar*,UserCtx *userctx,PetscScalar *b,PetscScalar *x);

 37: int main(int argc,char **args)
 38: {
 39:   UserCtx        userctx;
 41:   PetscInt       m = 6,n = 7,t,tmax = 2,i,Ii,j,N;
 42:   PetscScalar    *userx,*rho,*solution,*userb,hx,hy,x,y;
 43:   PetscReal      enorm;

 45:   /*
 46:      Initialize the PETSc libraries
 47:   */
 48:   PetscInitialize(&argc,&args,(char*)0,help);if (ierr) return ierr;
 49:   /*
 50:      The next two lines are for testing only; these allow the user to
 51:      decide the grid size at runtime.
 52:   */
 53:   PetscOptionsGetInt(NULL,NULL,"-m",&m,NULL);
 54:   PetscOptionsGetInt(NULL,NULL,"-n",&n,NULL);

 56:   /*
 57:      Create the empty sparse matrix and linear solver data structures
 58:   */
 59:   UserInitializeLinearSolver(m,n,&userctx);
 60:   N    = m*n;

 62:   /*
 63:      Allocate arrays to hold the solution to the linear system.
 64:      This is not normally done in PETSc programs, but in this case,
 65:      since we are calling these routines from a non-PETSc program, we
 66:      would like to reuse the data structures from another code. So in
 67:      the context of a larger application these would be provided by
 68:      other (non-PETSc) parts of the application code.
 69:   */
 70:   PetscMalloc1(N,&userx);
 71:   PetscMalloc1(N,&userb);
 72:   PetscMalloc1(N,&solution);

 74:   /*
 75:       Allocate an array to hold the coefficients in the elliptic operator
 76:   */
 77:   PetscMalloc1(N,&rho);

 79:   /*
 80:      Fill up the array rho[] with the function rho(x,y) = x; fill the
 81:      right-hand-side b[] and the solution with a known problem for testing.
 82:   */
 83:   hx = 1.0/(m+1);
 84:   hy = 1.0/(n+1);
 85:   y  = hy;
 86:   Ii = 0;
 87:   for (j=0; j<n; j++) {
 88:     x = hx;
 89:     for (i=0; i<m; i++) {
 90:       rho[Ii]      = x;
 91:       solution[Ii] = PetscSinScalar(2.*PETSC_PI*x)*PetscSinScalar(2.*PETSC_PI*y);
 92:       userb[Ii]    = -2*PETSC_PI*PetscCosScalar(2*PETSC_PI *x)*PetscSinScalar(2*PETSC_PI*y) +
 93:                      8*PETSC_PI*PETSC_PI*x*PetscSinScalar(2*PETSC_PI *x)*PetscSinScalar(2*PETSC_PI*y);
 94:       x += hx;
 95:       Ii++;
 96:     }
 97:     y += hy;
 98:   }

100:   /*
101:      Loop over a bunch of timesteps, setting up and solver the linear
102:      system for each time-step.

104:      Note this is somewhat artificial. It is intended to demonstrate how
105:      one may reuse the linear solver stuff in each time-step.
106:   */
107:   for (t=0; t<tmax; t++) {
108:      UserDoLinearSolver(rho,&userctx,userb,userx);

110:     /*
111:         Compute error: Note that this could (and usually should) all be done
112:         using the PETSc vector operations. Here we demonstrate using more
113:         standard programming practices to show how they may be mixed with
114:         PETSc.
115:     */
116:     enorm = 0.0;
117:     for (i=0; i<N; i++) enorm += PetscRealPart(PetscConj(solution[i]-userx[i])*(solution[i]-userx[i]));
118:     enorm *= PetscRealPart(hx*hy);
119:     PetscPrintf(PETSC_COMM_WORLD,"m %D n %D error norm %g\n",m,n,(double)enorm);
120:   }

122:   /*
123:      We are all finished solving linear systems, so we clean up the
124:      data structures.
125:   */
126:   PetscFree(rho);
127:   PetscFree(solution);
128:   PetscFree(userx);
129:   PetscFree(userb);
130:   UserFinalizeLinearSolver(&userctx);
131:   PetscFinalize();
132:   return ierr;
133: }

135: /* ------------------------------------------------------------------------*/
136: PetscErrorCode UserInitializeLinearSolver(PetscInt m,PetscInt n,UserCtx *userctx)
137: {
139:   PetscInt       N;

141:   /*
142:      Here we assume use of a grid of size m x n, with all points on the
143:      interior of the domain, i.e., we do not include the points corresponding
144:      to homogeneous Dirichlet boundary conditions.  We assume that the domain
145:      is [0,1]x[0,1].
146:   */
147:   userctx->m   = m;
148:   userctx->n   = n;
149:   userctx->hx2 = (m+1)*(m+1);
150:   userctx->hy2 = (n+1)*(n+1);
151:   N            = m*n;

153:   /*
154:      Create the sparse matrix. Preallocate 5 nonzeros per row.
155:   */
156:   MatCreateSeqAIJ(PETSC_COMM_SELF,N,N,5,0,&userctx->A);

158:   /*
159:      Create vectors. Here we create vectors with no memory allocated.
160:      This way, we can use the data structures already in the program
161:      by using VecPlaceArray() subroutine at a later stage.
162:   */
163:   VecCreateSeqWithArray(PETSC_COMM_SELF,1,N,NULL,&userctx->b);
164:   VecDuplicate(userctx->b,&userctx->x);

166:   /*
167:      Create linear solver context. This will be used repeatedly for all
168:      the linear solves needed.
169:   */
170:   KSPCreate(PETSC_COMM_SELF,&userctx->ksp);

172:   return 0;
173: }

175: /*
176:    Solves -div (rho grad psi) = F using finite differences.
177:    rho is a 2-dimensional array of size m by n, stored in Fortran
178:    style by columns. userb is a standard one-dimensional array.
179: */
180: /* ------------------------------------------------------------------------*/
181: PetscErrorCode UserDoLinearSolver(PetscScalar *rho,UserCtx *userctx,PetscScalar *userb,PetscScalar *userx)
182: {
184:   PetscInt       i,j,Ii,J,m = userctx->m,n = userctx->n;
185:   Mat            A = userctx->A;
186:   PC             pc;
187:   PetscScalar    v,hx2 = userctx->hx2,hy2 = userctx->hy2;

189:   /*
190:      This is not the most efficient way of generating the matrix
191:      but let's not worry about it. We should have separate code for
192:      the four corners, each edge and then the interior. Then we won't
193:      have the slow if-tests inside the loop.

195:      Computes the operator
196:              -div rho grad
197:      on an m by n grid with zero Dirichlet boundary conditions. The rho
198:      is assumed to be given on the same grid as the finite difference
199:      stencil is applied.  For a staggered grid, one would have to change
200:      things slightly.
201:   */
202:   Ii = 0;
203:   for (j=0; j<n; j++) {
204:     for (i=0; i<m; i++) {
205:       if (j>0) {
206:         J    = Ii - m;
207:         v    = -.5*(rho[Ii] + rho[J])*hy2;
208:         MatSetValues(A,1,&Ii,1,&J,&v,INSERT_VALUES);
209:       }
210:       if (j<n-1) {
211:         J    = Ii + m;
212:         v    = -.5*(rho[Ii] + rho[J])*hy2;
213:         MatSetValues(A,1,&Ii,1,&J,&v,INSERT_VALUES);
214:       }
215:       if (i>0) {
216:         J    = Ii - 1;
217:         v    = -.5*(rho[Ii] + rho[J])*hx2;
218:         MatSetValues(A,1,&Ii,1,&J,&v,INSERT_VALUES);
219:       }
220:       if (i<m-1) {
221:         J    = Ii + 1;
222:         v    = -.5*(rho[Ii] + rho[J])*hx2;
223:         MatSetValues(A,1,&Ii,1,&J,&v,INSERT_VALUES);
224:       }
225:       v    = 2.0*rho[Ii]*(hx2+hy2);
226:       MatSetValues(A,1,&Ii,1,&Ii,&v,INSERT_VALUES);
227:       Ii++;
228:     }
229:   }

231:   /*
232:      Assemble matrix
233:   */
234:   MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
235:   MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);

237:   /*
238:      Set operators. Here the matrix that defines the linear system
239:      also serves as the preconditioning matrix. Since all the matrices
240:      will have the same nonzero pattern here, we indicate this so the
241:      linear solvers can take advantage of this.
242:   */
243:   KSPSetOperators(userctx->ksp,A,A);

245:   /*
246:      Set linear solver defaults for this problem (optional).
247:      - Here we set it to use direct LU factorization for the solution
248:   */
249:   KSPGetPC(userctx->ksp,&pc);
250:   PCSetType(pc,PCLU);

252:   /*
253:      Set runtime options, e.g.,
254:         -ksp_type <type> -pc_type <type> -ksp_monitor -ksp_rtol <rtol>
255:      These options will override those specified above as long as
256:      KSPSetFromOptions() is called _after_ any other customization
257:      routines.

259:      Run the program with the option -help to see all the possible
260:      linear solver options.
261:   */
262:   KSPSetFromOptions(userctx->ksp);

264:   /*
265:      This allows the PETSc linear solvers to compute the solution
266:      directly in the user's array rather than in the PETSc vector.

268:      This is essentially a hack and not highly recommend unless you
269:      are quite comfortable with using PETSc. In general, users should
270:      write their entire application using PETSc vectors rather than
271:      arrays.
272:   */
273:   VecPlaceArray(userctx->x,userx);
274:   VecPlaceArray(userctx->b,userb);

276:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
277:                       Solve the linear system
278:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

280:   KSPSolve(userctx->ksp,userctx->b,userctx->x);

282:   /*
283:     Put back the PETSc array that belongs in the vector xuserctx->x
284:   */
285:   VecResetArray(userctx->x);
286:   VecResetArray(userctx->b);

288:   return 0;
289: }

291: /* ------------------------------------------------------------------------*/
292: PetscErrorCode UserFinalizeLinearSolver(UserCtx *userctx)
293: {
295:   /*
296:      We are all done and don't need to solve any more linear systems, so
297:      we free the work space.  All PETSc objects should be destroyed when
298:      they are no longer needed.
299:   */
300:   KSPDestroy(&userctx->ksp);
301:   VecDestroy(&userctx->x);
302:   VecDestroy(&userctx->b);
303:   MatDestroy(&userctx->A);
304:   return 0;
305: }