00001 #include "../copyright.h"
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 #include <math.h>
00027 #include <float.h>
00028 #include "../defs.h"
00029 #include "../athena.h"
00030 #include "../globals.h"
00031 #include "prototypes.h"
00032 #include "../prototypes.h"
00033
00034 #ifdef SELF_GRAVITY_USING_FFT
00035
00036 #ifndef FFT_ENABLED
00037 #error self gravity with FFT requires configure --enable-fft
00038 #endif
00039
00040
00041 static struct ath_2d_fft_plan *fplan2d, *bplan2d;
00042 static struct ath_3d_fft_plan *fplan3d, *bplan3d;
00043 static ath_fft_data *work=NULL;
00044
00045 #ifdef STATIC_MESH_REFINEMENT
00046 #error self gravity with FFT not yet implemented to work with SMR
00047 #endif
00048
00049
00050
00051
00052
00053
00054
00055
00056
00057 void selfg_fft_1d(DomainS *pD)
00058 {
00059 GridS *pG = (pD->Grid);
00060 int i, is = pG->is, ie = pG->ie;
00061 int js = pG->js;
00062 int ks = pG->ks;
00063 Real total_Phi=0.0,drho,dx_sq = (pG->dx1*pG->dx1);
00064
00065
00066
00067 for (i=is-nghost; i<=ie+nghost; i++){
00068 pG->Phi_old[ks][js][i] = pG->Phi[ks][js][i];
00069 }
00070
00071
00072
00073 pG->Phi[ks][js][is] = 0.0;
00074 for (i=is; i<=ie; i++) {
00075 drho = (pG->U[ks][js][i].d - grav_mean_rho);
00076 pG->Phi[ks][js][is] += ((float)(i-is+1))*four_pi_G*dx_sq*drho;
00077 }
00078 pG->Phi[ks][js][is] /= (float)(pG->Nx[0]);
00079
00080 drho = (pG->U[ks][js][is].d - grav_mean_rho);
00081 pG->Phi[ks][js][is+1] = 2.0*pG->Phi[ks][js][is] + four_pi_G*dx_sq*drho;
00082 for (i=is+2; i<=ie; i++) {
00083 drho = (pG->U[ks][js][i-1].d - grav_mean_rho);
00084 pG->Phi[ks][js][i] = four_pi_G*dx_sq*drho
00085 + 2.0*pG->Phi[ks][js][i-1] - pG->Phi[ks][js][i-2];
00086 }
00087
00088
00089
00090 for (i=is; i<=ie; i++) {
00091 total_Phi += pG->Phi[ks][js][i];
00092 }
00093 total_Phi /= (float)(pG->Nx[0]);
00094
00095 for (i=is; i<=ie; i++) {
00096 pG->Phi[ks][js][i] -= total_Phi;
00097 }
00098 }
00099
00100
00101
00102
00103
00104
00105 void selfg_fft_2d(DomainS *pD)
00106 {
00107 GridS *pG = (pD->Grid);
00108 int i, is = pG->is, ie = pG->ie;
00109 int j, js = pG->js, je = pG->je;
00110 int ks = pG->ks;
00111 Real dx1sq=(pG->dx1*pG->dx1),dx2sq=(pG->dx2*pG->dx2);
00112 Real dkx,dky,pcoeff;
00113
00114
00115
00116 for (j=js-nghost; j<=je+nghost; j++){
00117 for (i=is-nghost; i<=ie+nghost; i++){
00118 pG->Phi_old[ks][j][i] = pG->Phi[ks][j][i];
00119 }
00120 }
00121
00122
00123
00124 for (j=js; j<=je; j++){
00125 for (i=is; i<=ie; i++){
00126 work[F2DI(i-is,j-js,pG->Nx[0],pG->Nx[1])][0] =
00127 four_pi_G*(pG->U[ks][j][i].d - grav_mean_rho);
00128 work[F2DI(i-is,j-js,pG->Nx[0],pG->Nx[1])][1] = 0.0;
00129 }
00130 }
00131
00132 ath_2d_fft(fplan2d, work);
00133
00134
00135
00136
00137
00138 dkx = 2.0*PI/(double)(pD->Nx[0]);
00139 dky = 2.0*PI/(double)(pD->Nx[1]);
00140
00141 if ((pG->Disp[1])==0 && (pG->Disp[0])==0) {
00142 work[F2DI(0,0,pG->Nx[0],pG->Nx[1])][0] = 0.0;
00143 work[F2DI(0,0,pG->Nx[0],pG->Nx[1])][1] = 0.0;
00144 } else {
00145 pcoeff = 1.0/(((2.0*cos((pG->Disp[0])*dkx)-2.0)/dx1sq) +
00146 ((2.0*cos((pG->Disp[1])*dky)-2.0)/dx2sq));
00147 work[F2DI(0,0,pG->Nx[0],pG->Nx[1])][0] *= pcoeff;
00148 work[F2DI(0,0,pG->Nx[0],pG->Nx[1])][1] *= pcoeff;
00149 }
00150
00151 for (j=js+1; j<=je; j++){
00152 pcoeff = 1.0/(((2.0*cos(( pG->Disp[0] )*dkx)-2.0)/dx1sq) +
00153 ((2.0*cos(( (j-js)+pG->Disp[1] )*dky)-2.0)/dx2sq));
00154 work[F2DI(0,j-js,pG->Nx[0],pG->Nx[1])][0] *= pcoeff;
00155 work[F2DI(0,j-js,pG->Nx[0],pG->Nx[1])][1] *= pcoeff;
00156 }
00157
00158 for (i=is+1; i<=ie; i++){
00159 for (j=js; j<=je; j++){
00160 pcoeff = 1.0/(((2.0*cos(( (i-is)+pG->Disp[0] )*dkx)-2.0)/dx1sq) +
00161 ((2.0*cos(( (j-js)+pG->Disp[1] )*dky)-2.0)/dx2sq));
00162 work[F2DI(i-is,j-js,pG->Nx[0],pG->Nx[1])][0] *= pcoeff;
00163 work[F2DI(i-is,j-js,pG->Nx[0],pG->Nx[1])][1] *= pcoeff;
00164 }
00165 }
00166
00167
00168
00169 ath_2d_fft(bplan2d, work);
00170
00171 for (j=js; j<=je; j++){
00172 for (i=is; i<=ie; i++){
00173 pG->Phi[ks][j][i] = work[F2DI(i-is,j-js,pG->Nx[0],pG->Nx[1])][0]/
00174 bplan2d->gcnt;
00175 }
00176 }
00177
00178 return;
00179 }
00180
00181
00182
00183
00184
00185
00186 void selfg_fft_3d(DomainS *pD)
00187 {
00188 GridS *pG = (pD->Grid);
00189 int i, is = pG->is, ie = pG->ie;
00190 int j, js = pG->js, je = pG->je;
00191 int k, ks = pG->ks, ke = pG->ke;
00192 Real dx1sq=(pG->dx1*pG->dx1),dx2sq=(pG->dx2*pG->dx2),dx3sq=(pG->dx3*pG->dx3);
00193 Real dkx,dky,dkz,pcoeff;
00194
00195
00196
00197 for (k=ks-nghost; k<=ke+nghost; k++){
00198 for (j=js-nghost; j<=je+nghost; j++){
00199 for (i=is-nghost; i<=ie+nghost; i++){
00200 pG->Phi_old[k][j][i] = pG->Phi[k][j][i];
00201 }
00202 }}
00203
00204
00205
00206 for (k=ks; k<=ke; k++){
00207 for (j=js; j<=je; j++){
00208 for (i=is; i<=ie; i++){
00209 work[F3DI(i-is,j-js,k-ks,pG->Nx[0],pG->Nx[1],pG->Nx[2])][0] =
00210 four_pi_G*(pG->U[k][j][i].d - grav_mean_rho);
00211 work[F3DI(i-is,j-js,k-ks,pG->Nx[0],pG->Nx[1],pG->Nx[2])][1] = 0.0;
00212 }
00213 }}
00214
00215 ath_3d_fft(fplan3d, work);
00216
00217
00218
00219
00220
00221 dkx = 2.0*PI/(double)(pD->Nx[0]);
00222 dky = 2.0*PI/(double)(pD->Nx[1]);
00223 dkz = 2.0*PI/(double)(pD->Nx[2]);
00224
00225 if ((pG->Disp[2])==0 && (pG->Disp[1])==0 && (pG->Disp[0])==0) {
00226 work[F3DI(0,0,0,pG->Nx[0],pG->Nx[1],pG->Nx[2])][0] = 0.0;
00227 work[F3DI(0,0,0,pG->Nx[0],pG->Nx[1],pG->Nx[2])][1] = 0.0;
00228 } else {
00229 pcoeff = 1.0/(((2.0*cos((pG->Disp[0])*dkx)-2.0)/dx1sq) +
00230 ((2.0*cos((pG->Disp[1])*dky)-2.0)/dx2sq) +
00231 ((2.0*cos((pG->Disp[2])*dkz)-2.0)/dx3sq));
00232 work[F3DI(0,0,0,pG->Nx[0],pG->Nx[1],pG->Nx[2])][0] *= pcoeff;
00233 work[F3DI(0,0,0,pG->Nx[0],pG->Nx[1],pG->Nx[2])][1] *= pcoeff;
00234 }
00235
00236
00237 for (k=ks+1; k<=ke; k++){
00238 pcoeff = 1.0/(((2.0*cos(( pG->Disp[0] )*dkx)-2.0)/dx1sq) +
00239 ((2.0*cos(( pG->Disp[1] )*dky)-2.0)/dx2sq) +
00240 ((2.0*cos(( (k-ks)+pG->Disp[2] )*dkz)-2.0)/dx3sq));
00241 work[F3DI(0,0,k-ks,pG->Nx[0],pG->Nx[1],pG->Nx[2])][0] *= pcoeff;
00242 work[F3DI(0,0,k-ks,pG->Nx[0],pG->Nx[1],pG->Nx[2])][1] *= pcoeff;
00243 }
00244
00245 for (j=js+1; j<=je; j++){
00246 for (k=ks; k<=ke; k++){
00247 pcoeff = 1.0/(((2.0*cos(( pG->Disp[0] )*dkx)-2.0)/dx1sq) +
00248 ((2.0*cos(( (j-js)+pG->Disp[1] )*dky)-2.0)/dx2sq) +
00249 ((2.0*cos(( (k-ks)+pG->Disp[2] )*dkz)-2.0)/dx3sq));
00250 work[F3DI(0,j-js,k-ks,pG->Nx[0],pG->Nx[1],pG->Nx[2])][0] *= pcoeff;
00251 work[F3DI(0,j-js,k-ks,pG->Nx[0],pG->Nx[1],pG->Nx[2])][1] *= pcoeff;
00252 }
00253 }
00254
00255 for (i=is+1; i<=ie; i++){
00256 for (j=js; j<=je; j++){
00257 for (k=ks; k<=ke; k++){
00258 pcoeff = 1.0/(((2.0*cos(( (i-is)+pG->Disp[0] )*dkx)-2.0)/dx1sq) +
00259 ((2.0*cos(( (j-js)+pG->Disp[1] )*dky)-2.0)/dx2sq) +
00260 ((2.0*cos(( (k-ks)+pG->Disp[2] )*dkz)-2.0)/dx3sq));
00261 work[F3DI(i-is,j-js,k-ks,pG->Nx[0],pG->Nx[1],pG->Nx[2])][0] *= pcoeff;
00262 work[F3DI(i-is,j-js,k-ks,pG->Nx[0],pG->Nx[1],pG->Nx[2])][1] *= pcoeff;
00263 }
00264 }}
00265
00266
00267
00268
00269 ath_3d_fft(bplan3d, work);
00270
00271 for (k=ks; k<=ke; k++){
00272 for (j=js; j<=je; j++){
00273 for (i=is; i<=ie; i++){
00274 pG->Phi[k][j][i] =
00275 work[F3DI(i-is,j-js,k-ks,pG->Nx[0],pG->Nx[1],pG->Nx[2])][0]
00276 / bplan3d->gcnt;
00277 }
00278 }}
00279
00280 return;
00281 }
00282
00283
00284
00285
00286
00287
00288
00289 void selfg_fft_2d_init(MeshS *pM)
00290 {
00291 DomainS *pD;
00292 int nl,nd;
00293 for (nl=0; nl<(pM->NLevels); nl++){
00294 for (nd=0; nd<(pM->DomainsPerLevel[nl]); nd++){
00295 if (pM->Domain[nl][nd].Grid != NULL){
00296 pD = (DomainS*)&(pM->Domain[nl][nd]);
00297 fplan2d = ath_2d_fft_quick_plan(pD, NULL, ATH_FFT_FORWARD);
00298 bplan2d = ath_2d_fft_quick_plan(pD, NULL, ATH_FFT_BACKWARD);
00299 work = ath_2d_fft_malloc(fplan2d);
00300 }
00301 }
00302 }
00303 }
00304
00305
00306
00307
00308
00309
00310
00311 void selfg_fft_3d_init(MeshS *pM)
00312 {
00313 DomainS *pD;
00314 int nl,nd;
00315 for (nl=0; nl<(pM->NLevels); nl++){
00316 for (nd=0; nd<(pM->DomainsPerLevel[nl]); nd++){
00317 if (pM->Domain[nl][nd].Grid != NULL){
00318 pD = (DomainS*)&(pM->Domain[nl][nd]);
00319 fplan3d = ath_3d_fft_quick_plan(pD, NULL, ATH_FFT_FORWARD);
00320 bplan3d = ath_3d_fft_quick_plan(pD, NULL, ATH_FFT_BACKWARD);
00321 work = ath_3d_fft_malloc(fplan3d);
00322 }
00323 }
00324 }
00325 }
00326
00327 #endif