download the original source code.
1 /*
2 Example 12
3
4 Interface: Semi-Structured interface (SStruct)
5
6 Compile with: make ex12 (may need to edit HYPRE_DIR in Makefile)
7
8 Sample runs: mpirun -np 2 ex12 -pfmg
9 mpirun -np 2 ex12 -boomeramg
10
11 Description: The grid layout is the same as ex1, but with nodal unknowns. The
12 solver is PCG preconditioned with either PFMG or BoomerAMG,
13 selected on the command line.
14
15 We recommend viewing the Struct examples before viewing this
16 and the other SStruct examples. This is one of the simplest
17 SStruct examples, used primarily to demonstrate how to set up
18 non-cell-centered problems, and to demonstrate how easy it is
19 to switch between structured solvers (PFMG) and solvers
20 designed for more general settings (AMG).
21 */
22
23 #include <stdio.h>
24 #include <stdlib.h>
25 #include <string.h>
26
27 #include "HYPRE_sstruct_ls.h"
28 #include "HYPRE_parcsr_ls.h"
29 #include "HYPRE_krylov.h"
30
31 #include "vis.c"
32
33 int main (int argc, char *argv[])
34 {
35 int i, j, myid, num_procs;
36
37 int vis = 0;
38
39 HYPRE_SStructGrid grid;
40 HYPRE_SStructGraph graph;
41 HYPRE_SStructStencil stencil;
42 HYPRE_SStructMatrix A;
43 HYPRE_SStructVector b;
44 HYPRE_SStructVector x;
45
46 /* We only have one part and one variable */
47 int nparts = 1;
48 int nvars = 1;
49 int part = 0;
50 int var = 0;
51
52 int precond_id = 1;
53 int object_type = HYPRE_STRUCT;
54
55 /* Initialize MPI */
56 MPI_Init(&argc, &argv);
57 MPI_Comm_rank(MPI_COMM_WORLD, &myid);
58 MPI_Comm_size(MPI_COMM_WORLD, &num_procs);
59
60 if (num_procs != 2)
61 {
62 if (myid == 0) printf("Must run with 2 processors!\n");
63 exit(1);
64 }
65
66 /* Parse command line */
67 {
68 int arg_index = 0;
69 int print_usage = 0;
70
71 while (arg_index < argc)
72 {
73 if ( strcmp(argv[arg_index], "-pfmg") == 0 )
74 {
75 arg_index++;
76 precond_id = 1;
77 object_type = HYPRE_STRUCT;
78 }
79 else if ( strcmp(argv[arg_index], "-boomeramg") == 0 )
80 {
81 arg_index++;
82 precond_id = 2;
83 object_type = HYPRE_PARCSR;
84 }
85 else if ( strcmp(argv[arg_index], "-vis") == 0 )
86 {
87 arg_index++;
88 vis = 1;
89 }
90 else if ( strcmp(argv[arg_index], "-help") == 0 )
91 {
92 print_usage = 1;
93 break;
94 }
95 else
96 {
97 arg_index++;
98 }
99 }
100
101 if ((print_usage) && (myid == 0))
102 {
103 printf("\n");
104 printf("Usage: %s [<options>]\n", argv[0]);
105 printf("\n");
106 printf(" -pfmg : use the structured PFMG solver (default)\n");
107 printf(" -boomeramg : use the unstructured BoomerAMG solver\n");
108 printf(" -vis : save the solution for GLVis visualization\n");
109 printf("\n");
110 }
111
112 if (print_usage)
113 {
114 MPI_Finalize();
115 return (0);
116 }
117 }
118
119 /* 1. Set up the grid. Here we use only one part. Each processor describes
120 the piece of the grid that it owns. */
121 {
122 /* Create an empty 2D grid object */
123 HYPRE_SStructGridCreate(MPI_COMM_WORLD, 2, nparts, &grid);
124
125 /* Add boxes to the grid */
126 if (myid == 0)
127 {
128 int ilower[2]={-3,1}, iupper[2]={-1,2};
129 HYPRE_SStructGridSetExtents(grid, part, ilower, iupper);
130 }
131 else if (myid == 1)
132 {
133 int ilower[2]={0,1}, iupper[2]={2,4};
134 HYPRE_SStructGridSetExtents(grid, part, ilower, iupper);
135 }
136
137 /* Set the variable type and number of variables on each part. */
138 {
139 HYPRE_SStructVariable vartypes[1] = {HYPRE_SSTRUCT_VARIABLE_NODE};
140
141 HYPRE_SStructGridSetVariables(grid, part, nvars, vartypes);
142 }
143
144 /* This is a collective call finalizing the grid assembly.
145 The grid is now ``ready to be used'' */
146 HYPRE_SStructGridAssemble(grid);
147 }
148
149 /* 2. Define the discretization stencil */
150 {
151 /* Create an empty 2D, 5-pt stencil object */
152 HYPRE_SStructStencilCreate(2, 5, &stencil);
153
154 /* Define the geometry of the stencil. Each represents a relative offset
155 (in the index space). */
156 {
157 int entry;
158 int offsets[5][2] = {{0,0}, {-1,0}, {1,0}, {0,-1}, {0,1}};
159
160 /* Assign numerical values to the offsets so that we can easily refer
161 to them - the last argument indicates the variable for which we are
162 assigning this stencil */
163 for (entry = 0; entry < 5; entry++)
164 HYPRE_SStructStencilSetEntry(stencil, entry, offsets[entry], var);
165 }
166 }
167
168 /* 3. Set up the Graph - this determines the non-zero structure of the matrix
169 and allows non-stencil relationships between the parts */
170 {
171 /* Create the graph object */
172 HYPRE_SStructGraphCreate(MPI_COMM_WORLD, grid, &graph);
173
174 /* See MatrixSetObjectType below */
175 HYPRE_SStructGraphSetObjectType(graph, object_type);
176
177 /* Now we need to tell the graph which stencil to use for each variable on
178 each part (we only have one variable and one part) */
179 HYPRE_SStructGraphSetStencil(graph, part, var, stencil);
180
181 /* Here we could establish connections between parts if we had more than
182 one part using the graph. For example, we could use
183 HYPRE_GraphAddEntries() routine or HYPRE_GridSetNeighborPart() */
184
185 /* Assemble the graph */
186 HYPRE_SStructGraphAssemble(graph);
187 }
188
189 /* 4. Set up a SStruct Matrix */
190 {
191 /* Create an empty matrix object */
192 HYPRE_SStructMatrixCreate(MPI_COMM_WORLD, graph, &A);
193
194 /* Set the object type (by default HYPRE_SSTRUCT). This determines the
195 data structure used to store the matrix. For PFMG we need to use
196 HYPRE_STRUCT, and for BoomerAMG we need HYPRE_PARCSR (set above). */
197 HYPRE_SStructMatrixSetObjectType(A, object_type);
198
199 /* Get ready to set values */
200 HYPRE_SStructMatrixInitialize(A);
201
202 /* Set the matrix coefficients. Each processor assigns coefficients for
203 the boxes in the grid that it owns. Note that the coefficients
204 associated with each stencil entry may vary from grid point to grid
205 point if desired. Here, we first set the same stencil entries for each
206 grid point. Then we make modifications to grid points near the
207 boundary. Note that the ilower values are different from those used in
208 ex1 because of the way nodal variables are referenced. Also note that
209 some of the stencil values are set on both processor 0 and processor 1.
210 See the User and Reference manuals for more details. */
211 if (myid == 0)
212 {
213 int ilower[2]={-4,0}, iupper[2]={-1,2};
214 int stencil_indices[5] = {0,1,2,3,4}; /* labels for the stencil entries -
215 these correspond to the offsets
216 defined above */
217 int nentries = 5;
218 int nvalues = 60; /* 12 grid points, each with 5 stencil entries */
219 double values[60];
220
221 for (i = 0; i < nvalues; i += nentries)
222 {
223 values[i] = 4.0;
224 for (j = 1; j < nentries; j++)
225 values[i+j] = -1.0;
226 }
227
228 HYPRE_SStructMatrixSetBoxValues(A, part, ilower, iupper, var, nentries,
229 stencil_indices, values);
230 }
231 else if (myid == 1)
232 {
233 int ilower[2]={-1,0}, iupper[2]={2,4};
234 int stencil_indices[5] = {0,1,2,3,4};
235 int nentries = 5;
236 int nvalues = 100; /* 20 grid points, each with 5 stencil entries */
237 double values[100];
238
239 for (i = 0; i < nvalues; i += nentries)
240 {
241 values[i] = 4.0;
242 for (j = 1; j < nentries; j++)
243 values[i+j] = -1.0;
244 }
245
246 HYPRE_SStructMatrixSetBoxValues(A, part, ilower, iupper, var, nentries,
247 stencil_indices, values);
248 }
249
250 /* Set the coefficients reaching outside of the boundary to 0. Note that
251 * both ilower *and* iupper may be different from those in ex1. */
252 if (myid == 0)
253 {
254 double values[4];
255 for (i = 0; i < 4; i++)
256 values[i] = 0.0;
257 {
258 /* values below our box */
259 int ilower[2]={-4,0}, iupper[2]={-1,0};
260 int stencil_indices[1] = {3};
261 HYPRE_SStructMatrixSetBoxValues(A, part, ilower, iupper, var, 1,
262 stencil_indices, values);
263 }
264 {
265 /* values to the left of our box */
266 int ilower[2]={-4,0}, iupper[2]={-4,2};
267 int stencil_indices[1] = {1};
268 HYPRE_SStructMatrixSetBoxValues(A, part, ilower, iupper, var, 1,
269 stencil_indices, values);
270 }
271 {
272 /* values above our box */
273 int ilower[2]={-4,2}, iupper[2]={-2,2};
274 int stencil_indices[1] = {4};
275 HYPRE_SStructMatrixSetBoxValues(A, part, ilower, iupper, var, 1,
276 stencil_indices, values);
277 }
278 }
279 else if (myid == 1)
280 {
281 double values[5];
282 for (i = 0; i < 5; i++)
283 values[i] = 0.0;
284 {
285 /* values below our box */
286 int ilower[2]={-1,0}, iupper[2]={2,0};
287 int stencil_indices[1] = {3};
288 HYPRE_SStructMatrixSetBoxValues(A, part, ilower, iupper, var, 1,
289 stencil_indices, values);
290 }
291 {
292 /* values to the right of our box */
293 int ilower[2]={2,0}, iupper[2]={2,4};
294 int stencil_indices[1] = {2};
295 HYPRE_SStructMatrixSetBoxValues(A, part, ilower, iupper, var, 1,
296 stencil_indices, values);
297 }
298 {
299 /* values above our box */
300 int ilower[2]={-1,4}, iupper[2]={2,4};
301 int stencil_indices[1] = {4};
302 HYPRE_SStructMatrixSetBoxValues(A, part, ilower, iupper, var, 1,
303 stencil_indices, values);
304 }
305 {
306 /* values to the left of our box
307 (that do not border the other box on proc. 0) */
308 int ilower[2]={-1,3}, iupper[2]={-1,4};
309 int stencil_indices[1] = {1};
310 HYPRE_SStructMatrixSetBoxValues(A, part, ilower, iupper, var, 1,
311 stencil_indices, values);
312 }
313 }
314
315 /* This is a collective call finalizing the matrix assembly.
316 The matrix is now ``ready to be used'' */
317 HYPRE_SStructMatrixAssemble(A);
318 }
319
320 /* 5. Set up SStruct Vectors for b and x. */
321 {
322 /* Create an empty vector object */
323 HYPRE_SStructVectorCreate(MPI_COMM_WORLD, grid, &b);
324 HYPRE_SStructVectorCreate(MPI_COMM_WORLD, grid, &x);
325
326 /* As with the matrix, set the appropriate object type for the vectors */
327 HYPRE_SStructVectorSetObjectType(b, object_type);
328 HYPRE_SStructVectorSetObjectType(x, object_type);
329
330 /* Indicate that the vector coefficients are ready to be set */
331 HYPRE_SStructVectorInitialize(b);
332 HYPRE_SStructVectorInitialize(x);
333
334 /* Set the vector coefficients. Again, note that the ilower values are
335 different from those used in ex1, and some of the values are set on
336 both processors. */
337 if (myid == 0)
338 {
339 int ilower[2]={-4,0}, iupper[2]={-1,2};
340 double values[12]; /* 12 grid points */
341
342 for (i = 0; i < 12; i ++)
343 values[i] = 1.0;
344 HYPRE_SStructVectorSetBoxValues(b, part, ilower, iupper, var, values);
345
346 for (i = 0; i < 12; i ++)
347 values[i] = 0.0;
348 HYPRE_SStructVectorSetBoxValues(x, part, ilower, iupper, var, values);
349 }
350 else if (myid == 1)
351 {
352 int ilower[2]={0,1}, iupper[2]={2,4};
353 double values[20]; /* 20 grid points */
354
355 for (i = 0; i < 20; i ++)
356 values[i] = 1.0;
357 HYPRE_SStructVectorSetBoxValues(b, part, ilower, iupper, var, values);
358
359 for (i = 0; i < 20; i ++)
360 values[i] = 0.0;
361 HYPRE_SStructVectorSetBoxValues(x, part, ilower, iupper, var, values);
362 }
363
364 /* This is a collective call finalizing the vector assembly.
365 The vectors are now ``ready to be used'' */
366 HYPRE_SStructVectorAssemble(b);
367 HYPRE_SStructVectorAssemble(x);
368 }
369
370 /* 6. Set up and use a solver (See the Reference Manual for descriptions
371 of all of the options.) */
372 if (precond_id == 1) /* PFMG */
373 {
374 HYPRE_StructMatrix sA;
375 HYPRE_StructVector sb;
376 HYPRE_StructVector sx;
377
378 HYPRE_StructSolver solver;
379 HYPRE_StructSolver precond;
380
381 /* Because we are using a struct solver, we need to get the
382 object of the matrix and vectors to pass in to the struct solvers */
383 HYPRE_SStructMatrixGetObject(A, (void **) &sA);
384 HYPRE_SStructVectorGetObject(b, (void **) &sb);
385 HYPRE_SStructVectorGetObject(x, (void **) &sx);
386
387 /* Create an empty PCG Struct solver */
388 HYPRE_StructPCGCreate(MPI_COMM_WORLD, &solver);
389
390 /* Set PCG parameters */
391 HYPRE_StructPCGSetTol(solver, 1.0e-06);
392 HYPRE_StructPCGSetPrintLevel(solver, 2);
393 HYPRE_StructPCGSetMaxIter(solver, 50);
394
395 /* Create the Struct PFMG solver for use as a preconditioner */
396 HYPRE_StructPFMGCreate(MPI_COMM_WORLD, &precond);
397
398 /* Set PFMG parameters */
399 HYPRE_StructPFMGSetMaxIter(precond, 1);
400 HYPRE_StructPFMGSetTol(precond, 0.0);
401 HYPRE_StructPFMGSetZeroGuess(precond);
402 HYPRE_StructPFMGSetNumPreRelax(precond, 2);
403 HYPRE_StructPFMGSetNumPostRelax(precond, 2);
404 /* non-Galerkin coarse grid (more efficient for this problem) */
405 HYPRE_StructPFMGSetRAPType(precond, 1);
406 /* R/B Gauss-Seidel */
407 HYPRE_StructPFMGSetRelaxType(precond, 2);
408 /* skip relaxation on some levels (more efficient for this problem) */
409 HYPRE_StructPFMGSetSkipRelax(precond, 1);
410
411
412 /* Set preconditioner and solve */
413 HYPRE_StructPCGSetPrecond(solver, HYPRE_StructPFMGSolve,
414 HYPRE_StructPFMGSetup, precond);
415 HYPRE_StructPCGSetup(solver, sA, sb, sx);
416 HYPRE_StructPCGSolve(solver, sA, sb, sx);
417
418 /* Free memory */
419 HYPRE_StructPCGDestroy(solver);
420 HYPRE_StructPFMGDestroy(precond);
421 }
422 else if (precond_id == 2) /* BoomerAMG */
423 {
424 HYPRE_ParCSRMatrix parA;
425 HYPRE_ParVector parb;
426 HYPRE_ParVector parx;
427
428 HYPRE_Solver solver;
429 HYPRE_Solver precond;
430
431 /* Because we are using a struct solver, we need to get the
432 object of the matrix and vectors to pass in to the struct solvers */
433 HYPRE_SStructMatrixGetObject(A, (void **) &parA);
434 HYPRE_SStructVectorGetObject(b, (void **) &parb);
435 HYPRE_SStructVectorGetObject(x, (void **) &parx);
436
437 /* Create an empty PCG Struct solver */
438 HYPRE_ParCSRPCGCreate(MPI_COMM_WORLD, &solver);
439
440 /* Set PCG parameters */
441 HYPRE_ParCSRPCGSetTol(solver, 1.0e-06);
442 HYPRE_ParCSRPCGSetPrintLevel(solver, 2);
443 HYPRE_ParCSRPCGSetMaxIter(solver, 50);
444
445 /* Create the BoomerAMG solver for use as a preconditioner */
446 HYPRE_BoomerAMGCreate(&precond);
447
448 /* Set BoomerAMG parameters */
449 HYPRE_BoomerAMGSetMaxIter(precond, 1);
450 HYPRE_BoomerAMGSetTol(precond, 0.0);
451 HYPRE_BoomerAMGSetPrintLevel(precond, 1); /* print amg solution info */
452 HYPRE_BoomerAMGSetRelaxType(precond, 6); /* Sym G.S./Jacobi hybrid */
453 HYPRE_BoomerAMGSetNumSweeps(precond, 1);
454
455 /* Set preconditioner and solve */
456 HYPRE_ParCSRPCGSetPrecond(solver, HYPRE_BoomerAMGSolve,
457 HYPRE_BoomerAMGSetup, precond);
458 HYPRE_ParCSRPCGSetup(solver, parA, parb, parx);
459 HYPRE_ParCSRPCGSolve(solver, parA, parb, parx);
460
461 /* Free memory */
462 HYPRE_ParCSRPCGDestroy(solver);
463 HYPRE_BoomerAMGDestroy(precond);
464 }
465
466 /* Save the solution for GLVis visualization, see vis/glvis-ex12.sh */
467 if (vis)
468 {
469 /* Gather the solution vector */
470 HYPRE_SStructVectorGather(x);
471
472 GLVis_PrintSStructGrid(grid, "vis/ex12.mesh", myid, NULL, NULL);
473 GLVis_PrintSStructVector(x, 0, "vis/ex12.sol", myid);
474 GLVis_PrintData("vis/ex12.data", myid, num_procs);
475 }
476
477 /* Free memory */
478 HYPRE_SStructGridDestroy(grid);
479 HYPRE_SStructStencilDestroy(stencil);
480 HYPRE_SStructGraphDestroy(graph);
481 HYPRE_SStructMatrixDestroy(A);
482 HYPRE_SStructVectorDestroy(b);
483 HYPRE_SStructVectorDestroy(x);
484
485 /* Finalize MPI */
486 MPI_Finalize();
487
488 return (0);
489 }
syntax highlighted by Code2HTML, v. 0.9.1