/***********************************************************************************
Copyright 2008 Charles F. Gammie and Xiaoyue Guan
HAM version 1.0 (released Jan 15, 2008)
This file is part of HAM. HAM is a program that solves non-relativistic
hyperbolic partial differential equation in conservative form using
high-resolution shock-capturing techniques. This version of HAM has been
configured to solve the magnetohydrodynamic equations of motion in
axisymmetry to evolve a shearing box model.
You are morally obligated to cite the following paper in his/her
scientific literature that results from use of any part of HAM:
[1] Guan, X. \& Gammie, C.F. \ 2008, Astrophysical Journal, Supplement,
174, 145
The latest version of the HAM is available on our code distribution
website:
http://rainman.astro.uiuc.edu/codelib/
HAM is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation; either version 2 of the License, or
(at your option) any later version.
HAM is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with HAM; if not, write to the Free Software
Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
***********************************************************************************/
/*
* Function init
*
* Generate initial conditions for a 2D shearingbox.
*
*/
#include "decs.h"
#define MAX(a,b) ((a) > (b) ? (a) : (b))
void
init ()
{
int i, j;
void set_arrays ();
double X1, X2, R, rho0, amp;
double shear;
double X, Y, Xp, Yp, Aij, Aijp, Aipj, Aipjp, Afunc (double X, double Y);
double B0, lambda;
int smooth;
/* assume throughout that G = rho = L = 1 */
dt = 0.00001;
L1 = 2.;
L2 = 2.;
/* set up arrays,details */
dx1 = L1 / N1;
dx2 = L2 / N2;
dV = dx1 * dx2;
t = 0.;
set_arrays ();
/* more hardwired choices */
tf = 32. * M_PI;
DTd = 1. / 5.; /* dumping frequency */
DTi = tf / 400.; /* imaging frequency */
DTl = tf / 1000.; /* logfile frequency */
amp = 0.01; /* pertubation amplitude */
rho0 = 1.;
cour = 0.9;
shear = 0.5;
q = 1.5;
W = 1.;
cs = 1.;
csq = cs * cs;
B0 = sqrt (15.) / (32. * M_PI); /* initial B field strength */
lambda = 2.;
LOOP
{
X1 = (i + 0.5) * dx1 - 0.5 * L1;
X2 = (j + 0.5) * dx2 - 0.5 * L2;
p[i][j][U1] = 0. + amp * (ranc (0) - 0.5);
p[i][j][U2] = 0. / (2. * M_PI) + amp * (ranc (0) - 0.5);
p[i][j][U3] = q * W * X1;
p[i][j][B1] = 0.;
p[i][j][B2] = B0 * sin (2. * M_PI * X1 / lambda);
p[i][j][B3] = 0.;
p[i][j][RHO] = rho0;
}
/* enforce boundary conditions */
bounds (p);
}
double
Afunc (double X, double Y)
{
double r;
r = sqrt (X * X + Y * Y);
return (1.e-5 * MAX (0.3 - r, 0.));
}