download the original source code.
1 /*
2 Example 1
3
4 Interface: Structured interface (Struct)
5
6 Compile with: make ex1 (may need to edit HYPRE_DIR in Makefile)
7
8 Sample run: mpirun -np 2 ex1
9
10 Description: This is a two processor example. Each processor owns one
11 box in the grid. For reference, the two grid boxes are those
12 in the example diagram in the struct interface chapter
13 of the User's Manual. Note that in this example code, we have
14 used the two boxes shown in the diagram as belonging
15 to processor 0 (and given one box to each processor). The
16 solver is PCG with no preconditioner.
17
18 We recommend viewing examples 1-4 sequentially for
19 a nice overview/tutorial of the struct interface.
20 */
21
22 #include <stdio.h>
23
24 /* Struct linear solvers header */
25 #include "HYPRE_struct_ls.h"
26
27 #include "vis.c"
28
29 int main (int argc, char *argv[])
30 {
31 int i, j, myid, num_procs;
32
33 int vis = 0;
34
35 HYPRE_StructGrid grid;
36 HYPRE_StructStencil stencil;
37 HYPRE_StructMatrix A;
38 HYPRE_StructVector b;
39 HYPRE_StructVector x;
40 HYPRE_StructSolver solver;
41
42 /* Initialize MPI */
43 MPI_Init(&argc, &argv);
44 MPI_Comm_rank(MPI_COMM_WORLD, &myid);
45 MPI_Comm_size(MPI_COMM_WORLD, &num_procs);
46
47 if (num_procs != 2)
48 {
49 if (myid == 0) printf("Must run with 2 processors!\n");
50 MPI_Finalize();
51
52 return(0);
53 }
54
55 /* Parse command line */
56 {
57 int arg_index = 0;
58 int print_usage = 0;
59
60 while (arg_index < argc)
61 {
62 if ( strcmp(argv[arg_index], "-vis") == 0 )
63 {
64 arg_index++;
65 vis = 1;
66 }
67 else if ( strcmp(argv[arg_index], "-help") == 0 )
68 {
69 print_usage = 1;
70 break;
71 }
72 else
73 {
74 arg_index++;
75 }
76 }
77
78 if ((print_usage) && (myid == 0))
79 {
80 printf("\n");
81 printf("Usage: %s [<options>]\n", argv[0]);
82 printf("\n");
83 printf(" -vis : save the solution for GLVis visualization\n");
84 printf("\n");
85 }
86
87 if (print_usage)
88 {
89 MPI_Finalize();
90 return (0);
91 }
92 }
93
94 /* 1. Set up a grid. Each processor describes the piece
95 of the grid that it owns. */
96 {
97 /* Create an empty 2D grid object */
98 HYPRE_StructGridCreate(MPI_COMM_WORLD, 2, &grid);
99
100 /* Add boxes to the grid */
101 if (myid == 0)
102 {
103 int ilower[2]={-3,1}, iupper[2]={-1,2};
104 HYPRE_StructGridSetExtents(grid, ilower, iupper);
105 }
106 else if (myid == 1)
107 {
108 int ilower[2]={0,1}, iupper[2]={2,4};
109 HYPRE_StructGridSetExtents(grid, ilower, iupper);
110 }
111
112 /* This is a collective call finalizing the grid assembly.
113 The grid is now ``ready to be used'' */
114 HYPRE_StructGridAssemble(grid);
115 }
116
117 /* 2. Define the discretization stencil */
118 {
119 /* Create an empty 2D, 5-pt stencil object */
120 HYPRE_StructStencilCreate(2, 5, &stencil);
121
122 /* Define the geometry of the stencil. Each represents a
123 relative offset (in the index space). */
124 {
125 int entry;
126 int offsets[5][2] = {{0,0}, {-1,0}, {1,0}, {0,-1}, {0,1}};
127
128 /* Assign each of the 5 stencil entries */
129 for (entry = 0; entry < 5; entry++)
130 HYPRE_StructStencilSetElement(stencil, entry, offsets[entry]);
131 }
132 }
133
134 /* 3. Set up a Struct Matrix */
135 {
136 /* Create an empty matrix object */
137 HYPRE_StructMatrixCreate(MPI_COMM_WORLD, grid, stencil, &A);
138
139 /* Indicate that the matrix coefficients are ready to be set */
140 HYPRE_StructMatrixInitialize(A);
141
142 /* Set the matrix coefficients. Each processor assigns coefficients
143 for the boxes in the grid that it owns. Note that the coefficients
144 associated with each stencil entry may vary from grid point to grid
145 point if desired. Here, we first set the same stencil entries for
146 each grid point. Then we make modifications to grid points near
147 the boundary. */
148 if (myid == 0)
149 {
150 int ilower[2]={-3,1}, iupper[2]={-1,2};
151 int stencil_indices[5] = {0,1,2,3,4}; /* labels for the stencil entries -
152 these correspond to the offsets
153 defined above */
154 int nentries = 5;
155 int nvalues = 30; /* 6 grid points, each with 5 stencil entries */
156 double values[30];
157
158 /* We have 6 grid points, each with 5 stencil entries */
159 for (i = 0; i < nvalues; i += nentries)
160 {
161 values[i] = 4.0;
162 for (j = 1; j < nentries; j++)
163 values[i+j] = -1.0;
164 }
165
166 HYPRE_StructMatrixSetBoxValues(A, ilower, iupper, nentries,
167 stencil_indices, values);
168 }
169 else if (myid == 1)
170 {
171 int ilower[2]={0,1}, iupper[2]={2,4};
172 int stencil_indices[5] = {0,1,2,3,4};
173 int nentries = 5;
174 int nvalues = 60; /* 12 grid points, each with 5 stencil entries */
175 double values[60];
176
177 for (i = 0; i < nvalues; i += nentries)
178 {
179 values[i] = 4.0;
180 for (j = 1; j < nentries; j++)
181 values[i+j] = -1.0;
182 }
183
184 HYPRE_StructMatrixSetBoxValues(A, ilower, iupper, nentries,
185 stencil_indices, values);
186 }
187
188 /* Set the coefficients reaching outside of the boundary to 0 */
189 if (myid == 0)
190 {
191 double values[3];
192 for (i = 0; i < 3; i++)
193 values[i] = 0.0;
194 {
195 /* values below our box */
196 int ilower[2]={-3,1}, iupper[2]={-1,1};
197 int stencil_indices[1] = {3};
198 HYPRE_StructMatrixSetBoxValues(A, ilower, iupper, 1,
199 stencil_indices, values);
200 }
201 {
202 /* values to the left of our box */
203 int ilower[2]={-3,1}, iupper[2]={-3,2};
204 int stencil_indices[1] = {1};
205 HYPRE_StructMatrixSetBoxValues(A, ilower, iupper, 1,
206 stencil_indices, values);
207 }
208 {
209 /* values above our box */
210 int ilower[2]={-3,2}, iupper[2]={-1,2};
211 int stencil_indices[1] = {4};
212 HYPRE_StructMatrixSetBoxValues(A, ilower, iupper, 1,
213 stencil_indices, values);
214 }
215 }
216 else if (myid == 1)
217 {
218 double values[4];
219 for (i = 0; i < 4; i++)
220 values[i] = 0.0;
221 {
222 /* values below our box */
223 int ilower[2]={0,1}, iupper[2]={2,1};
224 int stencil_indices[1] = {3};
225 HYPRE_StructMatrixSetBoxValues(A, ilower, iupper, 1,
226 stencil_indices, values);
227 }
228 {
229 /* values to the right of our box */
230 int ilower[2]={2,1}, iupper[2]={2,4};
231 int stencil_indices[1] = {2};
232 HYPRE_StructMatrixSetBoxValues(A, ilower, iupper, 1,
233 stencil_indices, values);
234 }
235 {
236 /* values above our box */
237 int ilower[2]={0,4}, iupper[2]={2,4};
238 int stencil_indices[1] = {4};
239 HYPRE_StructMatrixSetBoxValues(A, ilower, iupper, 1,
240 stencil_indices, values);
241 }
242 {
243 /* values to the left of our box
244 (that do not border the other box on proc. 0) */
245 int ilower[2]={0,3}, iupper[2]={0,4};
246 int stencil_indices[1] = {1};
247 HYPRE_StructMatrixSetBoxValues(A, ilower, iupper, 1,
248 stencil_indices, values);
249 }
250 }
251
252 /* This is a collective call finalizing the matrix assembly.
253 The matrix is now ``ready to be used'' */
254 HYPRE_StructMatrixAssemble(A);
255 }
256
257 /* 4. Set up Struct Vectors for b and x. Each processor sets the vectors
258 corresponding to its boxes. */
259 {
260 /* Create an empty vector object */
261 HYPRE_StructVectorCreate(MPI_COMM_WORLD, grid, &b);
262 HYPRE_StructVectorCreate(MPI_COMM_WORLD, grid, &x);
263
264 /* Indicate that the vector coefficients are ready to be set */
265 HYPRE_StructVectorInitialize(b);
266 HYPRE_StructVectorInitialize(x);
267
268 /* Set the vector coefficients */
269 if (myid == 0)
270 {
271 int ilower[2]={-3,1}, iupper[2]={-1,2};
272 double values[6]; /* 6 grid points */
273
274 for (i = 0; i < 6; i ++)
275 values[i] = 1.0;
276 HYPRE_StructVectorSetBoxValues(b, ilower, iupper, values);
277
278 for (i = 0; i < 6; i ++)
279 values[i] = 0.0;
280 HYPRE_StructVectorSetBoxValues(x, ilower, iupper, values);
281 }
282 else if (myid == 1)
283 {
284 int ilower[2]={0,1}, iupper[2]={2,4};
285 double values[12]; /* 12 grid points */
286
287 for (i = 0; i < 12; i ++)
288 values[i] = 1.0;
289 HYPRE_StructVectorSetBoxValues(b, ilower, iupper, values);
290
291 for (i = 0; i < 12; i ++)
292 values[i] = 0.0;
293 HYPRE_StructVectorSetBoxValues(x, ilower, iupper, values);
294 }
295
296 /* This is a collective call finalizing the vector assembly.
297 The vectors are now ``ready to be used'' */
298 HYPRE_StructVectorAssemble(b);
299 HYPRE_StructVectorAssemble(x);
300 }
301
302 /* 5. Set up and use a solver (See the Reference Manual for descriptions
303 of all of the options.) */
304 {
305 /* Create an empty PCG Struct solver */
306 HYPRE_StructPCGCreate(MPI_COMM_WORLD, &solver);
307
308 /* Set some parameters */
309 HYPRE_StructPCGSetTol(solver, 1.0e-06); /* convergence tolerance */
310 HYPRE_StructPCGSetPrintLevel(solver, 2); /* amount of info. printed */
311
312 /* Setup and solve */
313 HYPRE_StructPCGSetup(solver, A, b, x);
314 HYPRE_StructPCGSolve(solver, A, b, x);
315 }
316
317 /* Save the solution for GLVis visualization, see vis/glvis-ex1.sh */
318 if (vis)
319 {
320 GLVis_PrintStructGrid(grid, "vis/ex1.mesh", myid, NULL, NULL);
321 GLVis_PrintStructVector(x, "vis/ex1.sol", myid);
322 GLVis_PrintData("vis/ex1.data", myid, num_procs);
323 }
324
325 /* Free memory */
326 HYPRE_StructGridDestroy(grid);
327 HYPRE_StructStencilDestroy(stencil);
328 HYPRE_StructMatrixDestroy(A);
329 HYPRE_StructVectorDestroy(b);
330 HYPRE_StructVectorDestroy(x);
331 HYPRE_StructPCGDestroy(solver);
332
333 /* Finalize MPI */
334 MPI_Finalize();
335
336 return (0);
337 }
syntax highlighted by Code2HTML, v. 0.9.1