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