00001 #include "copyright.h"
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013 #include <math.h>
00014 #include <stdio.h>
00015 #include <stdlib.h>
00016 #include "defs.h"
00017 #include "athena.h"
00018 #include "globals.h"
00019 #include "prototypes.h"
00020
00021 #ifndef MHD
00022 #error: The cshock1d test only works for mhd.
00023 #endif
00024
00025 #ifndef ISOTHERMAL
00026 #error: The cshock1d test only works for isothermal equation of state.
00027 #endif
00028
00029 #ifndef RESISTIVITY
00030 #error: The cshock1d test only works with RESISTIVITY.
00031 #endif
00032
00033
00034
00035
00036 Real RK4(Real D, Real A, Real M, Real theta, Real h);
00037
00038 Real Dprime(Real D, Real A, Real M, Real theta);
00039
00040
00041 static ConsS *RootSoln=NULL;
00042
00043
00044
00045
00046 void problem(DomainS *pDomain)
00047 {
00048 GridS *pGrid=(pDomain->Grid);
00049 ConsS *Soln;
00050 int i, is = pGrid->is, ie = pGrid->ie;
00051 int j, js = pGrid->js, je = pGrid->je;
00052 int k, ks = pGrid->ks, ke = pGrid->ke;
00053 int i0, nx1, N_L;
00054 Real x1min,x1max,Lx,x1,x2,x3,xs,xe,h,x01,x02;
00055 Real Mach,Alfv,theta,d0,v0,vA,B0,Bx0,By0;
00056 Real Ls,Ns,D01,D02,myD;
00057
00058 nx1 = (ie-is)+1 + 2*nghost;
00059
00060
00061
00062 if ((Soln = (ConsS*)calloc_1d_array(nx1,sizeof(ConsS)))==NULL)
00063 ath_error("[problem]: Error allocating memory for Soln\n");
00064 if (pDomain->Level == 0){
00065 if ((RootSoln = (ConsS*)calloc_1d_array(nx1,sizeof(ConsS)))
00066 == NULL) ath_error("[problem]: Error alloc memory for RootSoln\n");
00067 }
00068
00069
00070 Mach = par_getd("problem","Mach");
00071 Alfv = par_getd("problem","Alfv");
00072
00073 theta = par_getd("problem","theta");
00074 theta = theta * PI / 180.0;
00075
00076
00077 d0 = 1.0;
00078 v0 = Mach * Iso_csound;
00079
00080 vA = (Mach/Alfv) * Iso_csound;
00081 B0 = sqrt(SQR(vA)*d0);
00082
00083 Bx0 = B0 * cos(theta);
00084 By0 = B0 * sin(theta);
00085
00086
00087 x1min = pDomain->RootMinX[0];
00088 x1max = pDomain->RootMaxX[0];
00089 Lx = x1max - x1min;
00090
00091
00092
00093
00094 Q_AD = 1.0/vA;
00095 d_ind= 0.0;
00096
00097
00098 Ls = par_getd_def("problem","Ls",20.0);
00099
00100 Ns = par_getd_def("problem","Ns",5e3);
00101
00102 if (Ls > Lx)
00103 ath_error("Domain size in x1 is shorter than the C-Shosck thickness!\n");
00104
00105 xs = x1min + 0.5*(Lx - Ls);
00106 xe = xs + Ls;
00107 h = (xe - xs)/Ns;
00108
00109
00110
00111
00112
00113
00114 i = is;
00115 cc_pos(pGrid,i,js,ks,&x1,&x2,&x3);
00116 while (x1 < xs){
00117
00118 Soln[i].d = d0;
00119 Soln[i].M1 = d0 * v0;
00120 Soln[i].M2 = 0.0;
00121 Soln[i].M3 = 0.0;
00122 Soln[i].B1c= Bx0;
00123 Soln[i].B2c= By0;
00124 Soln[i].B3c= 0.0;
00125
00126 i += 1;
00127 cc_pos(pGrid,i,js,ks,&x1,&x2,&x3);
00128 }
00129
00130
00131 x01 = xs; x02 = x01+h;
00132 D01 = d0+1.0e-6;
00133 while (x02 <= xe)
00134 {
00135 D02 = RK4(D01, Alfv, Mach, theta, h);
00136 if ((x1 >= x01) && (x1 < x02))
00137 {
00138 myD = (D01*(x02-x1)+D02*(x1-x01))/h;
00139 Soln[i].d = myD;
00140 Soln[i].B1c= Bx0;
00141 Soln[i].B2c= sqrt(SQR(By0)+2.0*SQR(Alfv*B0)*
00142 (myD-1)*(1.0/myD-SQR(1.0/Mach)));
00143 Soln[i].B3c= 0.0;
00144 Soln[i].M1 = d0 * v0;
00145 Soln[i].M2 = myD* SQR(vA)/v0*cos(theta)*(Soln[i].B2c/B0-sin(theta));
00146 Soln[i].M3 = 0.0;
00147
00148 i += 1;
00149 cc_pos(pGrid,i,js,ks,&x1,&x2,&x3);
00150 }
00151 x01 = x02;
00152 x02 += h;
00153 D01 = D02;
00154 }
00155 i0 = i-1;
00156
00157
00158 while (i <= ie+1)
00159 {
00160 Soln[i] = Soln[i0];
00161 i += 1;
00162 }
00163
00164
00165 if (pDomain->Level == 0){
00166 for (i=is; i<=ie+1; i++){
00167 RootSoln[i].d = Soln[i].d;
00168 RootSoln[i].M1 = Soln[i].M1;
00169 RootSoln[i].M2 = Soln[i].M2;
00170 RootSoln[i].M3 = Soln[i].M3;
00171 RootSoln[i].B1c= Soln[i].B1c;
00172 RootSoln[i].B2c= Soln[i].B2c;
00173 RootSoln[i].B3c= Soln[i].B3c;
00174 }
00175 }
00176
00177
00178 ie = ie+1;
00179 if (pGrid->Nx[1] > 1) je = je+1;
00180 if (pGrid->Nx[2] > 1) ke = ke+1;
00181
00182 for (k=ks; k<=ke; k++) {
00183 for (j=js; j<=je; j++) {
00184 for (i=is; i<=ie; i++) {
00185
00186 pGrid->U[k][j][i].d = Soln[i].d;
00187 pGrid->U[k][j][i].M1 = d0 * v0;
00188 pGrid->U[k][j][i].M2 = Soln[i].M2;
00189 pGrid->U[k][j][i].M3 = 0.0;
00190
00191 pGrid->U[k][j][i].B1c = pGrid->B1i[k][j][i] = Bx0;
00192 pGrid->U[k][j][i].B2c = pGrid->B2i[k][j][i] = Soln[i].B2c;
00193 pGrid->U[k][j][i].B3c = pGrid->B3i[k][j][i] = 0.0;
00194 }}}
00195
00196 return;
00197 }
00198
00199
00200
00201
00202
00203
00204
00205
00206
00207
00208
00209
00210 void problem_write_restart(MeshS *pM, FILE *fp)
00211 {
00212 return;
00213 }
00214
00215 void problem_read_restart(MeshS *pM, FILE *fp)
00216 {
00217 return;
00218 }
00219
00220 ConsFun_t get_usr_expr(const char *expr)
00221 {
00222 return NULL;
00223 }
00224
00225 VOutFun_t get_usr_out_fun(const char *name){
00226 return NULL;
00227 }
00228
00229 #ifdef RESISTIVITY
00230 void get_eta_user(GridS *pG, int i, int j, int k,
00231 Real *eta_O, Real *eta_H, Real *eta_A)
00232 {
00233 *eta_O = 0.0;
00234 *eta_H = 0.0;
00235 *eta_A = 0.0;
00236
00237 return;
00238 }
00239 #endif
00240
00241 void Userwork_in_loop(MeshS *pM)
00242 {
00243 }
00244
00245
00246
00247
00248
00249
00250
00251 void Userwork_after_loop(MeshS *pM)
00252 {
00253 }
00254
00255
00256
00257
00258
00259
00260 Real RK4(Real D, Real A, Real M, Real theta, Real h)
00261 {
00262 Real k1, k2, k3, k4;
00263
00264 k1 = Dprime(D,A,M,theta);
00265 k2 = Dprime(D+0.5*h*k1,A,M,theta);
00266 k3 = Dprime(D+0.5*h*k2,A,M,theta);
00267 k4 = Dprime(D+h*k3,A,M,theta);
00268
00269 return D + h/6.0*(k1+2.0*k2+2.0*k3+k4);
00270
00271 }
00272
00273
00274
00275 Real Dprime(Real D, Real A, Real M, Real theta)
00276 {
00277 Real sintheta, costheta, sintheta2, costheta2;
00278 Real M21, b, b2;
00279
00280 sintheta = sin(theta);
00281 costheta = cos(theta);
00282 sintheta2 = SQR(sintheta);
00283 costheta2 = SQR(costheta);
00284 M21 = 1.0/SQR(M);
00285
00286 b2 = sintheta2 + 2*SQR(A)*(D-1.0)*(1.0/D-M21);
00287 b = sqrt(b2);
00288
00289 return b/A*(b-D*((b-sintheta)/SQR(A)*costheta2+sintheta))/(b2+costheta2)/(1/SQR(D)-M21);
00290 }
00291