Athena: prob/field_loop.c Source File

Go to the documentation of this file.
00001 #include "copyright.h"
00002 /*============================================================================*/
00003 /*! \file field_loop.c
00004  *  \brief Problem generator for advection of a field loop test. 
00005  *
00006  * PURPOSE: Problem generator for advection of a field loop test.  Can only
00007  *   be run in 2D or 3D.  Input parameters are:
00008  *   -  problem/rad   = radius of field loop
00009  *   -  problem/amp   = amplitude of vector potential (and therefore B)
00010  *   -  problem/vflow = flow velocity
00011  *   -  problem/drat  = density ratio in loop.  Enables density advection and
00012  *                      thermal conduction tests.
00013  *   The flow is automatically set to run along the diagonal. 
00014  *
00015  *   Various test cases are possible:
00016  *   - (iprob=1): field loop in x1-x2 plane (cylinder in 3D)
00017  *   - (iprob=2): field loop in x2-x3 plane (cylinder in 3D)
00018  *   - (iprob=3): field loop in x3-x1 plane (cylinder in 3D) 
00019  *   - (iprob=4): rotated cylindrical field loop in 3D.
00020  *   - (iprob=5): spherical field loop in rotated plane
00021  *
00022  *   A sphere of passive scalar can be added to test advection of scalars.
00023  *
00024  *   The temperature in the loop can be changed using drat to test conduction.
00025  *
00026  * REFERENCE: T. Gardiner & J.M. Stone, "An unsplit Godunov method for ideal MHD
00027  *   via constrined transport", JCP, 205, 509 (2005)                          */
00028 /*============================================================================*/
00029 
00030 #include <math.h>
00031 #include <stdio.h>
00032 #include <stdlib.h>
00033 #include <string.h>
00034 #include "defs.h"
00035 #include "athena.h"
00036 #include "globals.h"
00037 #include "prototypes.h"
00038 
00039 /*----------------------------------------------------------------------------*/
00040 /* problem:   */
00041 
00042 void problem(DomainS *pDomain)
00043 {
00044   GridS *pGrid = pDomain->Grid;
00045   int i=0,j=0,k=0;
00046   int is,ie,js,je,ks,ke,nx1,nx2,nx3,iprob;
00047   Real x1c,x2c,x3c,x1f,x2f,x3f;       /* cell- and face-centered coordinates */
00048   Real x1size,x2size,x3size,lambda=0.0,ang_2=0.0,sin_a2=0.0,cos_a2=1.0,x,y;
00049   Real rad,amp,vflow,drat,diag;
00050   Real ***az,***ay,***ax;
00051 #ifdef MHD
00052   int ku;
00053 #endif
00054 #if (NSCALARS > 0)
00055   int n;
00056 #endif
00057 
00058   is = pGrid->is; ie = pGrid->ie;
00059   js = pGrid->js; je = pGrid->je;
00060   ks = pGrid->ks; ke = pGrid->ke;
00061   nx1 = (ie-is)+1 + 2*nghost;
00062   nx2 = (je-js)+1 + 2*nghost;
00063   nx3 = (ke-ks)+1 + 2*nghost;
00064 
00065   if (((je-js) == 0)) {
00066     ath_error("[field_loop]: This problem can only be run in 2D or 3D\n");
00067   }
00068 
00069   if ((ay = (Real***)calloc_3d_array(nx3, nx2, nx1, sizeof(Real))) == NULL) {
00070     ath_error("[field_loop]: Error allocating memory for vector pot\n");
00071   }
00072   if ((az = (Real***)calloc_3d_array(nx3, nx2, nx1, sizeof(Real))) == NULL) {
00073     ath_error("[field_loop]: Error allocating memory for vector pot\n");
00074   }
00075   if ((ax = (Real***)calloc_3d_array(nx3, nx2, nx1, sizeof(Real))) == NULL) {
00076     ath_error("[field_loop]: Error allocating memory for vector pot\n");
00077   }
00078 
00079 /* Read initial conditions, diffusion coefficients (if needed) */
00080 
00081   rad = par_getd("problem","rad");
00082   amp = par_getd("problem","amp");
00083   vflow = par_getd("problem","vflow");
00084   drat = par_getd_def("problem","drat",1.0);
00085   iprob = par_getd("problem","iprob");
00086 #ifdef OHMIC
00087   eta_R = par_getd("problem","eta");
00088 #endif
00089 #ifdef ISOTROPIC_CONDUCTION
00090   kappa_T = par_getd("problem","kappa");
00091 #endif
00092 #ifdef ANISOTROPIC_CONDUCTION
00093   chi_C = par_getd("problem","chi");
00094 #endif
00095 
00096 /* For (iprob=4) -- rotated cylinder in 3D -- set up rotation angle and
00097  * wavelength of cylinder */
00098 
00099   if(iprob == 4){
00100     x1size = pDomain->RootMaxX[0] - pDomain->RootMinX[0];
00101     x3size = pDomain->RootMaxX[2] - pDomain->RootMinX[2];
00102 
00103 /* We put 1 wavelength in each direction.  Hence the wavelength
00104  *     lambda = x1size*cos_a;
00105  *     AND   lambda = x3size*sin_a;
00106  *     are both satisfied. */
00107 
00108     if(x1size == x3size){
00109       ang_2 = PI/4.0;
00110       cos_a2 = sin_a2 = sqrt(0.5);
00111     }
00112     else{
00113       ang_2 = atan(x1size/x3size);
00114       sin_a2 = sin(ang_2);
00115       cos_a2 = cos(ang_2);
00116     }
00117 /* Use the larger angle to determine the wavelength */
00118     if (cos_a2 >= sin_a2) {
00119       lambda = x1size*cos_a2;
00120     } else {
00121       lambda = x3size*sin_a2;
00122     }
00123   }
00124 
00125 /* Use vector potential to initialize field loop */
00126 
00127   for (k=ks; k<=ke+1; k++) {
00128   for (j=js; j<=je+1; j++) {
00129   for (i=is; i<=ie+1; i++) {
00130     cc_pos(pGrid,i,j,k,&x1c,&x2c,&x3c);
00131     x1f = x1c - 0.5*pGrid->dx1;
00132     x2f = x2c - 0.5*pGrid->dx2;
00133     x3f = x3c - 0.5*pGrid->dx3;
00134      
00135 /* (iprob=1): field loop in x1-x2 plane (cylinder in 3D) */
00136 
00137     if(iprob==1) {  
00138       ax[k][j][i] = 0.0;
00139       ay[k][j][i] = 0.0;
00140       if ((x1f*x1f + x2f*x2f) < rad*rad) {
00141         az[k][j][i] = amp*(rad - sqrt(x1f*x1f + x2f*x2f));
00142       } else {
00143         az[k][j][i] = 0.0;
00144       }
00145     }
00146 
00147 /* (iprob=2): field loop in x2-x3 plane (cylinder in 3D) */
00148 
00149     if(iprob==2) {  
00150       if ((x2f*x2f + x3f*x3f) < rad*rad) {
00151         ax[k][j][i] = amp*(rad - sqrt(x2f*x2f + x3f*x3f));
00152       } else {
00153         ax[k][j][i] = 0.0;
00154       }
00155       ay[k][j][i] = 0.0;
00156       az[k][j][i] = 0.0;
00157     }
00158 
00159 /* (iprob=3): field loop in x3-x1 plane (cylinder in 3D) */
00160 
00161     if(iprob==3) {  
00162       if ((x1f*x1f + x3f*x3f) < rad*rad) {
00163         ay[k][j][i] = amp*(rad - sqrt(x1f*x1f + x3f*x3f));
00164       } else {
00165         ay[k][j][i] = 0.0;
00166       }
00167       ax[k][j][i] = 0.0;
00168       az[k][j][i] = 0.0;
00169     }
00170 
00171 /* (iprob=4): rotated cylindrical field loop in 3D.  Similar to iprob=1
00172  * with a rotation about the x2-axis.  Define coordinate systems (x1,x2,x3)
00173  * and (x,y,z) with the following transformation rules:
00174  *    x =  x1*cos(ang_2) + x3*sin(ang_2)
00175  *    y =  x2
00176  *    z = -x1*sin(ang_2) + x3*cos(ang_2)
00177  * This inverts to:
00178  *    x1  = x*cos(ang_2) - z*sin(ang_2)
00179  *    x2  = y
00180  *    x3  = x*sin(ang_2) + z*cos(ang_2)
00181  */
00182 
00183     if(iprob==4) {
00184       x = x1c*cos_a2 + x3f*sin_a2;
00185       y = x2f;
00186 /* shift x back to the domain -0.5*lambda <= x <= 0.5*lambda */
00187       while(x >  0.5*lambda) x -= lambda;
00188       while(x < -0.5*lambda) x += lambda;
00189       if ((x*x + y*y) < rad*rad) {
00190         ax[k][j][i] = amp*(rad - sqrt(x*x + y*y))*(-sin_a2);
00191       } else {
00192         ax[k][j][i] = 0.0;
00193       }
00194 
00195       ay[k][j][i] = 0.0;
00196 
00197       x = x1f*cos_a2 + x3c*sin_a2;
00198       y = x2f;
00199 /* shift x back to the domain -0.5*lambda <= x <= 0.5*lambda */
00200       while(x >  0.5*lambda) x -= lambda;
00201       while(x < -0.5*lambda) x += lambda;
00202       if ((x*x + y*y) < rad*rad) {
00203         az[k][j][i] = amp*(rad - sqrt(x*x + y*y))*(cos_a2);
00204       } else {
00205         az[k][j][i] = 0.0;
00206       }
00207     }
00208 
00209 /* (iprob=5): spherical field loop in rotated plane */
00210 
00211     if(iprob==5) { 
00212       ax[k][j][i] = 0.0;
00213       if ((x1f*x1f + x2c*x2c + x3f*x3f) < rad*rad) {
00214         ay[k][j][i] = amp*(rad - sqrt(x1f*x1f + x2c*x2c + x3f*x3f));
00215       } else {
00216         ay[k][j][i] = 0.0;
00217       }
00218       if ((x1f*x1f + x2f*x2f + x3c*x3c) < rad*rad) {
00219         az[k][j][i] = amp*(rad - sqrt(x1f*x1f + x2f*x2f + x3c*x3c));
00220       } else {
00221         az[k][j][i] = 0.0;
00222       }
00223     }
00224 
00225   }}}
00226 
00227 /* Initialize density and momenta.  If drat != 1, then density and temperature
00228  * will be different inside loop than background values
00229  */
00230 
00231   x1size = pDomain->RootMaxX[0] - pDomain->RootMinX[0];
00232   x2size = pDomain->RootMaxX[1] - pDomain->RootMinX[1];
00233   x3size = pDomain->RootMaxX[2] - pDomain->RootMinX[2];
00234   diag = sqrt(x1size*x1size + x2size*x2size + x3size*x3size);
00235   for (k=ks; k<=ke; k++) {
00236   for (j=js; j<=je; j++) {
00237   for (i=is; i<=ie; i++) {
00238      pGrid->U[k][j][i].d = 1.0;
00239      pGrid->U[k][j][i].M1 = pGrid->U[k][j][i].d*vflow*x1size/diag;
00240      pGrid->U[k][j][i].M2 = pGrid->U[k][j][i].d*vflow*x2size/diag;
00241      pGrid->U[k][j][i].M3 = pGrid->U[k][j][i].d*vflow*x3size/diag;
00242 
00243      cc_pos(pGrid,i,j,k,&x1c,&x2c,&x3c);
00244      if ((x1c*x1c + x2c*x2c + x3c*x3c) < rad*rad) {
00245        pGrid->U[k][j][i].d = drat;
00246        pGrid->U[k][j][i].M1 = pGrid->U[k][j][i].d*vflow*x1size/diag;
00247        pGrid->U[k][j][i].M2 = pGrid->U[k][j][i].d*vflow*x2size/diag;
00248        pGrid->U[k][j][i].M3 = pGrid->U[k][j][i].d*vflow*x3size/diag;
00249      }
00250 #if (NSCALARS > 0)
00251      for (n=0; n<NSCALARS; n++) pGrid->U[k][j][i].s[n] = 0.0;
00252      if ((x1c*x1c + x2c*x2c + x3c*x3c) < rad*rad) {
00253        for (n=0; n<NSCALARS; n++)  pGrid->U[k][j][i].s[n] = 1.0;
00254      }
00255 #endif
00256   }}}
00257 
00258 /* boundary conditions on interface B */
00259 
00260 #ifdef MHD
00261   for (k=ks; k<=ke; k++) {
00262   for (j=js; j<=je; j++) {
00263   for (i=is; i<=ie+1; i++) {
00264     pGrid->B1i[k][j][i] = (az[k][j+1][i] - az[k][j][i])/pGrid->dx2 -
00265                           (ay[k+1][j][i] - ay[k][j][i])/pGrid->dx3;
00266   }}}
00267   for (k=ks; k<=ke; k++) {
00268   for (j=js; j<=je+1; j++) {
00269   for (i=is; i<=ie; i++) {
00270     pGrid->B2i[k][j][i] = (ax[k+1][j][i] - ax[k][j][i])/pGrid->dx3 -
00271                           (az[k][j][i+1] - az[k][j][i])/pGrid->dx1;
00272   }}}
00273   if (ke > ks) {
00274     ku = ke+1;
00275   } else {
00276     ku = ke;
00277   }
00278   for (k=ks; k<=ku; k++) {
00279   for (j=js; j<=je; j++) {
00280   for (i=is; i<=ie; i++) {
00281     pGrid->B3i[k][j][i] = (ay[k][j][i+1] - ay[k][j][i])/pGrid->dx1 -
00282                           (ax[k][j+1][i] - ax[k][j][i])/pGrid->dx2;
00283   }}}
00284 #endif
00285 
00286 /* initialize total energy and cell-centered B */
00287 
00288 #if defined MHD || !defined ISOTHERMAL
00289   for (k=ks; k<=ke; k++) {
00290     for (j=js; j<=je; j++) {
00291       for (i=is; i<=ie; i++) {
00292 #ifdef MHD
00293         pGrid->U[k][j][i].B1c = 0.5*(pGrid->B1i[k][j][i  ] + 
00294                                      pGrid->B1i[k][j][i+1]);
00295         pGrid->U[k][j][i].B2c = 0.5*(pGrid->B2i[k][j  ][i] +
00296                                      pGrid->B2i[k][j+1][i]);
00297         if (ke > ks)
00298           pGrid->U[k][j][i].B3c = 0.5*(pGrid->B3i[k  ][j][i] + 
00299                                        pGrid->B3i[k+1][j][i]);
00300         else
00301           pGrid->U[k][j][i].B3c = pGrid->B3i[k][j][i];
00302 #endif
00303 
00304 #ifndef ISOTHERMAL
00305         pGrid->U[k][j][i].E = 1.0/Gamma_1
00306 #ifdef MHD
00307           + 0.5*(SQR(pGrid->U[k][j][i].B1c) + SQR(pGrid->U[k][j][i].B2c)
00308                + SQR(pGrid->U[k][j][i].B3c))
00309 #endif
00310           + 0.5*(SQR(pGrid->U[k][j][i].M1) + SQR(pGrid->U[k][j][i].M2)
00311                + SQR(pGrid->U[k][j][i].M3))/pGrid->U[k][j][i].d;
00312 #endif /* ISOTHERMAL */
00313       }
00314     }
00315   }
00316 #endif
00317 
00318   free_3d_array((void***)az);
00319   free_3d_array((void***)ay);
00320   free_3d_array((void***)ax);
00321 }
00322 
00323 /*==============================================================================
00324  * PROBLEM USER FUNCTIONS:
00325  * problem_write_restart() - writes problem-specific user data to restart files
00326  * problem_read_restart()  - reads problem-specific user data from restart files
00327  * get_usr_expr()          - sets pointer to expression for special output data
00328  * get_usr_out_fun()       - returns a user defined output function pointer
00329  * get_usr_par_prop()      - returns a user defined particle selection function
00330  * Userwork_in_loop        - problem specific work IN     main loop
00331  * Userwork_after_loop     - problem specific work AFTER  main loop
00332  * current() - computes x3-component of current
00333  * Bp2()     - computes magnetic pressure (Bx2 + By2)
00334  * color()   - returns first passively advected scalar s[0]
00335  * Temperature() - returns temperature for conduction tests
00336  *----------------------------------------------------------------------------*/
00337 
00338 void problem_write_restart(MeshS *pM, FILE *fp)
00339 {
00340   return;
00341 }
00342 
00343 void problem_read_restart(MeshS *pM, FILE *fp)
00344 {
00345   return;
00346 }
00347 
00348 #ifdef MHD
00349 /*! \fn static Real current(const GridS *pG, const int i, const int j, const 
00350  *                         int k)
00351  *  \brief computes x3-component of current
00352  */
00353 static Real current(const GridS *pG, const int i, const int j, const int k)
00354 {
00355   return ((pG->B2i[k][j][i]-pG->B2i[k][j][i-1])/pG->dx1 - 
00356           (pG->B1i[k][j][i]-pG->B1i[k][j-1][i])/pG->dx2);
00357 }
00358 
00359 /*! \fn static Real Bp2(const GridS *pG, const int i, const int j, const int k)
00360  *  \brief computes magnetic pressure (Bx2 + By2) */
00361 static Real Bp2(const GridS *pG, const int i, const int j, const int k)
00362 {
00363   return (pG->U[k][j][i].B1c*pG->U[k][j][i].B1c + 
00364           pG->U[k][j][i].B2c*pG->U[k][j][i].B2c);
00365 }
00366 
00367 /*! \fn static Real divB(const GridS *pG, const int i, const int j, const int k)
00368  *  \brief  calculates div(B) */
00369 static Real divB(const GridS *pG, const int i, const int j, const int k)
00370 {
00371   Real qa;
00372   if (pG->Nx[2] > 1) {
00373     qa = (pG->B1i[k][j][i+1]-pG->B1i[k][j][i])/pG->dx1 + 
00374          (pG->B2i[k][j+1][i]-pG->B2i[k][j][i])/pG->dx2 +
00375          (pG->B3i[k+1][j][i]-pG->B3i[k][j][i])/pG->dx3;
00376   } else {
00377     qa = (pG->B1i[k][j][i+1]-pG->B1i[k][j][i])/pG->dx1 + 
00378          (pG->B2i[k][j+1][i]-pG->B2i[k][j][i])/pG->dx2;
00379   }
00380   return qa;
00381 }
00382 #endif
00383 
00384 #if (NSCALARS > 0)
00385 /*! \fn static Real color(const GridS *pG, const int i, const int j,const int k)
00386  *  \brief returns first passively advected scalar s[0]
00387 static Real color(const GridS *pG, const int i, const int j, const int k)
00388 {
00389   return pG->U[k][j][i].s[0]/pG->U[k][j][i].d;
00390 }
00391 #endif
00392 
00393 #ifndef BAROTROPIC
00394 /*! \fn static Real Temperature(const GridS *pG, const int i, const int j, 
00395  *                              const int k)
00396  *  \brief returns temperature for conduction tests */
00397 static Real Temperature(const GridS *pG, const int i, const int j, const int k)
00398 {
00399   Real Temp;
00400   Temp = pG->U[k][j][i].E - (0.5/pG->U[k][j][i].d)*(
00401     SQR(pG->U[k][j][i].M1) + SQR(pG->U[k][j][i].M2) + SQR(pG->U[k][j][i].M3));
00402 #ifdef MHD
00403   Temp -= 0.5*(SQR(pG->U[k][j][i].B1c) + SQR(pG->U[k][j][i].B2c) 
00404             + SQR(pG->U[k][j][i].B3c));
00405 #endif
00406   Temp *= Gamma_1/pG->U[k][j][i].d;
00407   return Temp;
00408 }
00409 #endif
00410 
00411 ConsFun_t get_usr_expr(const char *expr)
00412 {
00413 #ifdef MHD
00414   if(strcmp(expr,"J3")==0) return current;
00415   else if(strcmp(expr,"Bp2")==0) return Bp2;
00416   else if(strcmp(expr,"divB")==0) return divB;
00417 #endif
00418 #if (NSCALARS > 0)
00419   if(strcmp(expr,"color")==0) return color;
00420 #endif
00421   if(strcmp(expr,"Temperature")==0) return Temperature;
00422   return NULL;
00423 }
00424 
00425 VOutFun_t get_usr_out_fun(const char *name){
00426   return NULL;
00427 }
00428 
00429 void Userwork_in_loop(MeshS *pM)
00430 {
00431 }
00432 
00433 void Userwork_after_loop(MeshS *pM)
00434 {
00435 }