Actual source code: ex28.c
petsc-3.13.6 2020-09-29
1: static char help[] = "Test sequential USFFT interface on a 3-dof field over a uniform DMDA and compares to the result of FFTW acting on a split version of the field\n\n";
3: /*
4: Compiling the code:
5: This code uses the complex numbers version of PETSc and the FFTW package, so configure
6: must be run to enable these.
8: */
10: #define DOF 3
12: #include <petscmat.h>
13: #include <petscdm.h>
14: #include <petscdmda.h>
15: PetscInt main(PetscInt argc,char **args)
16: {
17: typedef enum {RANDOM, CONSTANT, TANH, NUM_FUNCS} FuncType;
18: const char *funcNames[NUM_FUNCS] = {"random", "constant", "tanh"};
19: Mat A, AA;
20: PetscMPIInt size;
21: PetscInt N,i, stencil=1,dof=3;
22: PetscInt dim[3] = {10,10,10}, ndim = 3;
23: Vec coords,x,y,z,xx, yy, zz;
24: Vec xxsplit[DOF], yysplit[DOF], zzsplit[DOF];
25: PetscReal h[3];
26: PetscScalar s;
27: PetscRandom rdm;
28: PetscReal norm, enorm;
29: PetscInt func,ii;
30: FuncType function = TANH;
31: DM da, da1, coordsda;
32: PetscBool view_x = PETSC_FALSE, view_y = PETSC_FALSE, view_z = PETSC_FALSE;
35: PetscInitialize(&argc,&args,(char*)0,help);if (ierr) return ierr;
36: MPI_Comm_size(PETSC_COMM_WORLD, &size);
37: if (size != 1) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_SUP, "This is a uniprocessor example only!");
38: PetscOptionsBegin(PETSC_COMM_WORLD, NULL, "USFFT Options", "ex27");
39: PetscOptionsEList("-function", "Function type", "ex27", funcNames, NUM_FUNCS, funcNames[function], &func, NULL);
40: function = (FuncType) func;
41: PetscOptionsEnd();
42: PetscOptionsGetBool(NULL,NULL,"-view_x",&view_x,NULL);
43: PetscOptionsGetBool(NULL,NULL,"-view_y",&view_y,NULL);
44: PetscOptionsGetBool(NULL,NULL,"-view_z",&view_z,NULL);
45: PetscOptionsGetIntArray(NULL,NULL,"-dim",dim,&ndim,NULL);
47: /* DMDA with the correct fiber dimension */
48: DMDACreate3d(PETSC_COMM_SELF,DM_BOUNDARY_NONE,DM_BOUNDARY_NONE,DM_BOUNDARY_NONE,DMDA_STENCIL_STAR,dim[0],dim[1],dim[2],PETSC_DECIDE,PETSC_DECIDE,PETSC_DECIDE,
49: dof,stencil,NULL,NULL,NULL,&da);
50: DMSetFromOptions(da);
51: DMSetUp(da);
52: /* DMDA with fiber dimension 1 for split fields */
53: DMDACreate3d(PETSC_COMM_SELF,DM_BOUNDARY_NONE,DM_BOUNDARY_NONE,DM_BOUNDARY_NONE,DMDA_STENCIL_STAR,dim[0],dim[1],dim[2],PETSC_DECIDE,PETSC_DECIDE,PETSC_DECIDE,
54: 1,stencil,NULL,NULL,NULL,&da1);
55: DMSetFromOptions(da1);
56: DMSetUp(da1);
58: /* Coordinates */
59: DMGetCoordinateDM(da,&coordsda);
60: DMGetGlobalVector(coordsda,&coords);
61: PetscObjectSetName((PetscObject) coords,"Grid coordinates");
62: for (i = 0, N = 1; i < 3; i++) {
63: h[i] = 1.0/dim[i];
64: PetscScalar *a;
65: VecGetArray(coords, &a);
66: PetscInt j,k,n = 0;
67: for (i = 0; i < 3; ++i) {
68: for (j = 0; j < dim[i]; ++j) {
69: for (k = 0; k < 3; ++k) {
70: a[n] = j*h[i]; /* coordinate along the j-th point in the i-th dimension */
71: ++n;
72: }
73: }
74: }
75: VecRestoreArray(coords, &a);
77: }
78: DMSetCoordinates(da, coords);
79: VecDestroy(&coords);
81: /* Work vectors */
82: DMGetGlobalVector(da, &x);
83: PetscObjectSetName((PetscObject) x, "Real space vector");
84: DMGetGlobalVector(da, &xx);
85: PetscObjectSetName((PetscObject) xx, "Real space vector");
86: DMGetGlobalVector(da, &y);
87: PetscObjectSetName((PetscObject) y, "USFFT frequency space vector");
88: DMGetGlobalVector(da, &yy);
89: PetscObjectSetName((PetscObject) yy, "FFTW frequency space vector");
90: DMGetGlobalVector(da, &z);
91: PetscObjectSetName((PetscObject) z, "USFFT reconstructed vector");
92: DMGetGlobalVector(da, &zz);
93: PetscObjectSetName((PetscObject) zz, "FFTW reconstructed vector");
94: /* Split vectors for FFTW */
95: for (ii = 0; ii < 3; ++ii) {
96: DMGetGlobalVector(da1, &xxsplit[ii]);
97: PetscObjectSetName((PetscObject) xxsplit[ii], "Real space split vector");
98: DMGetGlobalVector(da1, &yysplit[ii]);
99: PetscObjectSetName((PetscObject) yysplit[ii], "FFTW frequency space split vector");
100: DMGetGlobalVector(da1, &zzsplit[ii]);
101: PetscObjectSetName((PetscObject) zzsplit[ii], "FFTW reconstructed split vector");
102: }
105: PetscPrintf(PETSC_COMM_SELF, "%3-D: USFFT on vector of ");
106: for (i = 0, N = 1; i < 3; i++) {
107: PetscPrintf(PETSC_COMM_SELF, "dim[%d] = %d ",i,dim[i]);
108: N *= dim[i];
109: }
110: PetscPrintf(PETSC_COMM_SELF, "; total size %d \n",N);
113: if (function == RANDOM) {
114: PetscRandomCreate(PETSC_COMM_SELF, &rdm);
115: PetscRandomSetFromOptions(rdm);
116: VecSetRandom(x, rdm);
117: PetscRandomDestroy(&rdm);
118: } else if (function == CONSTANT) {
119: VecSet(x, 1.0);
120: } else if (function == TANH) {
121: PetscScalar *a;
122: VecGetArray(x, &a);
123: PetscInt j,k = 0;
124: for (i = 0; i < 3; ++i) {
125: for (j = 0; j < dim[i]; ++j) {
126: a[k] = tanh((j - dim[i]/2.0)*(10.0/dim[i]));
127: ++k;
128: }
129: }
130: VecRestoreArray(x, &a);
131: }
132: if (view_x) {
133: VecView(x, PETSC_VIEWER_STDOUT_WORLD);
134: }
135: VecCopy(x,xx);
136: /* Split xx */
137: VecStrideGatherAll(xx,xxsplit, INSERT_VALUES); /*YES! 'Gather' means 'split' (or maybe 'scatter'?)! */
139: VecNorm(x,NORM_2,&norm);
140: PetscPrintf(PETSC_COMM_SELF, "|x|_2 = %g\n",norm);
142: /* create USFFT object */
143: MatCreateSeqUSFFT(da,da,&A);
144: /* create FFTW object */
145: MatCreateSeqFFTW(PETSC_COMM_SELF,3,dim,&AA);
147: /* apply USFFT and FFTW FORWARD "preemptively", so the fftw_plans can be reused on different vectors */
148: MatMult(A,x,z);
149: for (ii = 0; ii < 3; ++ii) {
150: MatMult(AA,xxsplit[ii],zzsplit[ii]);
151: }
152: /* Now apply USFFT and FFTW forward several (3) times */
153: for (i=0; i<3; ++i) {
154: MatMult(A,x,y);
155: for (ii = 0; ii < 3; ++ii) {
156: MatMult(AA,xxsplit[ii],yysplit[ii]);
157: }
158: MatMultTranspose(A,y,z);
159: for (ii = 0; ii < 3; ++ii) {
160: MatMult(AA,yysplit[ii],zzsplit[ii]);
161: }
162: }
163: /* Unsplit yy */
164: VecStrideScatterAll(yysplit, yy, INSERT_VALUES); /*YES! 'Scatter' means 'collect' (or maybe 'gather'?)! */
165: /* Unsplit zz */
166: VecStrideScatterAll(zzsplit, zz, INSERT_VALUES); /*YES! 'Scatter' means 'collect' (or maybe 'gather'?)! */
168: if (view_y) {
169: PetscPrintf(PETSC_COMM_WORLD, "y = \n");
170: VecView(y, PETSC_VIEWER_STDOUT_WORLD);
171: PetscPrintf(PETSC_COMM_WORLD, "yy = \n");
172: VecView(yy, PETSC_VIEWER_STDOUT_WORLD);
173: }
175: if (view_z) {
176: PetscPrintf(PETSC_COMM_WORLD, "z = \n");
177: VecView(z, PETSC_VIEWER_STDOUT_WORLD);
178: PetscPrintf(PETSC_COMM_WORLD, "zz = \n");
179: VecView(zz, PETSC_VIEWER_STDOUT_WORLD);
180: }
182: /* compare x and z. USFFT computes an unnormalized DFT, thus z = N*x */
183: s = 1.0/(PetscReal)N;
184: VecScale(z,s);
185: VecAXPY(x,-1.0,z);
186: VecNorm(x,NORM_1,&enorm);
187: PetscPrintf(PETSC_COMM_SELF, "|x-z| = %g\n",enorm);
189: /* compare xx and zz. FFTW computes an unnormalized DFT, thus zz = N*x */
190: s = 1.0/(PetscReal)N;
191: VecScale(zz,s);
192: VecAXPY(xx,-1.0,zz);
193: VecNorm(xx,NORM_1,&enorm);
194: PetscPrintf(PETSC_COMM_SELF, "|xx-zz| = %g\n",enorm);
196: /* compare y and yy: USFFT and FFTW results*/
197: VecNorm(y,NORM_2,&norm);
198: VecAXPY(y,-1.0,yy);
199: VecNorm(y,NORM_1,&enorm);
200: PetscPrintf(PETSC_COMM_SELF, "|y|_2 = %g\n",norm);
201: PetscPrintf(PETSC_COMM_SELF, "|y-yy| = %g\n",enorm);
203: /* compare z and zz: USFFT and FFTW results*/
204: VecNorm(z,NORM_2,&norm);
205: VecAXPY(z,-1.0,zz);
206: VecNorm(z,NORM_1,&enorm);
207: PetscPrintf(PETSC_COMM_SELF, "|z|_2 = %g\n",norm);
208: PetscPrintf(PETSC_COMM_SELF, "|z-zz| = %g\n",enorm);
211: /* free spaces */
212: DMRestoreGlobalVector(da,&x);
213: DMRestoreGlobalVector(da,&xx);
214: DMRestoreGlobalVector(da,&y);
215: DMRestoreGlobalVector(da,&yy);
216: DMRestoreGlobalVector(da,&z);
217: DMRestoreGlobalVector(da,&zz);
219: PetscFinalize();
220: return ierr;
221: }