00001 #include "copyright.h"
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020 #include <math.h>
00021 #include <stdio.h>
00022 #include <stdlib.h>
00023 #include <string.h>
00024 #include "defs.h"
00025 #include "athena.h"
00026 #include "globals.h"
00027 #include "prototypes.h"
00028
00029
00030
00031
00032
00033
00034
00035
00036
00037 static Real r0,phi0,amp,rad,rho0,omega0,vz0;
00038 static Real grav_pot(const Real x1, const Real x2, const Real x3);
00039 static Real grav_acc(const Real x1, const Real x2, const Real x3);
00040 Real M2(const Real x1, const Real x2, const Real x3);
00041 Real A3(const Real x1, const Real x2, const Real x3);
00042 static ConsS ***RootSoln=NULL;
00043
00044
00045
00046
00047 void problem(DomainS *pDomain)
00048 {
00049 GridS *pG = pDomain->Grid;
00050 int i=0,j=0,k=0;
00051 int is,ie,js,je,ks,ke,nx1,nx2,nx3;
00052 int il,iu,jl,ju,kl,ku;
00053 Real x1,x2,x3,lsf,rsf;
00054
00055 is = pG->is; ie = pG->ie;
00056 js = pG->js; je = pG->je;
00057 ks = pG->ks; ke = pG->ke;
00058
00059 il = is-nghost; iu = ie+nghost;
00060 jl = js-nghost; ju = je+nghost;
00061 kl = ks-nghost*(ke>ks); ku = ke+nghost*(ke>ks);
00062
00063 nx1 = iu-il+1;
00064 nx2 = ju-jl+1;
00065 nx3 = ku-kl+1;
00066
00067 #ifndef MHD
00068 ath_error("[cylfieldloop]: This problem can only be run in MHD!\n");
00069 #endif
00070
00071 if ((ie-is)==0 || (je-js)==0) {
00072 ath_error("[cylfieldloop]: This problem can only be run in 2D or 3D\n");
00073 }
00074
00075
00076 if ((RootSoln = (ConsS***)calloc_3d_array(nx3,nx2,nx1,sizeof(ConsS))) == NULL)
00077 ath_error("[cylfieldloop]: Error allocating memory for solution\n");
00078
00079
00080 r0 = par_getd("problem","r0");
00081 phi0 = par_getd("problem","phi0");
00082 amp = par_getd("problem","amp");
00083 rad = par_getd("problem","rad");
00084 rho0 = par_getd("problem","rho0");
00085 omega0 = par_getd("problem","omega0");
00086 vz0 = par_getd("problem","vz0");
00087
00088
00089 for (k=kl; k<=ku; k++) {
00090 for (j=jl; j<=ju; j++) {
00091 for (i=il; i<=iu; i++) {
00092 memset(&(pG->U[k][j][i]),0.0,sizeof(ConsS));
00093
00094 pG->U[k][j][i].d = rho0;
00095 pG->U[k][j][i].M1 = 0.0;
00096 pG->U[k][j][i].M2 = avg1d(M2,pG,i,j,k);
00097 pG->U[k][j][i].M3 = pG->U[k][j][i].d*vz0;
00098 pG->B1i[k][j][i] = vecpot2b1i(NULL,A3,pG,i,j,k);
00099 pG->B2i[k][j][i] = vecpot2b2i(NULL,A3,pG,i,j,k);
00100 pG->B3i[k][j][i] = 0.0;
00101 }
00102 }
00103 }
00104
00105
00106 for (k=kl; k<=ku; k++) {
00107 for (j=jl; j<=ju; j++) {
00108 for (i=il; i<=iu; i++) {
00109 cc_pos(pG,i,j,k,&x1,&x2,&x3);
00110 lsf = (x1 - 0.5*pG->dx1)/x1;
00111 rsf = (x1 + 0.5*pG->dx1)/x1;
00112
00113 if (i < iu)
00114 pG->U[k][j][i].B1c = 0.5*(lsf*pG->B1i[k][j][i] + rsf*pG->B1i[k][j][i+1]);
00115 else
00116 pG->U[k][j][i].B1c = 0.5*(lsf*pG->B1i[k][j][i] + rsf*vecpot2b1i(NULL,A3,pG,i+1,j,k));
00117
00118 if (j < ju)
00119 pG->U[k][j][i].B2c = 0.5*(pG->B2i[k][j][i] + pG->B2i[k][j+1][i]);
00120 else
00121 pG->U[k][j][i].B2c = 0.5*(pG->B2i[k][j][i] + vecpot2b2i(NULL,A3,pG,i,j+1,k));
00122
00123 if (ke > ks)
00124 if (k < ku)
00125 pG->U[k][j][i].B3c = 0.5*(pG->B3i[k][j][i] + pG->B3i[k+1][j][i]);
00126 else
00127 pG->U[k][j][i].B3c = 0.0;
00128 else
00129 pG->U[k][j][i].B3c = pG->B3i[k][j][i];
00130
00131 #ifndef ISOTHERMAL
00132 pG->U[k][j][i].E = 1.0/Gamma_1
00133 + 0.5*(SQR(pG->U[k][j][i].B1c) + SQR(pG->U[k][j][i].B2c) + SQR(pG->U[k][j][i].B3c))
00134 + 0.5*(SQR(pG->U[k][j][i].M1) + SQR(pG->U[k][j][i].M2) + SQR(pG->U[k][j][i].M3))/pG->U[k][j][i].d;
00135 #endif
00136
00137
00138 RootSoln[k][j][i] = pG->U[k][j][i];
00139 }
00140 }
00141 }
00142
00143 printf("Initial Max divB = %1.10e\n", compute_div_b(pG));
00144
00145 StaticGravPot = grav_pot;
00146 x1GravAcc = grav_acc;
00147 bvals_mhd_fun(pDomain, left_x1, do_nothing_bc);
00148 bvals_mhd_fun(pDomain, right_x1, do_nothing_bc);
00149
00150 return;
00151 }
00152
00153
00154
00155
00156
00157
00158
00159
00160
00161
00162
00163 void problem_write_restart(MeshS *pM, FILE *fp)
00164 {
00165 return;
00166 }
00167
00168 void problem_read_restart(MeshS *pM, FILE *fp)
00169 {
00170 return;
00171 }
00172
00173 ConsFun_t get_usr_expr(const char *expr)
00174 {
00175 return NULL;
00176 }
00177
00178 VOutFun_t get_usr_out_fun(const char *name){
00179 return NULL;
00180 }
00181
00182 void Userwork_in_loop(MeshS *pM)
00183 {
00184 printf("Max divB = %1.10e\n", compute_div_b(pM->Domain[0][0].Grid));
00185 }
00186
00187 void Userwork_after_loop(MeshS *pM)
00188 {
00189 compute_l1_error("CylFieldLoop", pM, RootSoln, 0);
00190 }
00191
00192
00193
00194
00195
00196
00197 static Real grav_pot(const Real x1, const Real x2, const Real x3) {
00198 return 0.5*SQR(x1*omega0);
00199 }
00200
00201
00202
00203 static Real grav_acc(const Real x1, const Real x2, const Real x3) {
00204 return x1*SQR(omega0);
00205 }
00206
00207
00208
00209 Real M2(const Real x1, const Real x2, const Real x3) {
00210 return rho0*omega0*x1;
00211 }
00212
00213
00214
00215 Real A3(const Real x1, const Real x2, const Real x3) {
00216 Real X0,X,Y0,Y,dist;
00217
00218
00219 X = x1*cos(x2); Y = x1*sin(x2);
00220
00221
00222 X0 = r0*cos(phi0); Y0 = r0*sin(phi0);
00223
00224
00225 dist = sqrt(SQR(X-X0) + SQR(Y-Y0));
00226 if (dist < rad) {
00227 return amp*(rad - dist);
00228 }
00229
00230 return 0.0;
00231 }