00001 #include "copyright.h"
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015 #include <math.h>
00016 #include <stdio.h>
00017 #include <stdlib.h>
00018 #include "defs.h"
00019 #include "athena.h"
00020 #include "globals.h"
00021 #include "prototypes.h"
00022
00023
00024
00025
00026 void problem(DomainS *pDomain)
00027 {
00028 GridS *pGrid = (pDomain->Grid);
00029 int i=0,j=0,k=0;
00030 int is,ie,js,je,ks,ke,n;
00031 int kmax,jmax;
00032 Real amp,njeans;
00033 int kdir;
00034 Real d0,p0,u0,v0,w0,b0;
00035 Real x1,x2,x3,sinkx,coskx,xl;
00036 Real lambda,omega,omega2,cs,va,kwave;
00037 #ifdef MHD
00038 Real beta;
00039 #endif
00040 is = pGrid->is; ie = pGrid->ie;
00041 js = pGrid->js; je = pGrid->je;
00042 ks = pGrid->ks; ke = pGrid->ke;
00043
00044 #ifndef SELF_GRAVITY
00045 ath_error("[jeans]: this test only works with SELF_GRAVITY configured");
00046 #endif
00047
00048
00049 amp = par_getd("problem","amp");
00050 #ifdef MHD
00051 beta = par_getd("problem","beta");
00052 #endif
00053
00054 njeans = par_getd("problem","njeans");
00055
00056
00057 kdir =par_geti("problem","kdir");
00058
00059 switch(kdir){
00060 case 1:
00061 lambda = pDomain->Nx[0]*pGrid->dx1;
00062 break;
00063 case 2:
00064 lambda = pDomain->Nx[1]*pGrid->dx2;
00065 break;
00066 case 3:
00067 lambda = pDomain->Nx[2]*pGrid->dx3;
00068 }
00069
00070 d0 = 1.0;
00071
00072 p0 = 1.0;
00073
00074 u0 = 0.0;
00075 v0 = 0.0;
00076 w0 = 0.0;
00077 #ifdef MHD
00078 b0 = sqrt(2.*p0/beta);
00079 va=b0/sqrt(d0);
00080 #else
00081 va = 0.0;
00082 b0 = 0.0;
00083 #endif
00084
00085 #ifdef SELF_GRAVITY
00086
00087 four_pi_G = (4.0*Gamma*p0)*(PI*PI*njeans*njeans)/(d0*d0*lambda*lambda);
00088 grav_mean_rho = d0;
00089 #endif
00090
00091 #ifndef ISOTHERMAL
00092 cs=sqrt(Gamma*p0/d0);
00093 #else
00094 cs=Iso_csound;
00095 #endif
00096 kwave=2.0*PI/lambda;
00097
00098 omega2=kwave*kwave*(cs*cs + va*va - four_pi_G*d0/(kwave*kwave));
00099 omega=sqrt(fabs(omega2));
00100
00101 printf("4piG=%e, lambda=%e, period=%e\n",four_pi_G, lambda, (2.0*PI/omega));
00102
00103 kmax= (pGrid->Nx[2]>1)? ke+1 : ke;
00104 jmax= (pGrid->Nx[1]>1)? je+1 : je;
00105
00106
00107 for (k=ks; k<=kmax; k++) {
00108 for (j=js; j<=jmax; j++) {
00109 for (i=is; i<=ie+1; i++) {
00110 cc_pos(pGrid,i,j,k,&x1,&x2,&x3);
00111 switch(kdir){
00112 case 1:
00113 xl=x1-0.5*pGrid->dx1;
00114 sinkx = sin(x1*kwave);
00115 coskx = cos(x1*kwave);
00116 break;
00117 case 2:
00118 xl=x2-0.5*pGrid->dx2;
00119 sinkx = sin(x2*kwave);
00120 coskx = cos(x2*kwave);
00121 break;
00122 case 3:
00123 xl=x3-0.5*pGrid->dx3;
00124 sinkx = sin(x3*kwave);
00125 coskx = cos(x3*kwave);
00126 }
00127 pGrid->U[k][j][i].d = d0*(1.0 + amp*sinkx);
00128 #ifndef ISOTHERMAL
00129 pGrid->U[k][j][i].E = (p0/Gamma_1)*(1.0+Gamma*amp*sinkx) +
00130 b0*b0*(0.5+amp*sinkx);
00131 #endif
00132 switch(kdir){
00133 case 1:
00134 pGrid->U[k][j][i].M1 = (omega2<0.0)? d0*(omega/kwave)*amp*coskx: 0.0;
00135 pGrid->U[k][j][i].M2 = 0.0;
00136 pGrid->U[k][j][i].M3 = 0.0;
00137 break;
00138 case 2:
00139 pGrid->U[k][j][i].M1 = 0.0;
00140 pGrid->U[k][j][i].M2 = (omega2<0.0)? d0*(omega/kwave)*amp*coskx: 0.0;
00141 pGrid->U[k][j][i].M3 = 0.0;
00142 break;
00143 case 3:
00144 pGrid->U[k][j][i].M1 = 0.0;
00145 pGrid->U[k][j][i].M2 = 0.0;
00146 pGrid->U[k][j][i].M3 = (omega2<0.0)? d0*(omega/kwave)*amp*coskx: 0.0;
00147 }
00148 #ifdef MHD
00149 switch(kdir){
00150 case 1:
00151 pGrid->B1i[k][j][i] = 0.0;
00152 pGrid->B2i[k][j][i] = b0*(1.0+amp*sin(kwave*xl));
00153 pGrid->B3i[k][j][i] = 0.0;
00154 pGrid->U[k][j][i].B1c = 0.0;
00155 pGrid->U[k][j][i].B2c = b0*(1.0+amp*sinkx);
00156 pGrid->U[k][j][i].B3c = 0.0;
00157 break;
00158 case 2:
00159 pGrid->B1i[k][j][i] = b0*(1.0+amp*sin(kwave*xl));
00160 pGrid->B2i[k][j][i] = 0.0;
00161 pGrid->B3i[k][j][i] = 0.0;
00162 pGrid->U[k][j][i].B1c = b0*(1.0+amp*sinkx);
00163 pGrid->U[k][j][i].B2c = 0.0;
00164 pGrid->U[k][j][i].B3c = 0.0;
00165 break;
00166 case 3:
00167 pGrid->B1i[k][j][i] = b0*(1.0+amp*sin(kwave*xl));
00168 pGrid->B2i[k][j][i] = 0.0;
00169 pGrid->B3i[k][j][i] = 0.0;
00170 pGrid->U[k][j][i].B1c = b0*(1.0+amp*sinkx);
00171 pGrid->U[k][j][i].B2c = 0.0;
00172 pGrid->U[k][j][i].B3c = 0.0;
00173 }
00174 #endif
00175 #if (NSCALARS > 0)
00176 for (n=0; n<NSCALARS; n++)
00177 pGrid->U[k][j][i].s[n] = 0.0;
00178 #endif
00179 }}}
00180
00181 return;
00182 }
00183
00184
00185
00186
00187
00188
00189
00190
00191
00192
00193
00194
00195 void problem_write_restart(MeshS *pM, FILE *fp)
00196 {
00197 return;
00198 }
00199
00200 void problem_read_restart(MeshS *pM, FILE *fp)
00201 {
00202 return;
00203 }
00204
00205 #if (NSCALARS > 0)
00206
00207
00208 static Real color(const GridS *pG, const int i, const int j, const int k)
00209 {
00210 return pG->U[k][j][i].s[0]/pG->U[k][j][i].d;
00211 }
00212 #endif
00213
00214 ConsFun_t get_usr_expr(const char *expr)
00215 {
00216 #if (NSCALARS > 0)
00217 if(strcmp(expr,"color")==0) return color;
00218 #endif
00219 return NULL;
00220 }
00221
00222
00223 VOutFun_t get_usr_out_fun(const char *name){
00224 return NULL;
00225 }
00226
00227
00228 void Userwork_in_loop(MeshS *pM)
00229 {
00230 }
00231
00232
00233 void Userwork_after_loop(MeshS *pM)
00234 {
00235 }