Go to the documentation of this file.00001 #include "copyright.h"
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014 #include <math.h>
00015 #include <stdio.h>
00016 #include <stdlib.h>
00017 #include "defs.h"
00018 #include "athena.h"
00019 #include "globals.h"
00020 #include "prototypes.h"
00021
00022 #ifndef MHD
00023 #error : The rotor problem can only be run with MHD.
00024 #endif
00025 #ifdef ISOTHERMAL
00026 #error : The rotor problem can only be run with an ADIABATIC eos.
00027 #endif
00028
00029
00030
00031
00032
00033 void problem(DomainS *pDomain)
00034 {
00035 GridS *pGrid = pDomain->Grid;
00036 int i,j,k,is,ie,js,je,ks,ke;
00037 Real v0,p0,bx0,x1,x2,x3,rad,frac,r0,r1;
00038
00039
00040
00041 v0 = par_getd("problem","v0");
00042 p0 = par_getd("problem","p0");
00043 bx0 = par_getd("problem","bx0");
00044 r0 = par_getd("problem","r0");
00045 r1 = par_getd("problem","r1");
00046
00047
00048
00049
00050 is = pGrid->is; ie = pGrid->ie;
00051 js = pGrid->js; je = pGrid->je;
00052 ks = pGrid->ks; ke = pGrid->ke;
00053
00054 for (k=ks; k<=ke; k++) {
00055 for (j=js; j<=je; j++) {
00056 for (i=is; i<=ie; i++) {
00057 pGrid->U[k][j][i].d = 1.0;
00058 pGrid->U[k][j][i].M1 = 0.0;
00059 pGrid->U[k][j][i].M2 = 0.0;
00060 pGrid->U[k][j][i].M3 = 0.0;
00061 pGrid->B1i[k][j][i] = bx0;
00062 pGrid->B2i[k][j][i] = 0.0;
00063 pGrid->U[k][j][i].B1c = bx0;
00064 pGrid->U[k][j][i].B2c = 0.0;
00065 pGrid->U[k][j][i].B3c = 0.0;
00066
00067
00068
00069 cc_pos(pGrid,i,j,k,&x1,&x2,&x3);
00070 rad = sqrt(x1*x1 + x2*x2);
00071 if (rad <= r0) {
00072 pGrid->U[k][j][i].d = 10.0;
00073 pGrid->U[k][j][i].M1 = -100.0*v0*x2;
00074 pGrid->U[k][j][i].M2 = 100.0*v0*x1;
00075 } else {
00076
00077
00078
00079 if (rad <= r1) {
00080 frac = (0.115 - rad)/(0.015);
00081 pGrid->U[k][j][i].d = 1.0 + 9.0*frac;
00082 pGrid->U[k][j][i].M1 = -frac*100.0*v0*x2;
00083 pGrid->U[k][j][i].M2 = frac*100.0*v0*x1;
00084 }
00085 }
00086
00087 pGrid->U[k][j][i].E = p0/Gamma_1 + 0.5*bx0*bx0
00088 + 0.5*(SQR(pGrid->U[k][j][i].M1) + SQR(pGrid->U[k][j][i].M2))
00089 /pGrid->U[k][j][i].d;
00090 }
00091 }
00092 }
00093
00094 for (k=ks; k<=ke; k++) {
00095 for (j=js; j<=je; j++) {
00096 pGrid->B1i[k][j][ie+1] = bx0;
00097 }
00098 }
00099
00100 return;
00101 }
00102
00103
00104
00105
00106
00107
00108
00109
00110
00111
00112
00113
00114 void problem_write_restart(MeshS *pM, FILE *fp)
00115 {
00116 return;
00117 }
00118
00119 void problem_read_restart(MeshS *pM, FILE *fp)
00120 {
00121 return;
00122 }
00123
00124 ConsFun_t get_usr_expr(const char *expr)
00125 {
00126 return NULL;
00127 }
00128
00129 VOutFun_t get_usr_out_fun(const char *name){
00130 return NULL;
00131 }
00132
00133 void Userwork_in_loop(MeshS *pM)
00134 {
00135 }
00136
00137 void Userwork_after_loop(MeshS *pM)
00138 {
00139 }