/***********************************************************************************
Copyright 2006 Charles F. Gammie, Jonathan C. McKinney, Scott C. Noble,
Gabor Toth, and Luca Del Zanna
TWODSSM version 1.0 (released December 2, 2007)
This file is part of TWODSSM. TWODSSM is a program that solves
the equations of self-gravitating hydrodynamics in the shearing sheet,
with cooling. TWODSSM is configured to use a simple, exponential cooling
model, although it contains code for a more complicated (and perhaps
more realistic) cooling model based on a one-zone vertical model. The
complicated cooling model can be switched on using a flag in the file
decs.h. This version of TWODSSM uses linear interpolation to do the
orbital advection, rather than consistent transport.
TWODSSM is based on ZEUS,
** Stone, J.M., and Norman, M. L., 1992, ApJS, 80, 753
modified to use ``orbital advection,'' i.e. it splits off advection due to the
mean flow and solves that piece by interpolation. This results in greater
accuracy and larger timesteps. Orbital advection is discussed in
** Masset, F., 2000, A&AS, 141, 165.
and in the following papers, which you are morally obligated to cite in any
formal or informal scientific discussion that results from the use of any
part of TWODSSM:
[1] Gammie, C. F., 2001, ApJ, 553, 174
[2] Johnson, B. M., & Gammie, C.F., 2003, ApJ, 597, 131
Further, we strongly encourage you to obtain the latest version of
TWODSSM directly from our distribution website:
http://rainman.astro.uiuc.edu/codelib/
TWODSSM 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.
TWODSSM 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 TWODSSM; if not, write to the Free Software
Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
This is version 1.1:
- fixed a bug discovered independently by both Bryan Johnson and
Yan-Fei Jiang in sweepy.c
- updated the imaging routine so it no longer requires a colormap file.
- included a timer for zone cycles/sec
- udpated the makefile to incorporate modern compiler flags
***********************************************************************************/
#include "decs.h"
#include "defs.h"
#include
#define TIMER_NSTEP (512)
int main(int argc, char *argv[])
{
void init(), timestep(), stepvar();
void diag(int flag);
void potsolv(double(*S)[NY + 4], double(*pot)[NY + 4]);
double tnext, mass, mass0;
int i, j;
int nstep = 0;
time_t ti_t,tf_t ;
double zcps ;
/* this was used in Gammie 2001 */
#if SIMPLECOOL
if (argc < 2) {
fprintf(stderr,
"usage: twodssm tc0 Lx Ly tf pamp is_restart\n");
exit(0);
} else {
sscanf(argv[1], "%lf", &tc0);
sscanf(argv[2], "%lf", &Lx);
sscanf(argv[3], "%lf", &Ly);
sscanf(argv[4], "%lf", &tf);
sscanf(argv[5], "%lf", &pamp);
sscanf(argv[8], "%d", &is_restart);
}
cool_fac = 1.0;
#endif
/* this was used in Johnson & Gammie */
#if FULLCOOL
if (argc < 2) {
fprintf(stderr,
"usage: twodssm T0 Lx Ly tf pamp rho0cgs cool_fac is_restart\n");
exit(0);
} else {
sscanf(argv[1], "%lf", &T0);
sscanf(argv[2], "%lf", &Lx);
sscanf(argv[3], "%lf", &Ly);
sscanf(argv[4], "%lf", &tf);
sscanf(argv[5], "%lf", &pamp);
sscanf(argv[6], "%lf", &rho0cgs);
sscanf(argv[7], "%lf", &cool_fac);
sscanf(argv[8], "%d", &is_restart);
}
#endif
/* perform initializations */
init();
fprintf(stderr, "cooling time: %g\n", tc0);
tnext = t;
mass0 = 0.;
for (i = 0; i < NX; i++)
for (j = 0; j < NY; j++)
mass0 += S[i][j];
mass0 /= NX * NY;
/* find potential */
if (G > 0.)
potsolv(S, pot);
/* do initial diagnostics */
diag(0);
ti_t = time(NULL) ; /* start timer */
tf_t = ti_t ;
while (t < tf) {
/* find timestep */
timestep();
/* step variables forward in time */
stepvar();
/* perform diagnostics */
if (t > tnext) {
diag(1);
tnext += DTl;
}
nstep++;
mass = 0.;
for (i = 0; i < NX; i++)
for (j = 0; j < NY; j++)
mass += S[i][j];
mass /= NX * NY;
if(nstep%TIMER_NSTEP == 0) {
tf_t = time(NULL) ;
zcps = nstep*NX*NY/((double)(tf_t - ti_t)) ;
fprintf(stderr,"Zone cycles/second: %g\n", zcps) ;
}
}
/* do final diagnostics */
diag(2);
fprintf(stderr, "nstep: %d\n", nstep);
return (0);
}