00001 #include "../copyright.h"
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021 #include <math.h>
00022 #include <float.h>
00023 #include "../defs.h"
00024 #include "../athena.h"
00025 #include "../globals.h"
00026 #include "prototypes.h"
00027 #include "../prototypes.h"
00028
00029 #ifdef SELF_GRAVITY_USING_FFT_OBC
00030
00031 #ifndef FFT_ENABLED
00032 #error self gravity with FFT requires configure --enable-fft
00033 #endif
00034
00035
00036 static struct ath_3d_fft_plan *fplan3d, *bplan3d;
00037 static ath_fft_data *work=NULL;
00038
00039
00040 #ifdef STATIC_MESH_REFINEMENT
00041 #error self gravity with FFT not yet implemented to work with SMR
00042 #endif
00043
00044
00045
00046
00047
00048
00049
00050 void selfg_fft_obc_3d(DomainS *pD)
00051 {
00052 GridS *pG = (pD->Grid);
00053 int i, is = pG->is, ie = pG->ie;
00054 int j, js = pG->js, je = pG->je;
00055 int k, ks = pG->ks, ke = pG->ke;
00056 int ip, jp, kp;
00057 int hNx1=pD->Nx[0]/2,hNx2=pD->Nx[1]/2,hNx3=pD->Nx[2];
00058 Real dx1sq=(pG->dx1*pG->dx1),dx2sq=(pG->dx2*pG->dx2),dx3sq=(pG->dx3*pG->dx3);
00059 Real dkx,dky,dkz;
00060 Real den;
00061 Real den_parent=0.;
00062
00063
00064 dkx = 2.0*PI/(double)(pD->Nx[0]);
00065 dky = 2.0*PI/(double)(pD->Nx[1]);
00066 dkz = 2.0*PI/(double)(pD->Nx[2]);
00067
00068
00069
00070
00071
00072 static int coeff_set=0;
00073 static Real ***Geee=NULL,***Gooo=NULL,***Goee=NULL,***Geoo=NULL,***Geoe=NULL, ***Goeo=NULL,***Geeo=NULL,***Gooe=NULL;
00074 if (!coeff_set){
00075
00076
00077
00078 if ((Geee = (Real***)calloc_3d_array(pG->Nx[0],pG->Nx[1],pG->Nx[2],sizeof(Real))) == NULL)
00079 ath_error("[poiss_coeff]: malloc returned a NULL pointer\n");
00080 if ((Gooo = (Real***)calloc_3d_array(pG->Nx[0],pG->Nx[1],pG->Nx[2],sizeof(Real))) == NULL)
00081 ath_error("[poiss_coeff]: malloc returned a NULL pointer\n");
00082 if ((Goee = (Real***)calloc_3d_array(pG->Nx[0],pG->Nx[1],pG->Nx[2],sizeof(Real))) == NULL)
00083 ath_error("[poiss_coeff]: malloc returned a NULL pointer\n");
00084 if ((Geoo = (Real***)calloc_3d_array(pG->Nx[0],pG->Nx[1],pG->Nx[2],sizeof(Real))) == NULL)
00085 ath_error("[poiss_coeff]: malloc returned a NULL pointer\n");
00086 if ((Geoe = (Real***)calloc_3d_array(pG->Nx[0],pG->Nx[1],pG->Nx[2],sizeof(Real))) == NULL)
00087 ath_error("[poiss_coeff]: malloc returned a NULL pointer\n");
00088 if ((Goeo = (Real***)calloc_3d_array(pG->Nx[0],pG->Nx[1],pG->Nx[2],sizeof(Real))) == NULL)
00089 ath_error("[poiss_coeff]: malloc returned a NULL pointer\n");
00090 if ((Geeo = (Real***)calloc_3d_array(pG->Nx[0],pG->Nx[1],pG->Nx[2],sizeof(Real))) == NULL)
00091 ath_error("[poiss_coeff]: malloc returned a NULL pointer\n");
00092 if ((Gooe = (Real***)calloc_3d_array(pG->Nx[0],pG->Nx[1],pG->Nx[2],sizeof(Real))) == NULL)
00093 ath_error("[poiss_coeff]: malloc returned a NULL pointer\n");
00094
00095
00096
00097
00098 for (i=is; i<=ie; i++){
00099 for (j=js; j<=je; j++){
00100 for (k=ks; k<=ke; k++){
00101 if ((k-ks+pG->Disp[2])==0 && (j-js+pG->Disp[1])==0 && (i-is+pG->Disp[0])==0)
00102 Geee[0][0][0] =0.0;
00103 else{
00104 Geee[i-is][j-js][k-ks] = 1.0/
00105 (((2.0*cos(( (i-is)+pG->Disp[0] )*dkx)-2.0)/dx1sq) +
00106 ((2.0*cos(( (j-js)+pG->Disp[1] )*dky)-2.0)/dx2sq) +
00107 ((2.0*cos(( (k-ks)+pG->Disp[2] )*dkz)-2.0)/dx3sq));
00108 }
00109 Gooo[i-is][j-js][k-ks] = 1.0/
00110 (((2.0*cos(( 0.5+(i-is)+pG->Disp[0] )*dkx)-2.0)/dx1sq) +
00111 ((2.0*cos(( 0.5+(j-js)+pG->Disp[1] )*dky)-2.0)/dx2sq) +
00112 ((2.0*cos(( 0.5+(k-ks)+pG->Disp[2] )*dkz)-2.0)/dx3sq));
00113 Goee[i-is][j-js][k-ks] = 1.0/
00114 (((2.0*cos(( 0.5+(i-is)+pG->Disp[0] )*dkx)-2.0)/dx1sq) +
00115 ((2.0*cos(( (j-js)+pG->Disp[1] )*dky)-2.0)/dx2sq) +
00116 ((2.0*cos(( (k-ks)+pG->Disp[2] )*dkz)-2.0)/dx3sq));
00117 Geoo[i-is][j-js][k-ks] = 1.0/
00118 (((2.0*cos(( (i-is)+pG->Disp[0] )*dkx)-2.0)/dx1sq) +
00119 ((2.0*cos(( 0.5+(j-js)+pG->Disp[1] )*dky)-2.0)/dx2sq) +
00120 ((2.0*cos(( 0.5+(k-ks)+pG->Disp[2] )*dkz)-2.0)/dx3sq));
00121 Geoe[i-is][j-js][k-ks] = 1.0/
00122 (((2.0*cos(( (i-is)+pG->Disp[0] )*dkx)-2.0)/dx1sq) +
00123 ((2.0*cos(( 0.5+(j-js)+pG->Disp[1] )*dky)-2.0)/dx2sq) +
00124 ((2.0*cos(( (k-ks)+pG->Disp[2] )*dkz)-2.0)/dx3sq));
00125 Goeo[i-is][j-js][k-ks] = 1.0/
00126 (((2.0*cos(( 0.5+(i-is)+pG->Disp[0] )*dkx)-2.0)/dx1sq) +
00127 ((2.0*cos(( (j-js)+pG->Disp[1] )*dky)-2.0)/dx2sq) +
00128 ((2.0*cos(( 0.5+(k-ks)+pG->Disp[2] )*dkz)-2.0)/dx3sq));
00129 Geeo[i-is][j-js][k-ks] = 1.0/
00130 (((2.0*cos(( (i-is)+pG->Disp[0] )*dkx)-2.0)/dx1sq) +
00131 ((2.0*cos(( (j-js)+pG->Disp[1] )*dky)-2.0)/dx2sq) +
00132 ((2.0*cos(( 0.5+(k-ks)+pG->Disp[2] )*dkz)-2.0)/dx3sq));
00133 Gooe[i-is][j-js][k-ks] = 1.0/
00134 (((2.0*cos(( 0.5+(i-is)+pG->Disp[0] )*dkx)-2.0)/dx1sq) +
00135 ((2.0*cos(( 0.5+(j-js)+pG->Disp[1] )*dky)-2.0)/dx2sq) +
00136 ((2.0*cos(( (k-ks)+pG->Disp[2] )*dkz)-2.0)/dx3sq));
00137 }
00138 }
00139 }
00140 coeff_set=1;
00141 }
00142
00143
00144
00145 for (k=ks-nghost; k<=ke+nghost; k++){
00146 for (j=js-nghost; j<=je+nghost; j++){
00147 for (i=is-nghost; i<=ie+nghost; i++){
00148 pG->Phi_old[k][j][i] = pG->Phi[k][j][i];
00149 }
00150 }
00151 }
00152
00153
00154
00155
00156
00157
00158
00159
00160
00161
00162
00163
00164
00165
00166
00167
00168
00169
00170
00171
00172
00173
00174
00175
00176
00177
00178
00179
00180
00181
00182
00183
00184
00185
00186 for (k=ks; k<=ke; k++){
00187 for (j=js; j<=je; j++){
00188 for (i=is; i<=ie; i++){
00189 work[F3DI(i-is,j-js,k-ks,pG->Nx[0],pG->Nx[1],pG->Nx[2])][0] =
00190 (pG->U[k][j][i].d - den_parent);
00191 work[F3DI(i-is,j-js,k-ks,pG->Nx[0],pG->Nx[1],pG->Nx[2])][1] = 0.0;
00192 }
00193 }
00194 }
00195
00196 ath_3d_fft(fplan3d, work);
00197
00198
00199 for (i=is; i<=ie; i++){
00200 for (j=js; j<=je; j++){
00201 for (k=ks; k<=ke; k++){
00202 work[F3DI(i-is,j-js,k-ks,pG->Nx[0],pG->Nx[1],pG->Nx[2])][0] *=
00203 Geee[i-is][j-js][k-ks];
00204 work[F3DI(i-is,j-js,k-ks,pG->Nx[0],pG->Nx[1],pG->Nx[2])][1] *=
00205 Geee[i-is][j-js][k-ks];
00206 }
00207 }
00208 }
00209
00210 ath_3d_fft(bplan3d, work);
00211
00212
00213
00214
00215 for (k=ks; k<=ke; k++){
00216 for (j=js; j<=je; j++){
00217 for (i=is; i<=ie; i++){
00218 pG->Phi[k][j][i] =
00219 work[F3DI(i-is,j-js,k-ks,pG->Nx[0],pG->Nx[1],pG->Nx[2])][0]
00220 *four_pi_G/(8.*bplan3d->gcnt);
00221 }
00222 }
00223 }
00224
00225
00226
00227
00228
00229 for (k=ks; k<=ke; k++){
00230 for (j=js; j<=je; j++){
00231 for (i=is; i<=ie; i++){
00232 work[F3DI(i-is,j-js,k-ks,pG->Nx[0],pG->Nx[1],pG->Nx[2])][0] =
00233 cos(0.5*(((i-is)+pG->Disp[0])*dkx+
00234 ((j-js)+pG->Disp[1])*dky+
00235 ((k-ks)+pG->Disp[2])*dkz))*
00236 (pG->U[k][j][i].d - den_parent);
00237 work[F3DI(i-is,j-js,k-ks,pG->Nx[0],pG->Nx[1],pG->Nx[2])][1] =
00238 sin(0.5*(((i-is)+pG->Disp[0])*dkx+
00239 ((j-js)+pG->Disp[1])*dky+
00240 ((k-ks)+pG->Disp[2])*dkz))*
00241 (pG->U[k][j][i].d - den_parent);
00242 }
00243 }
00244 }
00245
00246 ath_3d_fft(fplan3d, work);
00247
00248
00249 for (i=is; i<=ie; i++){
00250 for (j=js; j<=je; j++){
00251 for (k=ks; k<=ke; k++){
00252 work[F3DI(i-is,j-js,k-ks,pG->Nx[0],pG->Nx[1],pG->Nx[2])][0] *=
00253 Gooo[i-is][j-js][k-ks];
00254 work[F3DI(i-is,j-js,k-ks,pG->Nx[0],pG->Nx[1],pG->Nx[2])][1] *=
00255 Gooo[i-is][j-js][k-ks];
00256 }
00257 }
00258 }
00259
00260 ath_3d_fft(bplan3d, work);
00261
00262
00263
00264
00265 for (k=ks; k<=ke; k++){
00266 for (j=js; j<=je; j++){
00267 for (i=is; i<=ie; i++){
00268 pG->Phi[k][j][i] += (
00269 work[F3DI(i-is,j-js,k-ks,pG->Nx[0],pG->Nx[1],pG->Nx[2])][0]
00270 *cos(0.5*(((i-is)+pG->Disp[0])*dkx+
00271 ((j-js)+pG->Disp[1])*dky+
00272 ((k-ks)+pG->Disp[2])*dkz))
00273 + work[F3DI(i-is,j-js,k-ks,pG->Nx[0],pG->Nx[1],pG->Nx[2])][1]
00274 *sin(0.5*(((i-is)+pG->Disp[0])*dkx+
00275 ((j-js)+pG->Disp[1])*dky+
00276 ((k-ks)+pG->Disp[2])*dkz))
00277 )*four_pi_G/(8.*bplan3d->gcnt);
00278
00279 }
00280 }
00281 }
00282
00283
00284
00285
00286 for (k=ks; k<=ke; k++){
00287 for (j=js; j<=je; j++){
00288 for (i=is; i<=ie; i++){
00289 work[F3DI(i-is,j-js,k-ks,pG->Nx[0],pG->Nx[1],pG->Nx[2])][0] =
00290 cos(0.5*((i-is)+pG->Disp[0])*dkx)*
00291 (pG->U[k][j][i].d - den_parent);
00292 work[F3DI(i-is,j-js,k-ks,pG->Nx[0],pG->Nx[1],pG->Nx[2])][1] =
00293 sin(0.5*((i-is)+pG->Disp[0])*dkx)*
00294 (pG->U[k][j][i].d - den_parent);
00295 }
00296 }
00297 }
00298
00299 ath_3d_fft(fplan3d, work);
00300
00301
00302 for (i=is; i<=ie; i++){
00303 for (j=js; j<=je; j++){
00304 for (k=ks; k<=ke; k++){
00305 work[F3DI(i-is,j-js,k-ks,pG->Nx[0],pG->Nx[1],pG->Nx[2])][0] *=
00306 Goee[i-is][j-js][k-ks];
00307 work[F3DI(i-is,j-js,k-ks,pG->Nx[0],pG->Nx[1],pG->Nx[2])][1] *=
00308 Goee[i-is][j-js][k-ks];
00309 }
00310 }
00311 }
00312
00313 ath_3d_fft(bplan3d, work);
00314
00315
00316
00317
00318 for (k=ks; k<=ke; k++){
00319 for (j=js; j<=je; j++){
00320 for (i=is; i<=ie; i++){
00321 pG->Phi[k][j][i] += (
00322 work[F3DI(i-is,j-js,k-ks,pG->Nx[0],pG->Nx[1],pG->Nx[2])][0]
00323 *cos(0.5*((i-is)+pG->Disp[0])*dkx)
00324 + work[F3DI(i-is,j-js,k-ks,pG->Nx[0],pG->Nx[1],pG->Nx[2])][1]
00325 *sin(0.5*((i-is)+pG->Disp[0])*dkx)
00326 )*four_pi_G/(8.*bplan3d->gcnt);
00327
00328 }
00329 }
00330 }
00331
00332
00333
00334
00335 for (k=ks; k<=ke; k++){
00336 for (j=js; j<=je; j++){
00337 for (i=is; i<=ie; i++){
00338 work[F3DI(i-is,j-js,k-ks,pG->Nx[0],pG->Nx[1],pG->Nx[2])][0] =
00339 cos(0.5*(((j-js)+pG->Disp[1])*dky+
00340 ((k-ks)+pG->Disp[2])*dkz))*
00341 (pG->U[k][j][i].d - den_parent);
00342 work[F3DI(i-is,j-js,k-ks,pG->Nx[0],pG->Nx[1],pG->Nx[2])][1] =
00343 sin(0.5*(((j-js)+pG->Disp[1])*dky+
00344 ((k-ks)+pG->Disp[2])*dkz))*
00345 (pG->U[k][j][i].d - den_parent);
00346 }
00347 }
00348 }
00349
00350 ath_3d_fft(fplan3d, work);
00351
00352
00353 for (i=is; i<=ie; i++){
00354 for (j=js; j<=je; j++){
00355 for (k=ks; k<=ke; k++){
00356 work[F3DI(i-is,j-js,k-ks,pG->Nx[0],pG->Nx[1],pG->Nx[2])][0] *=
00357 Geoo[i-is][j-js][k-ks];
00358 work[F3DI(i-is,j-js,k-ks,pG->Nx[0],pG->Nx[1],pG->Nx[2])][1] *=
00359 Geoo[i-is][j-js][k-ks];
00360 }
00361 }
00362 }
00363
00364 ath_3d_fft(bplan3d, work);
00365
00366
00367
00368
00369 for (k=ks; k<=ke; k++){
00370 for (j=js; j<=je; j++){
00371 for (i=is; i<=ie; i++){
00372 pG->Phi[k][j][i] += (
00373 work[F3DI(i-is,j-js,k-ks,pG->Nx[0],pG->Nx[1],pG->Nx[2])][0]
00374 *cos(0.5*(((j-js)+pG->Disp[1])*dky+
00375 ((k-ks)+pG->Disp[2])*dkz))
00376 + work[F3DI(i-is,j-js,k-ks,pG->Nx[0],pG->Nx[1],pG->Nx[2])][1]
00377 *sin(0.5*(((j-js)+pG->Disp[1])*dky+
00378 ((k-ks)+pG->Disp[2])*dkz))
00379 )*four_pi_G/(8.*bplan3d->gcnt);
00380
00381 }
00382 }
00383 }
00384
00385
00386
00387
00388
00389 for (k=ks; k<=ke; k++){
00390 for (j=js; j<=je; j++){
00391 for (i=is; i<=ie; i++){
00392 work[F3DI(i-is,j-js,k-ks,pG->Nx[0],pG->Nx[1],pG->Nx[2])][0] =
00393 cos(0.5*((j-js)+pG->Disp[1])*dky)*
00394 (pG->U[k][j][i].d - den_parent);
00395 work[F3DI(i-is,j-js,k-ks,pG->Nx[0],pG->Nx[1],pG->Nx[2])][1] =
00396 sin(0.5*((j-js)+pG->Disp[1])*dky)*
00397 (pG->U[k][j][i].d - den_parent);
00398 }
00399 }
00400 }
00401
00402 ath_3d_fft(fplan3d, work);
00403
00404
00405 for (i=is; i<=ie; i++){
00406 for (j=js; j<=je; j++){
00407 for (k=ks; k<=ke; k++){
00408 work[F3DI(i-is,j-js,k-ks,pG->Nx[0],pG->Nx[1],pG->Nx[2])][0] *=
00409 Geoe[i-is][j-js][k-ks];
00410 work[F3DI(i-is,j-js,k-ks,pG->Nx[0],pG->Nx[1],pG->Nx[2])][1] *=
00411 Geoe[i-is][j-js][k-ks];
00412 }
00413 }
00414 }
00415
00416 ath_3d_fft(bplan3d, work);
00417
00418
00419
00420
00421 for (k=ks; k<=ke; k++){
00422 for (j=js; j<=je; j++){
00423 for (i=is; i<=ie; i++){
00424 pG->Phi[k][j][i] += (
00425 work[F3DI(i-is,j-js,k-ks,pG->Nx[0],pG->Nx[1],pG->Nx[2])][0]
00426 *cos(0.5*((j-js)+pG->Disp[1])*dky)
00427 + work[F3DI(i-is,j-js,k-ks,pG->Nx[0],pG->Nx[1],pG->Nx[2])][1]
00428 *sin(0.5*((j-js)+pG->Disp[1])*dky)
00429 )*four_pi_G/(8.*bplan3d->gcnt);
00430
00431 }
00432 }
00433 }
00434
00435
00436
00437
00438
00439 for (k=ks; k<=ke; k++){
00440 for (j=js; j<=je; j++){
00441 for (i=is; i<=ie; i++){
00442 work[F3DI(i-is,j-js,k-ks,pG->Nx[0],pG->Nx[1],pG->Nx[2])][0] =
00443 cos(0.5*(((i-is)+pG->Disp[0])*dkx+
00444 ((k-ks)+pG->Disp[2])*dkz))*
00445 (pG->U[k][j][i].d - den_parent);
00446 work[F3DI(i-is,j-js,k-ks,pG->Nx[0],pG->Nx[1],pG->Nx[2])][1] =
00447 sin(0.5*(((i-is)+pG->Disp[0])*dkx+
00448 ((k-ks)+pG->Disp[2])*dkz))*
00449 (pG->U[k][j][i].d - den_parent);
00450 }
00451 }
00452 }
00453
00454 ath_3d_fft(fplan3d, work);
00455
00456
00457 for (i=is; i<=ie; i++){
00458 for (j=js; j<=je; j++){
00459 for (k=ks; k<=ke; k++){
00460 work[F3DI(i-is,j-js,k-ks,pG->Nx[0],pG->Nx[1],pG->Nx[2])][0] *=
00461 Goeo[i-is][j-js][k-ks];
00462 work[F3DI(i-is,j-js,k-ks,pG->Nx[0],pG->Nx[1],pG->Nx[2])][1] *=
00463 Goeo[i-is][j-js][k-ks];
00464 }
00465 }
00466 }
00467
00468 ath_3d_fft(bplan3d, work);
00469
00470
00471
00472
00473 for (k=ks; k<=ke; k++){
00474 for (j=js; j<=je; j++){
00475 for (i=is; i<=ie; i++){
00476 pG->Phi[k][j][i] += (
00477 work[F3DI(i-is,j-js,k-ks,pG->Nx[0],pG->Nx[1],pG->Nx[2])][0]
00478 *cos(0.5*(((i-is)+pG->Disp[0])*dkx+
00479 ((k-ks)+pG->Disp[2])*dkz))
00480 + work[F3DI(i-is,j-js,k-ks,pG->Nx[0],pG->Nx[1],pG->Nx[2])][1]
00481 *sin(0.5*(((i-is)+pG->Disp[0])*dkx+
00482 ((k-ks)+pG->Disp[2])*dkz))
00483 )*four_pi_G/(8.*bplan3d->gcnt);
00484
00485 }
00486 }
00487 }
00488
00489
00490
00491
00492
00493 for (k=ks; k<=ke; k++){
00494 for (j=js; j<=je; j++){
00495 for (i=is; i<=ie; i++){
00496 work[F3DI(i-is,j-js,k-ks,pG->Nx[0],pG->Nx[1],pG->Nx[2])][0] =
00497 cos(0.5*((k-ks)+pG->Disp[2])*dkz)*
00498 (pG->U[k][j][i].d - den_parent);
00499 work[F3DI(i-is,j-js,k-ks,pG->Nx[0],pG->Nx[1],pG->Nx[2])][1] =
00500 sin(0.5*((k-ks)+pG->Disp[2])*dkz)*
00501 (pG->U[k][j][i].d - den_parent);
00502 }
00503 }
00504 }
00505
00506 ath_3d_fft(fplan3d, work);
00507
00508
00509 for (i=is; i<=ie; i++){
00510 for (j=js; j<=je; j++){
00511 for (k=ks; k<=ke; k++){
00512 work[F3DI(i-is,j-js,k-ks,pG->Nx[0],pG->Nx[1],pG->Nx[2])][0] *=
00513 Geeo[i-is][j-js][k-ks];
00514 work[F3DI(i-is,j-js,k-ks,pG->Nx[0],pG->Nx[1],pG->Nx[2])][1] *=
00515 Geeo[i-is][j-js][k-ks];
00516 }
00517 }
00518 }
00519
00520 ath_3d_fft(bplan3d, work);
00521
00522
00523
00524
00525 for (k=ks; k<=ke; k++){
00526 for (j=js; j<=je; j++){
00527 for (i=is; i<=ie; i++){
00528 pG->Phi[k][j][i] += (
00529 work[F3DI(i-is,j-js,k-ks,pG->Nx[0],pG->Nx[1],pG->Nx[2])][0]
00530 *cos(0.5*((k-ks)+pG->Disp[2])*dkz)
00531 + work[F3DI(i-is,j-js,k-ks,pG->Nx[0],pG->Nx[1],pG->Nx[2])][1]
00532 *sin(0.5*((k-ks)+pG->Disp[2])*dkz)
00533 )*four_pi_G/(8.*bplan3d->gcnt);
00534
00535 }
00536 }
00537 }
00538
00539
00540
00541
00542
00543 for (k=ks; k<=ke; k++){
00544 for (j=js; j<=je; j++){
00545 for (i=is; i<=ie; i++){
00546 work[F3DI(i-is,j-js,k-ks,pG->Nx[0],pG->Nx[1],pG->Nx[2])][0] =
00547 cos(0.5*(((i-is)+pG->Disp[0])*dkx+
00548 ((j-js)+pG->Disp[1])*dky))*
00549 (pG->U[k][j][i].d - den_parent);
00550 work[F3DI(i-is,j-js,k-ks,pG->Nx[0],pG->Nx[1],pG->Nx[2])][1] =
00551 sin(0.5*(((i-is)+pG->Disp[0])*dkx+
00552 ((j-js)+pG->Disp[1])*dky))*
00553 (pG->U[k][j][i].d - den_parent);
00554 }
00555 }
00556 }
00557
00558 ath_3d_fft(fplan3d, work);
00559
00560
00561 for (i=is; i<=ie; i++){
00562 for (j=js; j<=je; j++){
00563 for (k=ks; k<=ke; k++){
00564 work[F3DI(i-is,j-js,k-ks,pG->Nx[0],pG->Nx[1],pG->Nx[2])][0] *=
00565 Gooe[i-is][j-js][k-ks];
00566 work[F3DI(i-is,j-js,k-ks,pG->Nx[0],pG->Nx[1],pG->Nx[2])][1] *=
00567 Gooe[i-is][j-js][k-ks];
00568 }
00569 }
00570 }
00571
00572 ath_3d_fft(bplan3d, work);
00573
00574
00575
00576
00577 for (k=ks; k<=ke; k++){
00578 for (j=js; j<=je; j++){
00579 for (i=is; i<=ie; i++){
00580 pG->Phi[k][j][i] += (
00581 work[F3DI(i-is,j-js,k-ks,pG->Nx[0],pG->Nx[1],pG->Nx[2])][0]
00582 *cos(0.5*(((i-is)+pG->Disp[0])*dkx+
00583 ((j-js)+pG->Disp[1])*dky))
00584 + work[F3DI(i-is,j-js,k-ks,pG->Nx[0],pG->Nx[1],pG->Nx[2])][1]
00585 *sin(0.5*(((i-is)+pG->Disp[0])*dkx+
00586 ((j-js)+pG->Disp[1])*dky))
00587 )*four_pi_G/(8.*bplan3d->gcnt);
00588
00589 }
00590 }
00591 }
00592
00593 return;
00594 }
00595
00596
00597
00598
00599
00600
00601
00602
00603
00604
00605
00606
00607 void selfg_fft_3d_init(MeshS *pM)
00608 {
00609 DomainS *pD;
00610 int nl,nd;
00611 for (nl=0; nl<(pM->NLevels); nl++){
00612 for (nd=0; nd<(pM->DomainsPerLevel[nl]); nd++){
00613 if (pM->Domain[nl][nd].Grid != NULL){
00614 pD = (DomainS*)&(pM->Domain[nl][nd]);
00615 fplan3d = ath_3d_fft_quick_plan(pD, NULL, ATH_FFT_FORWARD);
00616 bplan3d = ath_3d_fft_quick_plan(pD, NULL, ATH_FFT_BACKWARD);
00617 work = ath_3d_fft_malloc(fplan3d);
00618 }
00619 }
00620 }
00621 }
00622
00623 #endif