Athena: prob/jet.c Source File

Go to the documentation of this file.
00001 #include "copyright.h"
00002 /*============================================================================*/
00003 /*! \file jet.c
00004  *  \brief Sets up a jet introduced through L-x1 boundary (left edge) */
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  * PRIVATE FUNCTION PROTOTYPES:
00017  * jet_iib() - sets BCs on L-x1 (left edge) of grid.
00018  *============================================================================*/
00019 
00020 void jet_iib(GridS *pGrid);
00021 
00022 /* Make radius of jet and jet variables static global so they can be accessed
00023    by BC functions */
00024 static Real rjet,Bxjet;
00025 static Prim1DS Wjet;
00026 static Cons1DS Ujet;
00027 static Real x1_mid,x2_mid,x3_mid;
00028 
00029 /*----------------------------------------------------------------------------*/
00030 /* problem */
00031 
00032 void problem(DomainS *pDomain){
00033    GridS *pGrid=(pDomain->Grid);
00034    int i, is = pGrid->is, ie = pGrid->ie;
00035    int j, js = pGrid->js, je = pGrid->je;
00036    int k, ks = pGrid->ks, ke = pGrid->ke;
00037    int il,iu,jl,ju,kl,ku;
00038    Prim1DS W;
00039    Cons1DS U;
00040    Real x1_min, x1_max;
00041    Real x2_min, x2_max;
00042    Real x3_min, x3_max;
00043    Real Bx=0.0;
00044    Bxjet = 0.0;
00045 
00046 /* read parameters from input file */
00047 
00048    W.d  = par_getd("problem", "d");
00049    W.P  = par_getd("problem", "p");
00050    W.Vx = par_getd("problem", "vx");
00051    W.Vy = par_getd("problem", "vy");
00052    W.Vz = par_getd("problem", "vz");
00053 #ifdef MHD
00054    Bx   = par_getd("problem", "bx");
00055    W.By = par_getd("problem", "by");
00056    W.Bz = par_getd("problem", "bz");
00057 #endif
00058 
00059    Wjet.d  = par_getd("problem", "djet");
00060    Wjet.P  = par_getd("problem", "pjet");
00061    Wjet.Vx = par_getd("problem", "vxjet");
00062    Wjet.Vy = par_getd("problem", "vyjet");
00063    Wjet.Vz = par_getd("problem", "vzjet");
00064 #ifdef MHD
00065    Bxjet   = par_getd("problem", "bxjet");
00066    Wjet.By = par_getd("problem", "byjet");
00067    Wjet.Bz = par_getd("problem", "bzjet");
00068 #endif
00069 
00070    rjet = par_getd("problem", "rjet");
00071    
00072    U = Prim1D_to_Cons1D(&W, &Bx);
00073    Ujet = Prim1D_to_Cons1D(&Wjet, &Bxjet);
00074 
00075    x1_min = pDomain->RootMinX[0];
00076    x1_max = pDomain->RootMaxX[0];
00077    x2_min = pDomain->RootMinX[1];
00078    x2_max = pDomain->RootMaxX[1];
00079    x3_min = pDomain->RootMinX[2];
00080    x3_max = pDomain->RootMaxX[2];
00081 
00082    x1_mid = 0.5 * (x1_max + x1_min);
00083    x2_mid = 0.5 * (x2_max + x2_min);
00084    x3_mid = 0.5 * (x3_max + x3_min);
00085 
00086 /* initialize index bounds assuming problem in xy plane */
00087    iu = ie + nghost;
00088    il = is - nghost;
00089    ju = je + nghost;
00090    jl = js - nghost;
00091    if(pGrid->Nx[2] > 1){
00092       ku = pGrid->ke + nghost;
00093       kl = pGrid->ks - nghost;
00094    }
00095    else{
00096       ku = pGrid->ke;
00097       kl = pGrid->ks;
00098    }
00099 
00100 /* initialize conserved variables */
00101    
00102    for(k=kl; k<=ku; k++){
00103       for(j=jl; j<=ju; j++){
00104          for(i=il; i<=iu; i++){
00105             pGrid->U[k][j][i].d  = U.d;
00106             pGrid->U[k][j][i].M1 = U.Mx;
00107             pGrid->U[k][j][i].M2 = U.My;
00108             pGrid->U[k][j][i].M3 = U.Mz;
00109             pGrid->U[k][j][i].E  = U.E;
00110 #ifdef MHD
00111             pGrid->U[k][j][i].B1c = Bx;
00112             pGrid->U[k][j][i].B2c = U.By;
00113             pGrid->U[k][j][i].B3c = U.Bz;
00114             pGrid->B1i[k][j][i] = Bx;
00115             pGrid->B2i[k][j][i] = U.By;
00116             pGrid->B3i[k][j][i] = U.Bz;
00117 #endif
00118          }
00119       }
00120    }
00121 
00122 /* Set boundary value function pointers */
00123 
00124    if (pDomain->Disp[0] == 0) bvals_mhd_fun(pDomain,left_x1,jet_iib);
00125 
00126 }
00127 
00128 
00129 /*==============================================================================
00130  * PROBLEM USER FUNCTIONS:
00131  * problem_write_restart() - writes problem-specific user data to restart files
00132  * problem_read_restart()  - reads problem-specific user data from restart files
00133  * get_usr_expr()          - sets pointer to expression for special output data
00134  * get_usr_out_fun()       - returns a user defined output function pointer
00135  * get_usr_par_prop()      - returns a user defined particle selection function
00136  * Userwork_in_loop        - problem specific work IN     main loop
00137  * Userwork_after_loop     - problem specific work AFTER  main loop
00138  *----------------------------------------------------------------------------*/
00139 
00140 void problem_write_restart(MeshS *pM, FILE *fp)
00141 {
00142   return;
00143 }
00144 
00145 void problem_read_restart(MeshS *pM, FILE *fp)
00146 {
00147   int nl,nd;
00148 
00149   Wjet.d  = par_getd("problem", "djet");
00150   Wjet.P  = par_getd("problem", "pjet");
00151   Wjet.Vx = par_getd("problem", "vxjet");
00152   Wjet.Vy = par_getd("problem", "vyjet");
00153   Wjet.Vz = par_getd("problem", "vzjet");
00154 #ifdef MHD
00155   Bxjet   = par_getd("problem", "bxjet");
00156   Wjet.By = par_getd("problem", "byjet");
00157   Wjet.Bz = par_getd("problem", "bzjet");
00158 #endif
00159 
00160   rjet = par_getd("problem", "rjet");
00161   Ujet = Prim1D_to_Cons1D(&Wjet, &Bxjet);
00162 
00163   for (nl=0; nl<pM->NLevels; nl++){
00164   for (nd=0; nd<pM->DomainsPerLevel[nl]; nd++){
00165     if (pM->Domain[nl][nd].Grid != NULL) {
00166       if (pM->Domain[nl][nd].Disp[0] == 0) 
00167         bvals_mhd_fun(&(pM->Domain[nl][nd]),left_x1,jet_iib);
00168     }
00169   }}
00170 
00171   return;
00172 }
00173 
00174 ConsFun_t get_usr_expr(const char *expr)
00175 {
00176   return NULL;
00177 }
00178 
00179 VOutFun_t get_usr_out_fun(const char *name){
00180   return NULL;
00181 }
00182 
00183 void Userwork_in_loop(MeshS *pM)
00184 {
00185   return;
00186 }
00187 
00188 void Userwork_after_loop(MeshS *pM)
00189 {
00190   return;
00191 }
00192 
00193 
00194 /*===================== PRIVATE FUNCTIONS ====================================*/
00195 
00196 /******************************************************************************/
00197 /*! \fn void jet_iib(GridS *pGrid) 
00198  *  \brief Sets ghost zones to either outflow BC or
00199  * if cell is within jet radius, to jet values */
00200 /******************************************************************************/
00201 
00202 void jet_iib(GridS *pGrid){
00203   int i, is = pGrid->is;
00204   int j, js = pGrid->js, je = pGrid->je;
00205   int k, ks = pGrid->ks, ke = pGrid->ke;
00206   Real rad,x1,x2,x3;
00207 
00208   for(k=ks; k<=ke; k++){
00209     for(j=js; j<=je; j++){
00210       for(i=1; i<=nghost; i++){
00211         cc_pos(pGrid,(is-i),j,k,&x1,&x2,&x3);
00212         rad = sqrt(SQR(x2 - x2_mid) + SQR(x3 - x3_mid));
00213             
00214         if(rad <= rjet){
00215           pGrid->U[k][j][is-i].d  = Ujet.d;
00216           pGrid->U[k][j][is-i].M1 = Ujet.Mx;
00217           pGrid->U[k][j][is-i].M2 = Ujet.My;
00218           pGrid->U[k][j][is-i].M3 = Ujet.Mz;
00219           pGrid->U[k][j][is-i].E  = Ujet.E;
00220 #ifdef MHD
00221           pGrid->U[k][j][is-i].B1c = Bxjet;
00222           pGrid->U[k][j][is-i].B2c = Ujet.By;
00223           pGrid->U[k][j][is-i].B3c = Ujet.Bz;
00224           pGrid->B1i[k][j][is-i] = Bxjet;
00225           pGrid->B2i[k][j][is-i] = Ujet.By;
00226           pGrid->B3i[k][j][is-i] = Ujet.Bz;
00227 #endif
00228         } else{
00229           pGrid->U[k][j][is-i] = pGrid->U[k][j][is+(i-1)];
00230           pGrid->U[k][j][is-i].M1 = -pGrid->U[k][j][is-i].M1;
00231 #ifdef MHD
00232           pGrid->U[k][j][is-i].B2c = -pGrid->U[k][j][is-i].B2c;
00233           pGrid->U[k][j][is-i].B3c = -pGrid->U[k][j][is-i].B3c;
00234           pGrid->B1i[k][j][is-i] =  pGrid->B1i[k][j][is+i];
00235           pGrid->B2i[k][j][is-i] = -pGrid->B2i[k][j][is+(i-1)];
00236           pGrid->B3i[k][j][is-i] = -pGrid->B3i[k][j][is+(i-1)];
00237 #endif
00238         }
00239       }
00240     }
00241   }
00242 }