source: CIVL/examples/mpi-omp/AMG2013/parcsr_ls/par_vardifconv.c

main
Last change on this file was ea777aa, checked in by Alex Wilton <awilton@…>, 3 years ago

Moved examples, include, build_default.properties, common.xml, and README out from dev.civl.com into the root of the repo.

git-svn-id: ​svn://vsl.cis.udel.edu/civl/trunk@5704 fb995dde-84ed-4084-dfe6-e5aef3e2452c
  • Property mode set to 100644
File size: 16.1 KB
RevLine 
[2aa6644]1/*BHEADER**********************************************************************
2 * Copyright (c) 2008, Lawrence Livermore National Security, LLC.
3 * Produced at the Lawrence Livermore National Laboratory.
4 * This file is part of HYPRE. See file COPYRIGHT for details.
5 *
6 * HYPRE is free software; you can redistribute it and/or modify it under the
7 * terms of the GNU Lesser General Public License (as published by the Free
8 * Software Foundation) version 2.1 dated February 1999.
9 *
10 * $Revision: 2.4 $
11 ***********************************************************************EHEADER*/
12
13
14
15
16#include "headers.h"
17
18/*--------------------------------------------------------------------------
19 * hypre_GenerateVarDifConv
20 *--------------------------------------------------------------------------*/
21
22HYPRE_ParCSRMatrix
23GenerateVarDifConv( MPI_Comm comm,
24 HYPRE_BigInt nx,
25 HYPRE_BigInt ny,
26 HYPRE_BigInt nz,
27 int P,
28 int Q,
29 int R,
30 int p,
31 int q,
32 int r,
33 double eps,
34 HYPRE_ParVector *rhs_ptr,
35 HYPRE_ParVector *x_ptr)
36{
37 hypre_ParCSRMatrix *A;
38 hypre_CSRMatrix *diag;
39 hypre_CSRMatrix *offd;
40 hypre_ParVector *par_rhs;
41 hypre_Vector *rhs;
42 hypre_ParVector *par_x;
43 hypre_Vector *x;
44 double *rhs_data;
45 double *x_data;
46
47 int *diag_i;
48 int *diag_j;
49 double *diag_data;
50
51 int *offd_i;
52 int *offd_j = NULL;
53 double *offd_data = NULL;
54
55 HYPRE_BigInt *global_part;
56 HYPRE_BigInt *global_part_rhs;
57 HYPRE_BigInt *global_part_x;
58 int *tmp_j = NULL;
59 int ix, iy, iz;
60 int cnt, o_cnt;
61 int local_num_rows;
62 HYPRE_BigInt *col_map_offd = NULL;
63 int row_index;
64 int i, j;
65 int ip, iq, ir;
66
67 int nx_local, ny_local, nz_local;
68 int num_cols_offd;
69 HYPRE_BigInt grid_size;
70
71 HYPRE_BigInt *nx_part;
72 HYPRE_BigInt *ny_part;
73 HYPRE_BigInt *nz_part;
74
75 int num_procs, my_id;
76 int P_busy, Q_busy, R_busy;
77
78 double hhx, hhy, hhz;
79 double xx, yy, zz;
80 double afp, afm, bfp, bfm, cfp, cfm, df, ef, ff, gf;
81
82#ifdef HYPRE_NO_GLOBAL_PARTITION
83 MPI_Comm_size(comm,&num_procs);
84 MPI_Comm_rank(comm,&my_id);
85
86 hypre_GenerateLocalPartitioning(nx,P,p,&nx_part);
87 hypre_GenerateLocalPartitioning(ny,Q,q,&ny_part);
88 hypre_GenerateLocalPartitioning(nz,R,r,&nz_part);
89
90 nx_local = (int)(nx_part[1] - nx_part[0]);
91 ny_local = (int)(ny_part[1] - ny_part[0]);
92 nz_local = (int)(nz_part[1] - nz_part[0]);
93
94 local_num_rows = nx_local*ny_local*nz_local;
95
96 global_part = hypre_CTAlloc(HYPRE_BigInt,2);
97 global_part_rhs = hypre_CTAlloc(HYPRE_BigInt,2);
98 global_part_x = hypre_CTAlloc(HYPRE_BigInt,2);
99
100 global_part[0] = (HYPRE_BigInt)my_id*(HYPRE_BigInt)local_num_rows;
101 global_part[1] = global_part[0]+(HYPRE_BigInt)local_num_rows;
102 global_part_rhs[0] = global_part[0];
103 global_part_rhs[1] = global_part[1];
104 global_part_x[0] = global_part[0];
105 global_part_x[1] = global_part[1];
106
107 ip = p;
108 iq = q;
109 ir = r;
110#else
111 int nx_size, ny_size, nz_size;
112 MPI_Comm_size(comm,&num_procs);
113 MPI_Comm_rank(comm,&my_id);
114
115 hypre_GeneratePartitioning(nx,P,&nx_part);
116 hypre_GeneratePartitioning(ny,Q,&ny_part);
117 hypre_GeneratePartitioning(nz,R,&nz_part);
118
119 global_part = hypre_CTAlloc(HYPRE_BigInt,P*Q*R+1);
120 global_part_rhs = hypre_CTAlloc(HYPRE_BigInt,P*Q*R+1);
121 global_part_x = hypre_CTAlloc(HYPRE_BigInt,P*Q*R+1);
122
123 global_part[0] = 0;
124 global_part_rhs[0] = 0;
125 global_part_x[0] = 0;
126 cnt = 1;
127 for (iz = 0; iz < R; iz++)
128 {
129 nz_size = (int)(nz_part[iz+1]-nz_part[iz]);
130 for (iy = 0; iy < Q; iy++)
131 {
132 ny_size = (int)(ny_part[iy+1]-ny_part[iy]);
133 for (ix = 0; ix < P; ix++)
134 {
135 nx_size = (int)(nx_part[ix+1] - nx_part[ix]);
136 global_part[cnt] = global_part[cnt-1];
137 global_part[cnt] += (HYPRE_BigInt)(nx_size*ny_size*nz_size);
138 global_part_rhs[cnt] = global_part[cnt];
139 global_part_x[cnt] = global_part[cnt];
140 cnt++;
141 }
142 }
143 }
144
145 nx_local = (int)(nx_part[p+1] - nx_part[p]);
146 ny_local = (int)(ny_part[q+1] - ny_part[q]);
147 nz_local = (int)(nz_part[r+1] - nz_part[r]);
148
149 local_num_rows = nx_local*ny_local*nz_local;
150 ip = p;
151 iq = q;
152 ir = r;
153#endif
154
155 grid_size = nx*ny*nz;
156
157 diag_i = hypre_CTAlloc(int, local_num_rows+1);
158 offd_i = hypre_CTAlloc(int, local_num_rows+1);
159 rhs_data = hypre_CTAlloc(double, local_num_rows);
160 x_data = hypre_CTAlloc(double, local_num_rows);
161 for (i=0; i < local_num_rows; i++)
162 x_data[i] = 0.0;
163
164 P_busy = hypre_min(nx,P);
165 Q_busy = hypre_min(ny,Q);
166 R_busy = hypre_min(nz,R);
167
168 num_cols_offd = 0;
169 if (p) num_cols_offd += ny_local*nz_local;
170 if (p < P_busy-1) num_cols_offd += ny_local*nz_local;
171 if (q) num_cols_offd += nx_local*nz_local;
172 if (q < Q_busy-1) num_cols_offd += nx_local*nz_local;
173 if (r) num_cols_offd += nx_local*ny_local;
174 if (r < R_busy-1) num_cols_offd += nx_local*ny_local;
175
176 if (!local_num_rows) num_cols_offd = 0;
177
178 if (num_cols_offd)
179 col_map_offd = hypre_CTAlloc(HYPRE_BigInt, num_cols_offd);
180
181 hhx = 1.0/(double)(nx+1);
182 hhy = 1.0/(double)(ny+1);
183 hhz = 1.0/(double)(nz+1);
184
185 cnt = 1;
186 o_cnt = 1;
187 diag_i[0] = 0;
188 offd_i[0] = 0;
189 for (iz = 0; iz < nz_local; iz++)
190 {
191 for (iy = 0; iy < ny_local; iy++)
192 {
193 for (ix = 0; ix < nx_local; ix++)
194 {
195 diag_i[cnt] = diag_i[cnt-1];
196 offd_i[o_cnt] = offd_i[o_cnt-1];
197 diag_i[cnt]++;
198 if (iz)
199 diag_i[cnt]++;
200 else
201 {
202 if (ir)
203 {
204 offd_i[o_cnt]++;
205 }
206 }
207 if (iy)
208 diag_i[cnt]++;
209 else
210 {
211 if (iq)
212 {
213 offd_i[o_cnt]++;
214 }
215 }
216 if (ix)
217 diag_i[cnt]++;
218 else
219 {
220 if (ip)
221 {
222 offd_i[o_cnt]++;
223 }
224 }
225 if (ix+1 < nx_local)
226 diag_i[cnt]++;
227 else
228 {
229 if (ip+1 < P)
230 {
231 offd_i[o_cnt]++;
232 }
233 }
234 if (iy+1 < ny_local)
235 diag_i[cnt]++;
236 else
237 {
238 if (iq+1 < Q)
239 {
240 offd_i[o_cnt]++;
241 }
242 }
243 if (iz+1 < nz_local)
244 diag_i[cnt]++;
245 else
246 {
247 if (ir+1 < R)
248 {
249 offd_i[o_cnt]++;
250 }
251 }
252 cnt++;
253 o_cnt++;
254 }
255 }
256 }
257
258 diag_j = hypre_CTAlloc(int, diag_i[local_num_rows]);
259 diag_data = hypre_CTAlloc(double, diag_i[local_num_rows]);
260
261 if (num_procs > 1)
262 {
263 offd_j = hypre_CTAlloc(int, offd_i[local_num_rows]);
264 offd_data = hypre_CTAlloc(double, offd_i[local_num_rows]);
265 }
266
267 row_index = 0;
268 cnt = 0;
269 o_cnt = 0;
270 for (iz = 0; iz < nz_local; iz++)
271 {
272 zz = (double)(nz_part[ir]+iz+1)*hhz;
273 for (iy = 0; iy < ny_local; iy++)
274 {
275 yy = (double)(ny_part[iq]+iy+1)*hhy;
276 for (ix = 0; ix < nx_local; ix++)
277 {
278 xx = (double)(nx_part[ip]+ix+1)*hhx;
279 afp = eps*afun(xx+0.5*hhx,yy,zz)/hhx/hhx;
280 afm = eps*afun(xx-0.5*hhx,yy,zz)/hhx/hhx;
281 bfp = eps*bfun(xx,yy+0.5*hhy,zz)/hhy/hhy;
282 bfm = eps*bfun(xx,yy-0.5*hhy,zz)/hhy/hhy;
283 cfp = eps*cfun(xx,yy,zz+0.5*hhz)/hhz/hhz;
284 cfm = eps*cfun(xx,yy,zz-0.5*hhz)/hhz/hhz;
285 df = dfun(xx,yy,zz)/hhx;
286 ef = efun(xx,yy,zz)/hhy;
287 ff = ffun(xx,yy,zz)/hhz;
288 gf = gfun(xx,yy,zz);
289 diag_j[cnt] = row_index;
290 diag_data[cnt++] = afp+afm+bfp+bfm+cfp+cfm+gf-df-ef-ff;
291 rhs_data[row_index] = rfun(xx,yy,zz);
292 if (ix == 0 && ip == 0) rhs_data[row_index] += afm*bndfun(0,yy,zz);
293 if (iy == 0 && iq == 0) rhs_data[row_index] += bfm*bndfun(xx,0,zz);
294 if (iz == 0 && ir == 0) rhs_data[row_index] += cfm*bndfun(xx,yy,0);
295 if (ix+1 == nx_local && ip+1 == P)
296 rhs_data[row_index] += (afp-df)*bndfun(1.0,yy,zz);
297 if (iy+1 == ny_local && iq+1 == Q)
298 rhs_data[row_index] += (bfp-ef)*bndfun(xx,1.0,zz);
299 if (iz+1 == nz_local && ir+1 == R)
300 rhs_data[row_index] += (cfp-ff)*bndfun(xx,yy,1.0);
301 if (iz)
302 {
303 diag_j[cnt] = row_index-nx_local*ny_local;
304 diag_data[cnt++] = -cfm;
305 }
306 else
307 {
308 if (ir)
309 {
310 offd_j[o_cnt] = o_cnt;
311 hypre_map(ix,iy,iz-1,p,q,r-1,P,Q,R,
312 nx_part,ny_part,nz_part,global_part,
313 &col_map_offd[o_cnt]);
314 offd_data[o_cnt++] = -cfm;
315 }
316 }
317 if (iy)
318 {
319 diag_j[cnt] = row_index-nx_local;
320 diag_data[cnt++] = -bfm;
321 }
322 else
323 {
324 if (iq)
325 {
326 offd_j[o_cnt] = o_cnt;
327 hypre_map(ix,iy-1,iz,p,q-1,r,P,Q,R,
328 nx_part,ny_part,nz_part,global_part,
329 &col_map_offd[o_cnt]);
330 offd_data[o_cnt++] = -bfm;
331 }
332 }
333 if (ix)
334 {
335 diag_j[cnt] = row_index-1;
336 diag_data[cnt++] = -afm;
337 }
338 else
339 {
340 if (ip)
341 {
342 offd_j[o_cnt] = o_cnt;
343 hypre_map(ix-1,iy,iz,p-1,q,r,P,Q,R,
344 nx_part,ny_part,nz_part,global_part,
345 &col_map_offd[o_cnt]);
346 offd_data[o_cnt++] = -afm;
347 }
348 }
349 if (ix+1 < nx_local)
350 {
351 diag_j[cnt] = row_index+1;
352 diag_data[cnt++] = -afp+df;
353 }
354 else
355 {
356 if (ip+1 < P)
357 {
358 offd_j[o_cnt] = o_cnt;
359 hypre_map(ix+1,iy,iz,p+1,q,r,P,Q,R,
360 nx_part,ny_part,nz_part,global_part,
361 &col_map_offd[o_cnt]);
362 offd_data[o_cnt++] = -afp+df;
363 }
364 }
365 if (iy+1 < ny_local)
366 {
367 diag_j[cnt] = row_index+nx_local;
368 diag_data[cnt++] = -bfp +ef;
369 }
370 else
371 {
372 if (iq+1 < Q)
373 {
374 offd_j[o_cnt] = o_cnt;
375 hypre_map(ix,iy+1,iz,p,q+1,r,P,Q,R,
376 nx_part,ny_part,nz_part,global_part,
377 &col_map_offd[o_cnt]);
378 offd_data[o_cnt++] = -bfp+ef;
379 }
380 }
381 if (iz+1 < nz_local)
382 {
383 diag_j[cnt] = row_index+nx_local*ny_local;
384 diag_data[cnt++] = -cfp+ff;
385 }
386 else
387 {
388 if (ir+1 < R)
389 {
390 offd_j[o_cnt] = o_cnt;
391 hypre_map(ix,iy,iz+1,p,q,r+1,P,Q,R,
392 nx_part,ny_part,nz_part,global_part,
393 &col_map_offd[o_cnt]);
394 offd_data[o_cnt++] = -cfp+ff;
395 }
396 }
397 row_index++;
398 }
399 }
400 }
401
402 if (num_cols_offd) tmp_j = hypre_CTAlloc(int, num_cols_offd);
403 if (num_procs > 1)
404 {
405 for (i=0; i < num_cols_offd; i++)
406 tmp_j[i] = offd_j[i];
407
408 hypre_BigQsortbi(col_map_offd, tmp_j, 0, num_cols_offd-1);
409
410 for (i=0; i < num_cols_offd; i++)
411 for (j=0; j < num_cols_offd; j++)
412 if (offd_j[i] == tmp_j[j])
413 {
414 offd_j[i] = j;
415 break;
416 }
417 }
418
419 hypre_TFree(tmp_j);
420
421 par_rhs = hypre_ParVectorCreate(comm, grid_size, global_part_rhs);
422 rhs = hypre_ParVectorLocalVector(par_rhs);
423 hypre_VectorData(rhs) = rhs_data;
424
425 par_x = hypre_ParVectorCreate(comm, grid_size, global_part_x);
426 x = hypre_ParVectorLocalVector(par_x);
427 hypre_VectorData(x) = x_data;
428
429 A = hypre_ParCSRMatrixCreate(comm, grid_size, grid_size,
430 global_part, global_part, num_cols_offd,
431 diag_i[local_num_rows],
432 offd_i[local_num_rows]);
433
434 hypre_ParCSRMatrixColMapOffd(A) = col_map_offd;
435
436 diag = hypre_ParCSRMatrixDiag(A);
437 hypre_CSRMatrixI(diag) = diag_i;
438 hypre_CSRMatrixJ(diag) = diag_j;
439 hypre_CSRMatrixData(diag) = diag_data;
440
441 offd = hypre_ParCSRMatrixOffd(A);
442 hypre_CSRMatrixI(offd) = offd_i;
443 if (num_cols_offd)
444 {
445 hypre_CSRMatrixJ(offd) = offd_j;
446 hypre_CSRMatrixData(offd) = offd_data;
447 }
448
449 hypre_TFree(nx_part);
450 hypre_TFree(ny_part);
451 hypre_TFree(nz_part);
452
453 *rhs_ptr = (HYPRE_ParVector) par_rhs;
454 *x_ptr = (HYPRE_ParVector) par_x;
455
456 return (HYPRE_ParCSRMatrix) A;
457}
458
459double afun(double xx, double yy, double zz)
460{
461 double value;
462 /* value = 1.0 + 1000.0*fabs(xx-yy); */
463 if ((xx < 0.1 && yy < 0.1 && zz < 0.1)
464 || (xx < 0.1 && yy < 0.1 && zz > 0.9)
465 || (xx < 0.1 && yy > 0.9 && zz < 0.1)
466 || (xx > 0.9 && yy < 0.1 && zz < 0.1)
467 || (xx > 0.9 && yy > 0.9 && zz < 0.1)
468 || (xx > 0.9 && yy < 0.1 && zz > 0.9)
469 || (xx < 0.1 && yy > 0.9 && zz > 0.9)
470 || (xx > 0.9 && yy > 0.9 && zz > 0.9))
471 value = 0.01;
472 else if (xx >= 0.1 && xx <= 0.9
473 && yy >= 0.1 && yy <= 0.9
474 && zz >= 0.1 && zz <= 0.9)
475 value = 1000.0;
476 else
477 value = 1.0 ;
478 /* double value, pi;
479 pi = 4.0 * atan(1.0);
480 value = cos(pi*xx)*cos(pi*yy); */
481 return value;
482}
483
484double bfun(double xx, double yy, double zz)
485{
486 double value;
487 /* value = 1.0 + 1000.0*fabs(xx-yy); */
488 if ((xx < 0.1 && yy < 0.1 && zz < 0.1)
489 || (xx < 0.1 && yy < 0.1 && zz > 0.9)
490 || (xx < 0.1 && yy > 0.9 && zz < 0.1)
491 || (xx > 0.9 && yy < 0.1 && zz < 0.1)
492 || (xx > 0.9 && yy > 0.9 && zz < 0.1)
493 || (xx > 0.9 && yy < 0.1 && zz > 0.9)
494 || (xx < 0.1 && yy > 0.9 && zz > 0.9)
495 || (xx > 0.9 && yy > 0.9 && zz > 0.9))
496 value = 0.01;
497 else if (xx >= 0.1 && xx <= 0.9
498 && yy >= 0.1 && yy <= 0.9
499 && zz >= 0.1 && zz <= 0.9)
500 value = 1000.0;
501 else
502 value = 1.0 ;
503 /* double value, pi;
504 pi = 4.0 * atan(1.0);
505 value = 1.0 - 2.0*xx;
506 value = cos(pi*xx)*cos(pi*yy); */
507 /* double value;
508 value = 1.0 + 1000.0 * fabs(xx-yy);
509 double value, x0, y0;
510 x0 = fabs(xx - 0.5);
511 y0 = fabs(yy - 0.5);
512 if (y0 > x0) x0 = y0;
513 if (x0 >= 0.125 && x0 <= 0.25)
514 value = 1.0;
515 else
516 value = 1000.0;*/
517 return value;
518}
519
520double cfun(double xx, double yy, double zz)
521{
522 double value;
523 if ((xx < 0.1 && yy < 0.1 && zz < 0.1)
524 || (xx < 0.1 && yy < 0.1 && zz > 0.9)
525 || (xx < 0.1 && yy > 0.9 && zz < 0.1)
526 || (xx > 0.9 && yy < 0.1 && zz < 0.1)
527 || (xx > 0.9 && yy > 0.9 && zz < 0.1)
528 || (xx > 0.9 && yy < 0.1 && zz > 0.9)
529 || (xx < 0.1 && yy > 0.9 && zz > 0.9)
530 || (xx > 0.9 && yy > 0.9 && zz > 0.9))
531 value = 0.01;
532 else if (xx >= 0.1 && xx <= 0.9
533 && yy >= 0.1 && yy <= 0.9
534 && zz >= 0.1 && zz <= 0.9)
535 value = 1000.0;
536 else
537 value = 1.0 ;
538 /*if (xx <= 0.75 && yy <= 0.75 && zz <= 0.75)
539 value = 0.1;
540 else if (xx > 0.75 && yy > 0.75 && zz > 0.75)
541 value = 100000;
542 else
543 value = 1.0 ;*/
544 return value;
545}
546
547double dfun(double xx, double yy, double zz)
548{
549 double value;
550 /*double pi;
551 pi = 4.0 * atan(1.0);
552 value = -sin(pi*xx)*cos(pi*yy);*/
553 value = 0;
554 return value;
555}
556
557double efun(double xx, double yy, double zz)
558{
559 double value;
560 /*double pi;
561 pi = 4.0 * atan(1.0);
562 value = sin(pi*yy)*cos(pi*xx);*/
563 value = 0;
564 return value;
565}
566
567double ffun(double xx, double yy, double zz)
568{
569 double value;
570 value = 0.0;
571 return value;
572}
573
574double gfun(double xx, double yy, double zz)
575{
576 double value;
577 value = 0.0;
578 return value;
579}
580
581double rfun(double xx, double yy, double zz)
582{
583 /* double value, pi;
584 pi = 4.0 * atan(1.0);
585 value = -4.0*pi*pi*sin(pi*xx)*sin(pi*yy)*cos(pi*xx)*cos(pi*yy); */
586 double value;
587 /* value = xx*(1.0-xx)*yy*(1.0-yy); */
588 value = 1.0;
589 return value;
590}
591
592double bndfun(double xx, double yy, double zz)
593{
594 double value;
595 /*double pi;
596 pi = 4.0 * atan(1.0);
597 value = sin(pi*xx)+sin(13*pi*xx)+sin(pi*yy)+sin(13*pi*yy);*/
598 value = 0.0;
599 return value;
600}
Note: See TracBrowser for help on using the repository browser.