00001 #include "copyright.h"
00002
00003
00004
00005
00006
00007 #include <math.h>
00008 #include <stdio.h>
00009 #include <stdlib.h>
00010 #include "defs.h"
00011 #include "athena.h"
00012 #include "globals.h"
00013 #include "prototypes.h"
00014
00015
00016
00017
00018
00019 static Real d0, d1, p0, vx0, vy0, vz0, bx0, by0, bz0, radius;
00020
00021
00022
00023
00024
00025
00026
00027 void problem(DomainS *pDomain)
00028 {
00029 GridS *pGrid = (pDomain->Grid);
00030 int i,is,ie,j,js,je,k,ks,ke,nx1,nx2,nx3;
00031 int ki, kj, kk, count, err, ndim;
00032 int iwhich = 0;
00033 Real x1,x2,x3,x1min,x1max,x2min,x2max,x3min,x3max;
00034 Real x1len,x2len,x3len,dx1,dx2,dx3,x10,x20,x30,r2;
00035 Real huge=1.0e60,radius=0.0,sig0=0.0;
00036
00037 #ifdef MPI_PARALLEL
00038 int myid = myID_Comm_world;
00039 #else
00040 int myid = 0;
00041 #endif
00042
00043 d0 = par_getd("problem","d0" );
00044 d1 = par_getd("problem","d1" );
00045 #ifndef ISOTHERMAL
00046 p0 = par_getd("problem","p0" );
00047 #endif
00048 vx0 = par_getd("problem","vx0" );
00049 vy0 = par_getd("problem","vy0" );
00050 vz0 = par_getd("problem","vz0" );
00051 #ifdef MHD
00052 bx0 = par_getd("problem","bx0");
00053 by0 = par_getd("problem","by0");
00054 bz0 = par_getd("problem","bz0");
00055 #endif
00056 x10 = par_getd("problem","x10" );
00057 x20 = par_getd("problem","x20" );
00058 x30 = par_getd("problem","x30" );
00059 radius = par_getd("problem","radius");
00060 sig0 = par_getd("problem","sig0");
00061 iwhich = par_geti("problem","iwhich");
00062 #ifdef SELF_GRAVITY
00063 four_pi_G= par_getd("problem","four_pi_G");
00064 #endif
00065 if (myid == 0) {
00066 fprintf(stdout,"[collapse3d]: d0 = %13.5e\n",d0);
00067 fprintf(stdout,"[collapse3d]: d1 = %13.5e\n",d1);
00068 #ifndef ISOTHERMAL
00069 fprintf(stdout,"[collapse3d]: p0 = %13.5e\n",p0);
00070 #endif
00071 #ifdef MHD
00072 fprintf(stdout,"[collapse3d]: bx0 = %13.5e\n",bx0);
00073 fprintf(stdout,"[collapse3d]: by0 = %13.5e\n",by0);
00074 fprintf(stdout,"[collapse3d]: bz0 = %13.5e\n",bz0);
00075 #endif
00076 fprintf(stdout,"[collapse3d]: x10 = %13.5e\n",x10);
00077 fprintf(stdout,"[collapse3d]: x20 = %13.5e\n",x20);
00078 fprintf(stdout,"[collapse3d]: x30 = %13.5e\n",x30);
00079 fprintf(stdout,"[collapse3d]: radius = %13.5e\n",radius);
00080 fprintf(stdout,"[collapse3d]: sig0 = %13.5e\n",sig0);
00081 #ifdef SELF_GRAVITY
00082 fprintf(stdout,"[collapse3d]: four_pi_G = %13.5e\n",four_pi_G);
00083 #endif
00084 }
00085
00086 is = pGrid->is; ie = pGrid->ie;
00087 js = pGrid->js; je = pGrid->je;
00088 ks = pGrid->ks; ke = pGrid->ke;
00089 nx1 = (ie-is)+1;
00090 nx2 = (je-js)+1;
00091 nx3 = (ke-ks)+1;
00092 ndim = 1;
00093 if (is == ie) {
00094 ath_error("[collapse3d]: This problem can only be run with Nx1>1\n");
00095 }
00096 nx1 += 2*nghost;
00097 if (nx2 > 1) {
00098 ndim++;
00099 nx2 += 2*nghost;
00100 }
00101 if (nx3 > 1) {
00102 ndim++;
00103 nx3 += 2*nghost;
00104 }
00105 if (myid == 0) {
00106 fprintf(stdout,"[collapse3d]: ndim = %13d\n",ndim);
00107 fprintf(stdout,"[collapse3d]: nx1,2,3 = %5d %5d %5d\n",nx1,nx2,nx3);
00108 fprintf(stdout,"[collapse3d]: is,js,ks = %5d %5d %5d\n",is,js,ks);
00109 fprintf(stdout,"[collapse3d]: ie,je,ke = %5d %5d %5d\n",ie,je,ke);
00110 }
00111
00112
00113 cc_pos(pGrid,is,js,ks,&x1min,&x2min,&x3min);
00114 cc_pos(pGrid,ie,je,ke,&x1max,&x2max,&x3max);
00115 dx1 = pGrid->dx1;
00116 dx2 = pGrid->dx2;
00117 dx3 = pGrid->dx3;
00118 x1min = par_getd("grid","x1min");
00119 x1max = par_getd("grid","x1max");
00120 x2min = par_getd("grid","x2min");
00121 x2max = par_getd("grid","x2max");
00122 x3min = par_getd("grid","x3min");
00123 x3max = par_getd("grid","x3max");
00124 x1len = x1max-x1min;
00125 x2len = x2max-x2min;
00126 x3len = x3max-x3min;
00127 if (myid == 0) {
00128 fprintf(stdout,"[collapse3d]: x1min = %13.5e, x2min = %13.5e, x3min = %13.5e\n",x1min,x2min,x3min);
00129 fprintf(stdout,"[collapse3d]: x1max = %13.5e, x2max = %13.5e, x3max = %13.5e\n",x1max,x2max,x3max);
00130 fprintf(stdout,"[collapse3d]: dx1 = %13.5e, dx2 = %13.5e, dx3 = %13.5e\n",dx1,dx2,dx3);
00131 fprintf(stdout,"[collapse3d]: x1len = %13.5e, x2len = %13.5e, x3len = %13.5e\n",x1len,x2len,x3len);
00132 }
00133
00134 if (ndim == 3) {
00135 for (k=ks; k<=ke; k++) {
00136 for (j=js; j<=je; j++) {
00137 for (i=is; i<=ie; i++) {
00138
00139 cc_pos(pGrid,i,j,k,&x1,&x2,&x3);
00140 r2 = sqrt(SQR(x1-x10)+SQR(x2-x20)+SQR(x3-x30));
00141 switch (iwhich) {
00142 case 0: {
00143 pGrid->U[k][j][i].d = d0 +
00144 0.5*(d1-d0)*(1.0-tanh((r2-radius)/(sig0*radius)));
00145 break;
00146 }
00147 case 1: {
00148 pGrid->U[k][j][i].d = d1 +
00149 (0.75*d0/(PI*radius*radius*radius))*pow((1.0+SQR(r2/radius)),-2.5);
00150 break;
00151 }
00152 default: ath_error("[collapse3d]: invalid iwhich\n");
00153 }
00154 pGrid->U[k][j][i].M1 = 0.0;
00155 pGrid->U[k][j][i].M2 = 0.0;
00156 pGrid->U[k][j][i].M3 = 0.0;
00157 #ifdef MHD
00158 pGrid->B1i[k][j][i] = bx0;
00159 pGrid->B2i[k][j][i] = by0;
00160 pGrid->B3i[k][j][i] = bz0;
00161 #endif
00162 }
00163 }
00164 }
00165 if (myid == 0)
00166 fprintf(stdout,"[collapse3d]: 3D setup finished\n");
00167 }
00168
00169
00170 #ifdef MHD
00171 for (k=ks; k<=ke; k++) {
00172 for (j=js; j<=je; j++) {
00173 pGrid->B1i[k][j][ie+1] = bx0;
00174 }
00175 }
00176 for (k=ks; k<=ke; k++) {
00177 for (i=is; i<=ie; i++) {
00178 pGrid->B2i[k][je+1][i] = by0;
00179 }
00180 }
00181 if (ndim == 3) {
00182 for (j=js; j<=je; j++) {
00183 for (i=is; i<=ie; i++) {
00184 pGrid->B3i[ke+1][j][i] = bz0;
00185 }
00186 }
00187 }
00188 #endif
00189
00190
00191
00192 if (ndim == 3) {
00193
00194 for (k=ks; k<=ke; k++) {
00195 for (j=js; j<=je; j++) {
00196 for (i=is; i<=ie; i++) {
00197 #ifdef MHD
00198 pGrid->U[k][j][i].B1c= 0.5*(pGrid->B1i[k][j][i]+pGrid->B1i[k][j][i+1]);
00199 pGrid->U[k][j][i].B2c= 0.5*(pGrid->B2i[k][j][i]+pGrid->B2i[k][j+1][i]);
00200 pGrid->U[k][j][i].B3c= 0.5*(pGrid->B3i[k][j][i]+pGrid->B3i[k+1][j][i]);
00201 #endif
00202 #ifndef ISOTHERMAL
00203 pGrid->U[k][j][i].E = p0/Gamma_1
00204 #ifdef MHD
00205 + 0.5*(SQR(pGrid->U[k][j][i].B1c) + SQR(pGrid->U[k][j][i].B2c)
00206 + SQR(pGrid->U[k][j][i].B3c))
00207 #endif
00208 + 0.5*(SQR(pGrid->U[k][j][i].M1) + SQR(pGrid->U[k][j][i].M2)
00209 + SQR(pGrid->U[k][j][i].M3))/pGrid->U[k][j][i].d;
00210 #endif
00211 }
00212 }
00213 }
00214 }
00215
00216 return;
00217 }
00218
00219
00220
00221
00222
00223
00224
00225
00226
00227
00228 void problem_write_restart(MeshS *pM, FILE *fp)
00229 {
00230 return;
00231 }
00232
00233
00234
00235
00236 void problem_read_restart(MeshS *pM, FILE *fp)
00237 {
00238 return;
00239 }
00240
00241 ConsFun_t get_usr_expr(const char *expr)
00242 {
00243 return NULL;
00244 }
00245
00246 VOutFun_t get_usr_out_fun(const char *name){
00247 return NULL;
00248 }
00249
00250 void Userwork_in_loop(MeshS *pM)
00251 {
00252 return;
00253 }
00254
00255 void Userwork_after_loop(MeshS *pM)
00256 {
00257 return;
00258 }