00001 #include "../copyright.h" 00002 /*============================================================================*/ 00003 /*! \file get_eta.c 00004 * \brief Functions to calculate the diffusion coefficients for resistivity. 00005 * 00006 * We provide two standard prescriptions, divided into 3 cases: 00007 * 00008 * - A. Single ion prescription (eq. 21 of Balbus & Terquem 2001, ApJ): 00009 * E = eta_Ohm*J + (JxB)/(4pi*n_e*e/c) - [(JxB)xB]/(4pi*gamma*rho_i*rho) 00010 * 00011 * -# CASE 1: eta_Ohm is constant and n_e and rho_i are proportional to gas 00012 * density rho. The user needs to provide eta_Ohm and two proportional 00013 * coefficients Q_Hall and Q_AD, which are defined in global.h. 00014 * 00015 * -# CASE 2: The user provides his/her own function to evaluate the etas, 00016 * e.g., spatially variable diffussivity, or adopting some chemical model. 00017 * The (blank) function is given in the problem generator. 00018 * 00019 *- B. Generalized prescription (eq. 25-31 of Wardle 2007, ApSS): 00020 * E = eta_Ohm*J + eta_Hall*(JxB)/B + eta_AD*J_perp 00021 * where eta_Ohm, eta_Hall and eta_AD are determined by sigma_O, sigma_H and 00022 * sigma_P. 00023 * 00024 * -# CASE 3: The user provides his/her own function to calculate the sigmas, 00025 * e.g., from non-equilibrium chemistry to evaluate them self-consistently. 00026 * The (blank) function is given in the problem generator. 00027 * 00028 * 00029 * CONTAINS PUBLIC FUNCTIONS: 00030 *- get_eta() - main function call to get magnetic diffusivities 00031 *- eta_single_const() - constant diffusivities for single ion prescription 00032 *- eta_single_user() - user defined diffusivities for single ion prescription 00033 *- eta_general_user() - user defined diffusivities for general prescription */ 00034 /*============================================================================*/ 00035 00036 #include <math.h> 00037 #include <float.h> 00038 #include "../defs.h" 00039 #include "../athena.h" 00040 #include "../globals.h" 00041 #include "prototypes.h" 00042 #include "../prototypes.h" 00043 00044 #ifdef RESISTIVITY 00045 00046 #ifdef HYDRO 00047 #error : resistivity only works for MHD. 00048 #endif /* HYDRO */ 00049 00050 /*============================================================================== 00051 * PRIVATE FUNCTION PROTOTYPES: 00052 * convert_diffusion() - convert conductivities to diffusion coefficients 00053 *============================================================================*/ 00054 00055 void convert_diffusion(Real sigma_O, Real sigma_H, Real sigma_P, 00056 Real *eta_O, Real *eta_H, Real *eta_A ); 00057 00058 /*=========================== PUBLIC FUNCTIONS ===============================*/ 00059 /*----------------------------------------------------------------------------*/ 00060 /*! \fn void get_eta(GridS *pG) 00061 * \brief Get magnetic diffusivities 00062 */ 00063 void get_eta(GridS *pG) 00064 { 00065 int i, il, iu, is = pG->is, ie = pG->ie; 00066 int j, jl, ju, js = pG->js, je = pG->je; 00067 int k, kl, ku, ks = pG->ks, ke = pG->ke; 00068 00069 il = is - 2; 00070 iu = ie + 2; 00071 if (pG->Nx[1] > 1){ 00072 jl = js - 2; 00073 ju = je + 2; 00074 } else { 00075 jl = js; 00076 ju = je; 00077 } 00078 if (pG->Nx[2] > 1){ 00079 kl = ks - 2; 00080 ku = ke + 2; 00081 } else { 00082 kl = ks; 00083 ku = ke; 00084 } 00085 00086 for (k=kl; k<=ku; k++) { 00087 for (j=jl; j<=ju; j++) { 00088 for (i=il; i<=iu; i++) { 00089 00090 get_myeta(pG, i,j,k, &(pG->eta_Ohm[k][j][i]), 00091 &(pG->eta_Hall[k][j][i]), &(pG->eta_AD[k][j][i])); 00092 00093 }}} 00094 00095 return; 00096 } 00097 00098 /*----------------------------------------------------------------------------*/ 00099 /*! \fn void eta_single_const(GridS *pG, int i, int j, int k, 00100 * Real *eta_O, Real *eta_H, Real *eta_A) 00101 * \brief Single-ion prescription with constant diffusivities 00102 * 00103 * NEED global parameters: eta_Ohm, Q_Hall and Q_AD. 00104 * ONLY eta_Ohm has the dimension of diffusivity. 00105 * This correspond to CASE 1. 00106 */ 00107 00108 void eta_single_const(GridS *pG, int i, int j, int k, 00109 Real *eta_O, Real *eta_H, Real *eta_A) 00110 { 00111 Real Bsq, Bmag; 00112 00113 *eta_O = eta_Ohm; 00114 00115 if ((Q_Hall > 0.0) || (Q_AD > 0.0)) 00116 { 00117 Bsq = SQR(pG->U[k][j][i].B1c) + SQR(pG->U[k][j][i].B2c) 00118 + SQR(pG->U[k][j][i].B3c); 00119 Bmag = sqrt(Bsq); 00120 00121 *eta_H = Q_Hall * Bmag / pow(pG->U[k][j][i].d, d_ind); 00122 *eta_A = Q_AD * Bsq / pow(pG->U[k][j][i].d, 1.0+d_ind); 00123 } 00124 00125 return; 00126 } 00127 00128 /*----------------------------------------------------------------------------*/ 00129 /*! \fn void eta_single_user(GridS *pG, int i, int j, int k, 00130 * Real *eta_O, Real *eta_H, Real *eta_A) 00131 * \brief Single-ion prescription with user defined diffusivities 00132 * 00133 * This corresponds to CASE 2. 00134 */ 00135 00136 void eta_single_user(GridS *pG, int i, int j, int k, 00137 Real *eta_O, Real *eta_H, Real *eta_A) 00138 { 00139 00140 get_eta_user(pG, i,j,k, eta_O, eta_H, eta_A); 00141 00142 return; 00143 } 00144 00145 /*----------------------------------------------------------------------------*/ 00146 /*! \fn void eta_general_user(GridS *pG, int i, int j, int k, 00147 * Real *eta_O, Real *eta_H, Real *eta_A) 00148 * \brief Generalized prescription with user defined conductivities 00149 * 00150 * This corresponds to CASE 3. 00151 * Note the same get_eta_user() function is used to calculate conductivities 00152 * rather than diffusivities here. 00153 */ 00154 00155 void eta_general_user(GridS *pG, int i, int j, int k, 00156 Real *eta_O, Real *eta_H, Real *eta_A) 00157 { 00158 Real sigma_O, sigma_H, sigma_P; 00159 00160 get_eta_user(pG, i,j,k, &sigma_O, &sigma_H, &sigma_P); 00161 00162 convert_diffusion(sigma_O, sigma_H, sigma_P, eta_O, eta_H, eta_A); 00163 00164 return; 00165 } 00166 00167 00168 /*----------------------------------------------------------------------------*/ 00169 /*! \fn void convert_diffusion(Real sigma_O, Real sigma_H, Real sigma_P, 00170 * Real *eta_O, Real *eta_H, Real *eta_A ) 00171 * \brief Convert conductivities to diffusivities 00172 * 00173 * NOTE: c^2/4pi factor is NOT included, so (eta x sigma) is dimensionless. 00174 * The conductivity coefficients should have the right dimension. 00175 */ 00176 00177 void convert_diffusion(Real sigma_O, Real sigma_H, Real sigma_P, 00178 Real *eta_O, Real *eta_H, Real *eta_A ) 00179 { 00180 Real sigma_perp2; 00181 00182 sigma_perp2 = SQR(sigma_H) + SQR(sigma_P); 00183 00184 *eta_O = 1.0/sigma_O; 00185 *eta_H = sigma_H / sigma_perp2; 00186 *eta_A = sigma_P / sigma_perp2 - *eta_O; 00187 00188 return; 00189 } 00190 00191 #endif /* RESISTIVITY */
1.7.1