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 #include <math.h>
00026 #include <stdio.h>
00027 #include <stdlib.h>
00028 #include <string.h>
00029 #include "defs.h"
00030 #include "athena.h"
00031 #include "globals.h"
00032 #include "prototypes.h"
00033
00034 #ifndef ADIABATIC
00035 #error The Bernoulli solver assumes an adiabatic eos.
00036 #endif
00037
00038 #ifdef MHD
00039 #error The Bernoulli solver assumes a hydrodynamic fluid.
00040 #endif
00041
00042 static Real grav, psi;
00043 static Real H, S, Phi;
00044 static Real sin_a, cos_a;
00045 static Real lambda, k_par;
00046 static Real E0;
00047
00048 static int root;
00049 static Real ***d0=NULL;
00050
00051
00052
00053
00054
00055
00056
00057
00058
00059 static int rtbis(double (*pfun)(double), const double x1, const double x2,
00060 const double xacc, const int imax, double *prt);
00061 static Real grav_pot(const Real x1, const Real x2, const Real x3);
00062 static double Bfunc(double rho);
00063 static Real expr_drho(const Grid *pG, const int i, const int j, const int k);
00064
00065
00066
00067
00068
00069 void problem(Grid *pGrid, Domain *pDomain)
00070 {
00071 int i, is = pGrid->is, ie = pGrid->ie, nx1;
00072 int j, js = pGrid->js, je = pGrid->je, nx2;
00073 int k, ks = pGrid->ks, ke = pGrid->ke, nx3;
00074 Real x1,x2,x3;
00075 Real den,v_par,pres;
00076 Real rho_p,rho_e,rho_s;
00077 Real dt;
00078 double rho;
00079 double angle;
00080 nx1 = (ie-is)+1 + 2*nghost;
00081 nx2 = (je-js)+1 + 2*nghost;
00082 nx3 = (ke-ks)+1 + 2*nghost;
00083 if ((d0 = (Real***)calloc_3d_array(nx3,nx2,nx1,sizeof(Real))) == NULL)
00084 ath_error("[pgflow]: Error allocating memory\n");
00085
00086
00087 angle = par_getd("problem","angle");
00088
00089
00090 if(pGrid->Nx2 <= 1) angle = 0.0;
00091 if(pGrid->Nx1 <= 1) angle = 90.0;
00092
00093
00094 if (angle == 0.0) {
00095 sin_a = 0.0;
00096 cos_a = 1.0;
00097 lambda = pGrid->Nx1*pGrid->dx1;
00098 }
00099 else if (angle == 90.0) {
00100 sin_a = 1.0;
00101 cos_a = 0.0;
00102 lambda = pGrid->Nx2*pGrid->dx2;
00103 }
00104 else {
00105
00106
00107
00108
00109 if((pGrid->Nx1*pGrid->dx1) == (pGrid->Nx2*pGrid->dx2)){
00110 cos_a = sin_a = sqrt(0.5);
00111 }
00112 else{
00113 angle = atan((double)(pGrid->Nx1*pGrid->dx1)/(pGrid->Nx2*pGrid->dx2));
00114 sin_a = sin(angle);
00115 cos_a = cos(angle);
00116 }
00117
00118 if (cos_a >= sin_a) {
00119 lambda = pGrid->Nx1*pGrid->dx1*cos_a;
00120 } else {
00121 lambda = pGrid->Nx2*pGrid->dx2*sin_a;
00122 }
00123 }
00124
00125
00126 k_par = 2.0*PI/lambda;
00127
00128 E0 = 0.0;
00129 grav = par_getd("problem","grav");
00130 root = par_geti("problem","root");
00131
00132 den = par_getd("problem","den");
00133 pres = par_getd("problem","pres");
00134 v_par = par_getd("problem","v_par");
00135
00136
00137 Phi = den*v_par;
00138 S = pres/pow((double)den,(double)Gamma);
00139
00140 H = 0.5*v_par*v_par + Gamma*pres/(Gamma_1*den);
00141
00142 rho_e = pow((double)(Phi*Phi/(Gamma*S)),(double)(1.0/(Gamma+1.0)));
00143
00144 for (k=ks; k<=ke; k++) {
00145 for (j=js; j<=je; j++) {
00146 for (i=is; i<=ie; i++) {
00147 cc_pos(pGrid,i,j,k,&x1,&x2,&x3);
00148
00149
00150 psi = -grav*sin((double)k_par*(x1*cos_a + x2*sin_a))/k_par;
00151
00152 if(H <= psi)
00153 ath_error("[problem]: H < Psi -- No solution exists\n");
00154
00155 if(Bfunc((double)rho_e) < 0.0)
00156 ath_error("[problem]: Bfunc(rho_e) < 0.0 -- No solution exists\n");
00157
00158 if(root){
00159
00160 rho_s = pow((double)(Gamma_1*(H-psi)/(Gamma*S)),(double)(1.0/Gamma_1));
00161
00162 if(rtbis(Bfunc,rho_e,rho_s,1.0e-12*rho_e,100,&rho)){
00163 exit(1);
00164 }
00165 } else {
00166
00167 rho_p = fabs((double)Phi)/sqrt((double)(2.0*(H-psi)));
00168
00169 if(rtbis(Bfunc,rho_p,rho_e,1.0e-12*rho_e,100,&rho)){
00170 exit(1);
00171 }
00172 }
00173
00174 d0[k][j][i] = rho;
00175 pGrid->U[k][j][i].d = rho;
00176 pGrid->U[k][j][i].M1 = Phi*cos_a;
00177 pGrid->U[k][j][i].M2 = Phi*sin_a;
00178 pGrid->U[k][j][i].M3 = 0.0;
00179 pGrid->U[k][j][i].E = S*pow(rho,(double)Gamma)/Gamma_1
00180 + 0.5*Phi*Phi/rho;
00181
00182 E0 += pGrid->U[k][j][i].E + rho*psi;
00183 }
00184 }}
00185
00186
00187 E0 /= (Real)((ie - is + 1)*(je - js + 1)*(ke - ks + 1));
00188
00189
00190 StaticGravPot = grav_pot;
00191
00192 return;
00193 }
00194
00195
00196
00197
00198
00199
00200
00201
00202
00203
00204
00205
00206 void problem_write_restart(Grid *pG, Domain *pD, FILE *fp)
00207 {
00208 return;
00209 }
00210
00211 void problem_read_restart(Grid *pG, Domain *pD, FILE *fp)
00212 {
00213 return;
00214 }
00215
00216 Gasfun_t get_usr_expr(const char *expr)
00217 {
00218 if(strcmp(expr,"drho")==0) return expr_drho;
00219 return NULL;
00220 }
00221
00222 VGFunout_t get_usr_out_fun(const char *name){
00223 return NULL;
00224 }
00225
00226 void Userwork_in_loop(Grid *pGrid, Domain *pDomain)
00227 {
00228 }
00229
00230 void Userwork_after_loop(Grid *pGrid, Domain *pDomain)
00231 {
00232 }
00233
00234
00235
00236
00237
00238
00239
00240
00241
00242
00243
00244
00245
00246
00247
00248
00249
00250
00251 int rtbis(double (*pfun)(double), const double x1, const double x2,
00252 const double xacc, const int imax, double *prt)
00253 {
00254 int i;
00255 double xn, xm, xp, dx;
00256 double fn = (*pfun)(x1), fm, fp = (*pfun)(x2);
00257
00258 if(fn < 0.0 && fp > 0.0){
00259 xn = x1;
00260 xp = x2;
00261 }
00262 else if(fn > 0.0 && fp < 0.0){
00263 xn = x2;
00264 xp = x1;
00265
00266 fm = fn;
00267 fn = fp;
00268 fp = fm;
00269 }
00270 else if(fn == 0.0){ *prt = x1; return 0; }
00271 else if(fp == 0.0){ *prt = x2; return 0; }
00272 else{
00273 ath_perr(-1,"[rtbis]: Root must be bracketed for bisection\n");
00274 ath_perr(-1,"[rtbis]: x1=%g, x2=%g, F(x1)=%g, F(x2)=%g\n",
00275 x1,x2,fn,fp);
00276 return 1;
00277 }
00278
00279 dx = xp - xn;
00280
00281 for(i=0; i<imax; i++){
00282 dx *= 0.5;
00283 xm = xn + dx;
00284 fm = (*pfun)(xm);
00285
00286 if(fm < 0.0){ xn = xm; fn = fm; }
00287 else { xp = xm; fp = fm; }
00288
00289 if(fabs(dx) < xacc || fm == 0.0){
00290 *prt = fp < -fn ? xp : xn;
00291 return 0;
00292 }
00293 }
00294
00295 ath_perr(-1,"[rtbis]: Too many bisections\n");
00296
00297 *prt = fp < -fn ? xp : xn;
00298
00299 return 1;
00300 }
00301
00302
00303
00304
00305
00306
00307 static Real grav_pot(const Real x1, const Real x2, const Real x3)
00308 {
00309 return -grav*sin((double)k_par*(x1*cos_a + x2*sin_a))/k_par;
00310 }
00311
00312
00313
00314
00315
00316
00317 static double Bfunc(double rho){
00318 return (H - psi - 0.5*Phi*Phi/(rho*rho)
00319 - (Gamma*S/Gamma_1)*pow((double)rho,(double)Gamma_1));
00320 }
00321
00322
00323
00324
00325
00326
00327
00328 static Real expr_drho(const Grid *pG, const int i, const int j, const int k)
00329 {
00330 return (pG->U[k][j][i].d - d0[k][j][i]);
00331 }