00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027 #include <stdlib.h>
00028 #include "../defs.h"
00029 #include "../athena.h"
00030 #include "../prototypes.h"
00031 #include "../globals.h"
00032
00033 #ifdef FFT_ENABLED
00034
00035 #ifdef FFT_BLOCK_DECOMP
00036
00037 #include "fft_3d.h"
00038 #include "fft_2d.h"
00039 #else
00040
00041 #include "fftw3.h"
00042 #endif
00043
00044
00045
00046
00047
00048
00049
00050
00051
00052
00053
00054
00055
00056 struct ath_3d_fft_plan *ath_3d_fft_quick_plan(DomainS *pD,
00057 ath_fft_data *data, ath_fft_direction dir)
00058 {
00059 GridS *pGrid = (pD->Grid);
00060
00061 int gnx1 = pD->Nx[0];
00062 int gnx2 = pD->Nx[1];
00063 int gnx3 = pD->Nx[2];
00064
00065
00066 int gis = pD->Disp[0];
00067 int gie = pD->Disp[0] + pGrid->Nx[0];
00068 int gjs = pD->Disp[1];
00069 int gje = pD->Disp[1] + pGrid->Nx[1];
00070 int gks = pD->Disp[2];
00071 int gke = pD->Disp[2] + pGrid->Nx[2];
00072
00073
00074
00075 return ath_3d_fft_create_plan(gnx3, gnx2, gnx1, gks, gke, gjs, gje,
00076 gis, gie, data, 0, dir);
00077 }
00078
00079
00080
00081
00082
00083
00084
00085
00086
00087
00088
00089
00090
00091
00092
00093
00094
00095 struct ath_3d_fft_plan *ath_3d_fft_create_plan(int gnx3, int gnx2,
00096 int gnx1, int gks, int gke, int gjs, int gje,
00097 int gis, int gie, ath_fft_data *data, int al,
00098 ath_fft_direction dir)
00099 {
00100 int nbuf, tmp;
00101 struct ath_3d_fft_plan *ath_plan;
00102
00103 if ((dir != ATH_FFT_FORWARD) && (dir != ATH_FFT_BACKWARD)) {
00104 ath_error("Invalid Athena FFT direction.\n");
00105 }
00106
00107
00108 ath_plan = (struct ath_3d_fft_plan *)malloc(sizeof(struct ath_3d_fft_plan));
00109 if (ath_plan == NULL) {
00110 ath_error("Couldn't malloc for FFT plan.");
00111 }
00112
00113 ath_plan->dir = dir;
00114
00115 ath_plan->cnt = (gke-gks+1)*(gje-gjs+1)*(gie-gis+1);
00116 ath_plan->gcnt = gnx3*gnx2*gnx1;
00117
00118 tmp = (al==0 ? 1 : 0);
00119 if (data != NULL) tmp = 0;
00120
00121
00122 if (data == NULL)
00123 data = (ath_fft_data *)ath_3d_fft_malloc(ath_plan);
00124 if (data == NULL)
00125 ath_error("Couln't malloc for FFT plan data.");
00126
00127
00128 #ifdef FFT_BLOCK_DECOMP
00129
00130 ath_plan->plan = fft_3d_create_plan(MPI_COMM_WORLD, gnx3, gnx2, gnx1,
00131 gks, gke, gjs, gje, gis, gie,
00132 gks, gke, gjs, gje, gis, gie,
00133 0, 0, &nbuf);
00134 #else
00135 if (dir == ATH_FFT_FORWARD) {
00136 ath_plan->plan = fftw_plan_dft_3d(gnx1, gnx2, gnx3, data, data,
00137 FFTW_FORWARD, FFTW_MEASURE);
00138 } else {
00139 ath_plan->plan = fftw_plan_dft_3d(gnx1, gnx2, gnx3, data, data,
00140 FFTW_BACKWARD, FFTW_MEASURE);
00141 }
00142 #endif
00143
00144 if (tmp) ath_3d_fft_free(data);
00145
00146 return ath_plan;
00147 }
00148
00149
00150
00151
00152
00153 ath_fft_data *ath_3d_fft_malloc(struct ath_3d_fft_plan *ath_plan)
00154 {
00155 return (ath_fft_data *)fftw_malloc(sizeof(ath_fft_data) * ath_plan->cnt);
00156 }
00157
00158
00159
00160
00161
00162 void ath_3d_fft(struct ath_3d_fft_plan *ath_plan, ath_fft_data *data)
00163 {
00164 #ifdef FFT_BLOCK_DECOMP
00165 fft_3d(data, data, ath_plan->dir, ath_plan->plan);
00166 #else
00167
00168 fftw_execute_dft(ath_plan->plan, data, data);
00169 #endif
00170
00171 return;
00172 }
00173
00174
00175
00176
00177
00178 void ath_3d_fft_free(ath_fft_data *data)
00179 {
00180 if (data != NULL) fftw_free((void*)data);
00181
00182 return;
00183 }
00184
00185
00186
00187
00188
00189 void ath_3d_fft_destroy_plan(struct ath_3d_fft_plan *ath_plan)
00190 {
00191 if (ath_plan != NULL) {
00192 #ifdef FFT_BLOCK_DECOMP
00193 fft_3d_destroy_plan(ath_plan->plan);
00194 #else
00195 fftw_destroy_plan(ath_plan->plan);
00196 #endif
00197 free(ath_plan);
00198 }
00199
00200 return;
00201 }
00202
00203
00204
00205
00206
00207
00208
00209
00210
00211
00212
00213
00214
00215 struct ath_2d_fft_plan *ath_2d_fft_quick_plan(DomainS *pD,
00216 ath_fft_data *data, ath_fft_direction dir)
00217 {
00218 GridS *pGrid=(pD->Grid);
00219 if (pGrid->Nx[2] != 1) {
00220 ath_error("ath_2d_fft_quick_plan only works for Nx3=1.\n");
00221 }
00222
00223
00224 int gnx1 = pD->Nx[0];
00225 int gnx2 = pD->Nx[1];
00226
00227
00228 int gis = pD->Disp[0];
00229 int gie = pD->Disp[0] + pGrid->Nx[0];
00230 int gjs = pD->Disp[1];
00231 int gje = pD->Disp[1] + pGrid->Nx[1];
00232
00233
00234
00235 return ath_2d_fft_create_plan(gnx2, gnx1, gjs, gje, gis, gie, data, 0, dir);
00236 }
00237
00238
00239
00240
00241
00242
00243
00244
00245
00246
00247
00248
00249
00250
00251 struct ath_2d_fft_plan *ath_2d_fft_create_plan(int gnx2, int gnx1,
00252 int gjs, int gje, int gis, int gie,
00253 ath_fft_data *data, int al,
00254 ath_fft_direction dir)
00255 {
00256 int nbuf, tmp;
00257 struct ath_2d_fft_plan *ath_plan;
00258
00259 if ((dir != ATH_FFT_FORWARD) && (dir != ATH_FFT_BACKWARD)) {
00260 ath_error("Invalid Athena FFT direction.\n");
00261 }
00262
00263
00264 ath_plan = (struct ath_2d_fft_plan *)malloc(sizeof(struct ath_2d_fft_plan));
00265 if (ath_plan == NULL) {
00266 ath_error("Couldn't malloc for FFT plan.");
00267 }
00268
00269 ath_plan->dir = dir;
00270
00271 ath_plan->cnt = (gje-gjs+1)*(gie-gis+1);
00272 ath_plan->gcnt = gnx2*gnx1;
00273
00274 tmp = (al==0 ? 1 : 0);
00275 if (data != NULL) tmp = 0;
00276
00277
00278 if (data == NULL)
00279 data = (ath_fft_data *)ath_2d_fft_malloc(ath_plan);
00280 if (data == NULL)
00281 ath_error("Couln't malloc for FFT plan data.");
00282
00283
00284 #ifdef FFT_BLOCK_DECOMP
00285
00286 ath_plan->plan = fft_2d_create_plan(MPI_COMM_WORLD, gnx2, gnx1, gjs, gje,
00287 gis, gie, gjs, gje, gis, gie,
00288 0, 0, &nbuf);
00289 #else
00290 if (dir == ATH_FFT_FORWARD) {
00291 ath_plan->plan = fftw_plan_dft_2d(gnx1, gnx2, data, data, FFTW_FORWARD,
00292 FFTW_MEASURE);
00293 } else {
00294 ath_plan->plan = fftw_plan_dft_2d(gnx1, gnx2, data, data, FFTW_BACKWARD,
00295 FFTW_MEASURE);
00296 }
00297 #endif
00298
00299 if (tmp) ath_2d_fft_free(data);
00300
00301 return ath_plan;
00302 }
00303
00304
00305
00306
00307
00308 ath_fft_data *ath_2d_fft_malloc(struct ath_2d_fft_plan *ath_plan)
00309 {
00310 return (ath_fft_data *)fftw_malloc(sizeof(ath_fft_data) * ath_plan->cnt);
00311 }
00312
00313
00314
00315
00316
00317
00318 void ath_2d_fft(struct ath_2d_fft_plan *ath_plan, fftw_complex *data)
00319 {
00320 #ifdef FFT_BLOCK_DECOMP
00321 fft_2d(data, data, ath_plan->dir, ath_plan->plan);
00322 #else
00323
00324 fftw_execute_dft(ath_plan->plan, data, data);
00325 #endif
00326
00327 return;
00328 }
00329
00330
00331
00332
00333
00334 void ath_2d_fft_free(ath_fft_data *data)
00335 {
00336 if (data != NULL) fftw_free((void*)data);
00337
00338 return;
00339 }
00340
00341
00342
00343
00344
00345 void ath_2d_fft_destroy_plan(struct ath_2d_fft_plan *ath_plan)
00346 {
00347 if (ath_plan != NULL) {
00348 #ifdef FFT_BLOCK_DECOMP
00349 fft_2d_destroy_plan(ath_plan->plan);
00350 #else
00351 fftw_destroy_plan(ath_plan->plan);
00352 #endif
00353 free(ath_plan);
00354 }
00355
00356 return;
00357 }
00358
00359 #endif