pluto.h

Go to the documentation of this file.

/* ///////////////////////////////////////////////////////////////////// */
/* ///////////////////////////////////////////////////////////////////// */
#ifndef PLUTO_H
#define PLUTO_H

#define PLUTO_VERSION  "4.0"

#include <stdio.h>
#include <stdarg.h>
#include <string.h>
#include <math.h>
#include <stdlib.h>
#include <time.h>

#define MAX(a,b)  ( (a) >= (b) ? (a) : (b) ) 

#define MIN(a,b)  ( (a) <= (b) ? (a) : (b) ) 
      
#define ABS_MIN(a,b)  (fabs(a) < fabs(b) ? (a):(b)) 
                         
#define DSIGN(x)      ( (x) >= 0.0 ? (1.0) : (-1.0))

#define MINMOD(a,b)  ((a)*(b) > 0.0 ? (fabs(a) < fabs(b) ? (a):(b)):0.0)
#define VAN_LEER(a,b) ((a)*(b) > 0.0 ? 2.0*(a)*(b)/((a)+(b)):0.0)
#define MC(a,b) (MINMOD(0.5*((a)+(b)), 2.0*MINMOD((a),(b))))
#define SWAP_VAR(x) SwapEndian(&x, sizeof(x));

#define YES      1
#define NO       0
#define DEFAULT -1

         /*  ----  GEOMETRY LABELS ( > 0) ----  */

#define CARTESIAN    1
#define AXISYM       1   /*  should be the same as CARTESIAN in order 
                             to avoid complications                    */
#define CYLINDRICAL  2
#define POLAR        3
#define SPHERICAL    4

#define NEW_GEOM     1

#define UNIFORM_GRID             1
#define STRETCHED_GRID           2
#define LOGARITHMIC_INC_GRID     3
#define LOGARITHMIC_DEC_GRID     4

    /*  ----  SET LABELS FOR EQUATION OF STATE (EOS)  ----  */

#define IDEAL         1
#define TAUB          3
#define BAROTROPIC    4
#define ISOTHERMAL    5

         /*  ----  TIME_STEPPING LABELS  ----  */

#define EULER                      1
#define HANCOCK                    2
#define CHARACTERISTIC_TRACING     3
#define RK2                        5
#define RK3                        6
#define RK_MIDPOINT                7

#define EXPLICIT             1 /* -- just a number different from 0 !!!  -- */
#define SUPER_TIME_STEPPING  2 /* -- just a number different from EXPLICIT -- */ 
#define RK_CHEBYSHEV         4  

     /*  ----   Output labels:  ---- */

#define DBL_OUTPUT   1
#define FLT_OUTPUT   2
#define VTK_OUTPUT   3
#define DBL_H5_OUTPUT   4
#define FLT_H5_OUTPUT   5
#define TAB_OUTPUT   6
#define PPM_OUTPUT   7
#define PNG_OUTPUT   8

#define VTK_VECTOR  5  /* -- any number but NOT 1  -- */

#define MAX_OUTPUT_TYPES 11 

   /*  ----  SET LABELS FOR COOLING  ----  */

#define POWER_LAW    3
#define MINEq        4
#define SNEq         5
#define TABULATED    6
#define H2_COOL      7

  /*  ----  SET LABELS FOR PHYSICS MODULE  ----  */

#define HD     1
#define RHD    2
#define MHD    3
#define RMHD   5

    /*  ----  SET LABELS FOR DIV.B REMOVAL  ----  
        If you move them to the MHD header, 
        definitions.h (which is included before)
        cannot correctly use them                */
        
#define NONE                   0
#define EIGHT_WAVES            1
#define DIV_CLEANING           2
#define CONSTRAINED_TRANSPORT  3

   /*  ----  SET LABELS FOR BODY_FORCE  ----
      Please do not change them since they are
      used in bitwise operations                */
   
#define VECTOR     4   /* corresponds to  100 in binary  */
#define POTENTIAL  8   /* corresponds to 1000 in binary  */

   /*  ----  SET LABELS FOR BOUNDARIES  ----  */

#define OUTFLOW          1  /* any number except 0 !! */
#define REFLECTIVE       2 
#define AXISYMMETRIC     3
#define EQTSYMMETRIC     4
#define PERIODIC         5
#define SHEARING         6
#define USERDEF          7

#define X1_BEG           101
#define X1_END           102
#define X2_BEG           103
#define X2_END           104
#define X3_BEG           105
#define X3_END           106

    /* ---- LABELS FOR IMAGE SLICING ---- */

#define X12_PLANE       3  
#define X13_PLANE       5
#define X23_PLANE       6

#define FLAG_MINMOD      1  
#define FLAG_FLAT        2  
#define FLAG_HLL         4  
#define FLAG_ENTROPY     8  
#define FLAG_SPLIT_CELL  16 
#define FLAG_INTERNAL_BOUNDARY   32  
#define FLAG_BIT7        64
#define FLAG_BIT8       128  

#define PRS_FAIL  1
#define ENG_FAIL  2
#define RHO_FAIL  4

#define IDIR    0     /*   This sequence (0,1,2) should */
#define JDIR    1     /*   never be changed             */
#define KDIR    2     /*                                */
#define ALL_DIR -1

/* -- location of a variable inside the cell -- */

#define CENTER  0
#define X1FACE  1
#define X2FACE  2
#define X3FACE  3

#define CELL_CENTER    50  /* really needed ? */
#define FACE_CENTER    51
#define EDGE_CENTER    52

        /*  ----  SET LABELS FOR INTERPOLATION  ----   */

#define FLAT              1
#define LINEAR            2
#define CENO3             3
#define PARABOLIC         4
#define LINEAR_MULTID     5
#define MP5               6
#define LimO3             7
#define WENO3             8

#define WENO3_FD             103
#define WENO5_FD             105
#define WENOZ_FD             106
#define WENO7_FD             107
#define MP5_FD               125
#define PPM_FD               140
#define LIMO3_FD             300

#define ONED   1
#define MULTID 3

/* ---- limiter labels ---- */

#define FLAT_LIM          1
#define MINMOD_LIM        2 
#define VANALBADA_LIM     3
#define UMIST_LIM         4
#define VANLEER_LIM       5
#define MC_LIM            6
#define FOURTH_ORDER_LIM  7

/*

#define FLAT_LIMIT(dv, dp, dm) dv = 0.0;

#define MINMOD_LIMIT(dv, dp, dm)  dv = ((dm)*(dp) > 0.0 ? (fabs(dm) < fabs(dp) ? (dm):(dp)):0.0)

#define VANALBADA_LIMIT(dv, dp, dm)\
  if (dp*dm > 0.0) { \
    double _dpp= dp*dp, _dmm = dm_dm; \
    dv = dp*(_dmm + 1.e-18) + dm*(_dpp + 1.e-18))/(_dpp + _dmm + 1.e-18); \
  }else dv = 0.0;

#define VANLEER_LIMIT(dv, dp, dm)  if (dp*dm > 0.0) dv = 2.0*dp*dm/(dp+dm);\
                                 else             dv = 0.0;

#define UMIST_LIMIT(dv,dp,dm)\
  if (dp*dm > 0.0){ \
  double _ddp = 0.25*(dp + 3.0*dm), _ddm = 0.25*(dm + 3.0*dp); \
  double _d2  = 2.0*(fabs(dp) < fabs(dm) ? dp:dm);  \
  _d2 = (fabs(_d2) < fabs(_ddp) ? _d2:_ddp); \
  dv  = (fabs(_d2) < fabs(_ddm) ? _d2:_ddm);} else dv = 0.0;

#define MC_LIMIT(dv, dp, dm, dc)  \
  if (dp*dm > 0.0){  \
    double _d2 = 2.0*(fabs(dp) < fabs(dm)) ? dp:dm); \
    dv = 2.0*(fabs(dvp) < fabs(dvm) ? dvp:dvm); \
  }else dv = 0.0;
*/
 

/* ###############################################################

                 PROBLEM-DEPENDENT DECLARATIONS

   ############################################################### */

#include "definitions.h"

#ifdef PARALLEL
 #include <al.h>
#endif

/* -- include mpi.h for parallel Chombo, in order to use
      MPI_Abort function in the QUIT_PLUTO macro -- */
 
#ifdef CH_MPI
 #include <mpi.h>
#endif


/* ###################################################################

        SWITCHES USED FOR DEBUGGING PURPOSES ONLY

   ################################################################### */
 
/* -- CHECK_EIGENVECTORS: used in eigenv.c in HD/, MHD/, RHD/
      to check orthogonality and the correctness through 
      the relation the A = L*\Lambda*R  -- */

#define CHECK_EIGENVECTORS      NO

/* -- CHECK_ROE_MATRIX: used in HD/roe.c, MHD/roe.c to 
      check that the characteristic decomp. reproduces
      Roe matrix -- */

#define CHECK_ROE_MATRIX        NO

/* -- CHECK_CONSERVATIVE_VAR: used in RHD/mappers.c to 
      check that conservative vars are physical -- */

#define CHECK_CONSERVATIVE_VAR  NO

/* -- CHECK_DIVB_CONDITION: used in MHD/CT/ct.c 
      to check if div.B = 0 -- */

#ifndef CHECK_DIVB_CONDITION
 #define CHECK_DIVB_CONDITION   NO
#endif

/* ###############################################################

                 Define precision (useless now)
   
   ############################################################### */

typedef double real;

/* ###############################################################
    
     EXIT MACRO. For Chombo it is defined elsewhere.

   ############################################################### */

#ifdef PARALLEL
/* #define QUIT_PLUTO(e_code)   {MPI_Finalize(); exit(e_code);}  */
 #define QUIT_PLUTO(e_code)   \
        {MPI_Abort(MPI_COMM_WORLD, e_code);MPI_Finalize(); exit(e_code);}
#elif (defined CH_MPI)
 #define QUIT_PLUTO(e_code)   \
        {MPI_Abort(MPI_COMM_WORLD, e_code); exit(e_code);}
#else
 #define QUIT_PLUTO(e_code)   exit(e_code);
#endif
        
/* ##################################################################

           MACROS TO EXPAND DIMENSION-DEPENDENT LINES

   #################################################################  */

#if COMPONENTS == 1
 #define EXPAND(a,b,c) a
 #define SELECT(a,b,c) a
#endif

#if COMPONENTS == 2
 #define EXPAND(a,b,c) a b
 #define SELECT(a,b,c) b
#endif

#if COMPONENTS == 3
 #define EXPAND(a,b,c) a b c
 #define SELECT(a,b,c) c
#endif

#if DIMENSIONS == 1
 #define D_EXPAND(a,b,c)  a
 #define D_SELECT(a,b,c)  a
 #define IOFFSET 1   /* ** so we know when to loop in boundary zones ** */
 #define JOFFSET 0
 #define KOFFSET 0
#endif

#if DIMENSIONS == 2
 #define D_EXPAND(a,b,c) a b
 #define D_SELECT(a,b,c) b
 #define IOFFSET 1   /* ** so we know when to loop in boundary zones ** */
 #define JOFFSET 1
 #define KOFFSET 0
#endif

#if DIMENSIONS == 3
 #define D_EXPAND(a,b,c) a b c
 #define D_SELECT(a,b,c) c
 #define IOFFSET 1   /* ** so we know when to loop in boundary zones ** */
 #define JOFFSET 1
 #define KOFFSET 1
#endif

#if WARNING_MESSAGES == YES
 #define WARNING(a)  a
#else
 #define WARNING(a)
#endif

/* ################################################################# 

                       SPATIAL LOOP MACROS        

   ################################################################# */

#define VAR_LOOP(n)   for ((n) = NVAR; (n)--;    )
#define DIM_LOOP(d)   for ((d) = 0; (d) < DIMENSIONS; (d)++)

#define IBEG_LOOP(i)  for ((i) = IBEG; (i)--;    )
#define JBEG_LOOP(j)  for ((j) = JBEG; (j)--;    )
#define KBEG_LOOP(k)  for ((k) = KBEG; (k)--;    )

#define IEND_LOOP(i)  for ((i) = IEND + 1; (i) < NX1_TOT; (i)++)
#define JEND_LOOP(j)  for ((j) = JEND + 1; (j) < NX2_TOT; (j)++)
#define KEND_LOOP(k)  for ((k) = KEND + 1; (k) < NX3_TOT; (k)++)

#define IDOM_LOOP(i)  for ((i) = IBEG; (i) <= IEND; (i)++)
#define JDOM_LOOP(j)  for ((j) = JBEG; (j) <= JEND; (j)++)
#define KDOM_LOOP(k)  for ((k) = KBEG; (k) <= KEND; (k)++)

#define ITOT_LOOP(i)  for ((i) = 0; (i) < NX1_TOT; (i)++)
#define JTOT_LOOP(j)  for ((j) = 0; (j) < NX2_TOT; (j)++)
#define KTOT_LOOP(k)  for ((k) = 0; (k) < NX3_TOT; (k)++)

#define DOM_LOOP(k,j,i) KDOM_LOOP(k) JDOM_LOOP(j) IDOM_LOOP(i)

#define TOT_LOOP(k,j,i) KTOT_LOOP(k) JTOT_LOOP(j) ITOT_LOOP(i)

#define X1_BEG_LOOP(k,j,i) KTOT_LOOP(k) JTOT_LOOP(j) IBEG_LOOP(i)
#define X2_BEG_LOOP(k,j,i) KTOT_LOOP(k) JBEG_LOOP(j) ITOT_LOOP(i)
#define X3_BEG_LOOP(k,j,i) KBEG_LOOP(k) JTOT_LOOP(j) ITOT_LOOP(i)

#define X1_END_LOOP(k,j,i) KTOT_LOOP(k) JTOT_LOOP(j) IEND_LOOP(i)
#define X2_END_LOOP(k,j,i) KTOT_LOOP(k) JEND_LOOP(j) ITOT_LOOP(i)
#define X3_END_LOOP(k,j,i) KEND_LOOP(k) JTOT_LOOP(j) ITOT_LOOP(i)

#define TRANSVERSE_LOOP(indx, in, i,j,k) \
 if (g_dir == IDIR) {i = &in; j = &indx.t1; k = &indx.t2;} \
 if (g_dir == JDIR) {j = &in; i = &indx.t1; k = &indx.t2;} \
 if (g_dir == KDIR) {k = &in; i = &indx.t1; j = &indx.t2;} \
 g_i = i; g_j = j; g_k = k;\
 for (indx.t2 = indx.t2_beg; indx.t2 <= indx.t2_end; indx.t2++) \
 for (indx.t1 = indx.t1_beg; indx.t1 <= indx.t1_end; indx.t1++)

/* ********************************************************************* */
#define BOX_LOOP(B,k,j,i) \
 for ((B)->dk = ((k=(B)->kb) <= (B)->ke ? 1:-1); k != (B)->ke+(B)->dk; k += (B)->dk)\
 for ((B)->dj = ((j=(B)->jb) <= (B)->je ? 1:-1); j != (B)->je+(B)->dj; j += (B)->dj)\
 for ((B)->di = ((i=(B)->ib) <= (B)->ie ? 1:-1); i != (B)->ie+(B)->di; i += (B)->di)


/*
#define DIR_LOOP(grid,in,k,j,i) \
    if (g_dir == IDIR) {in = &i; j = *g_j; k = *g_k;} \
    if (g_dir == JDIR) {in = &j; i = *g_i; k = *g_k;} \
    if (g_dir == KDIR) {in = &k; i = *g_i; j = *g_j;} \
    for (*in = grid[g_dir].lbeg - 2; *in <= grid[g_dir].lend + 2; (*in)++)
*/

/*  not used for now... will serve in EMF_BOUNDARY...
#define ISTAG_LOOP(i)  for ((i) = IBEG - 1; (i) <= IEND; (i)++)
#define JSTAG_LOOP(j)  for ((j) = JBEG - 1; (j) <= JEND; (j)++)
#define KSTAG_LOOP(k)  for ((k) = KBEG - 1; (k) <= KEND; (k)++)
*/


/* ################################################################# 

                   VARIABLE LOOP MACROS        

   ################################################################# */

 /* ------------------------------------------------------------
     The LRVAR_LOOP and LRFLX_LOOP will speed up computation 
     of the left/right states at a given interface by excluding 
     the normal component of magnetic field, which is usually  
     assigned at the end of the reconstruction step.
    ------------------------------------------------------------ */
/*      
#ifdef STAGGERED_MHD
 #define LR_VAR_LOOP(i) for (i = NVAR; i--; i -= (i == BXn+1)) 
 #define LR_FLX_LOOP(i) for (i = NFLX; i--; i -= (i == BXn+1)) 
#else
 #define LR_VAR_LOOP(i) for (i = NVAR; i--; ) 
 #define LR_FLX_LOOP(i) for (i = NFLX; i--; ) 
#endif
*/
/* ################################################################# 

                  TIME INTEGRATOR DEFINITIONS

   ################################################################# */

#if (TIME_STEPPING == HANCOCK) || (TIME_STEPPING == CHARACTERISTIC_TRACING)
 #define SINGLE_STEP   1
 #if DIMENSIONAL_SPLITTING == NO
  #define CTU      1    /* -- Corner Transport Upwind method of Colella -- */
 #endif
#endif

/* ------------------------------------------------------------
    the GET_MAX_DT switch determines how the time step is
    computed. For pure advection, setting it to YES will force
    PLUTO to compute dt in the old way, i.e., by taking the
    maximum.
    When set to NO, the time step is computed only during the
    predictor step and, for UNSPLIT RK schemes, it will be
    calculated as the average over dimensions resulting in 
    slightly larger time increments.
   ------------------------------------------------------------ */

#if ((TIME_STEPPING == RK2) || (TIME_STEPPING == RK3) || \
     (TIME_STEPPING == RK_MIDPOINT)) && DIMENSIONAL_SPLITTING == NO

 #define GET_MAX_DT    NO
#else
 #define GET_MAX_DT    YES
#endif

/* ############################################################### 
  
          Check whether a separate source step has to be done

   ############################################################### */

#if (COOLING != NO) 
 #define INCLUDE_SPLIT_SOURCE   YES
#else
 #define INCLUDE_SPLIT_SOURCE   NO
#endif


/* #################################################################

    set to NO all undefined identifier that may appear in the code.
    (compile with -Wundef).
    Diffusion operators: PARABOLIC_FLUX is the bitwise OR
    combination of all operators, each being either one of 
    NO, EXPLICIT (1st bit), STS (2nd bit). 
    It can take the following values

      00   --> no diffusion operator is being used
      01   --> there's at least one explicit diffusion operator and
               no sts.
               
      10   --> there's at least one sts diffusion operator and
               no explicit one.
      11   --> mixed: there is at least one explicit and sts operator
     
   ################################################################# */

#ifndef RESISTIVE_MHD 
 #define RESISTIVE_MHD   NO
#endif

#ifndef THERMAL_CONDUCTION
 #define THERMAL_CONDUCTION NO
#endif

#ifndef VISCOSITY
 #define VISCOSITY NO
#endif

#ifndef INCLUDE_PARTICLES
 #define INCLUDE_PARTICLES NO
#endif

#ifndef UPDATE_VECTOR_POTENTIAL
 #define UPDATE_VECTOR_POTENTIAL  NO
#endif

#ifndef BACKGROUND_FIELD
 #define BACKGROUND_FIELD NO
#endif

#ifndef USE_FOUR_VELOCITY
 #define USE_FOUR_VELOCITY NO
#endif

#ifndef CHAR_LIMITING
 #define CHAR_LIMITING  NO
#endif

#ifndef ARTIFICIAL_VISCOSITY 
 #define ARTIFICIAL_VISCOSITY  NO
#endif

#ifndef ROTATING_FRAME
 #define ROTATING_FRAME NO
#endif

#ifndef INTERNAL_BOUNDARY
 #define INTERNAL_BOUNDARY NO
#endif

#define PARABOLIC_FLUX (RESISTIVE_MHD|THERMAL_CONDUCTION|VISCOSITY)

/* ################################################################# 

                       USEFUL CONSTANTS        

   ################################################################# */

#define ACCURACY    1.e-14          /* ZERO ACCURACY */

#define CONST_PI      3.14159265358979   
#define CONST_amu     1.66053886e-24     
#define CONST_mp      1.67262171e-24     
#define CONST_mn      1.67492728e-24     
#define CONST_me      9.1093826e-28      
#define CONST_mH      1.6733e-24         
#define CONST_kB      1.3806505e-16      
#define CONST_sigma   5.67051e-5         
#define CONST_sigmaT  6.6524e-25         
#define CONST_NA      6.0221367e23       
#define CONST_c       2.99792458e10      
#define CONST_Msun    2.e33              
#define CONST_Rsun    6.96e10            
#define CONST_Mearth  5.9736e27          
#define CONST_Rearth  6.378136e8         
#define CONST_G       6.6726e-8          
#define CONST_h       6.62606876e-27     
#define CONST_pc      3.0856775807e18    
#define CONST_ly      0.9461e18          
#define CONST_au      1.49597892e13      
#define CONST_eV      1.602176463158e-12 
/* ############################################################### 
 
                     Define Structures here     

   ############################################################### */

typedef struct CMD_LINE {
  int restart;       
  int h5restart;       
  int nrestart;      
  int makegrid;      
  int write;         
  int maxsteps;
  int jet;  /* -- follow jet evolution in a given direction -- */
  int parallel_dim[3];
  int nproc[3];  /* -- user supplied number of processors -- */
  int show_dec; /* -- show domain decomposition ? -- */
  int xres; /* -- change the resolution via command line -- */
  int fill; /* useless, it makes the struct a power of 2 */ 
} Cmd_Line;
   
/* ********************************************************************* */
typedef struct DATA{
  double ****Uc;
  double ****Vc;  
  double ****Vs;  
  double ****Vuser; 
  double ***Ax1; 
  double ***Ax2; 
  double ***Ax3; 
  unsigned char ***flag; 
  int fill; /* useless, it makes the struct a power of 2 */ 
} Data;
   
/* ********************************************************************* */
typedef struct STATE_1D{
  double **v;    
  double **vL;   
  double **vR;   
  double **vm;   
  double **vp;   
  double **uL;   
  double **uR;   
  double **um;   
  double **up;   
  double **flux;    
  double **par_flx; 
  double **pnt_flx;
  double **dff_flx;

  double ***Lp, ***Rp; 
  double **lambda;     
  double *lmax;   
  double *a2;     
  double *h;      
  double **vt;
  double **src;     
  double **par_src; 
  double **vh;      
  double **rhs;     
  double *press;    
  double *bn;       
  double *SL;       
  double *SR;       
  double q1;            
  unsigned char uc1;    
  unsigned char *flag;
} State_1D;

/* ********************************************************************* */
typedef struct GRID{
  double xi, xf;        
  double *x, *x_glob;   
  double *xr, *xr_glob; 
  double *xl, *xl_glob; 
  double *dx, *dx_glob; 
  double *xgc;          
  double *dV;           
  double *A;            
  double *r_1;          
  double *ct;           
  double *inv_dx;       
  double *inv_dxi;      
  double dl_min;      
  double *dfg;  
  double *df2g; 
  double *dfR;  
  double *dfL;  
  int np_tot_glob; 
  int np_int_glob; 
  int np_tot;      
  int np_int;      
  int nghost;      
  int lbound;      
  int rbound;      
  int gbeg;        
  int gend;        
  int beg;         
  int end;         
  int lbeg;        
  int lend;        
  int uniform;     /* = 1 when the grid is cartesian AND uniform everywhere  */
  int nproc;       
  int rank_coord;  
  int level;       
  char fill[84];   /* useless, just to make the structure size a power of 2 */
} Grid;

/* ********************************************************************* */
typedef struct TIME_STEP{
  double *cmax;     
  double inv_dta;   
  double inv_dtp;   
  double dt_cool;   
  double cfl;       
  double cfl_par;   
  double rmax_par;  
  int    Nsts;      
  int    Nrkc;      
  char  fill[24];   /* useless, just to make the structure size a power of 2 */
} Time_Step;


/* ********************************************************************* */
typedef struct OUTPUT{
  int    type;       
  int    nvar;       
  int    user_outs;  
  int    nfile;      
  int    dn;         
  int    *stag_var;  
  int    *dump_var;  
  char   mode[32];   
  char   **var_name; 
  char   ext[8];     
  double dt;         
  double dclock;     
  double ***V[64];   
  char   fill[168];  
} Output;

typedef struct INPUT{
  int    npoint[3];           
  int    lft_bound_side[3];       /* left  boundary type */
  int    rgt_bound_side[3];      /* right boundary type */
  int    grid_is_uniform[3];    /* = 1 when grid is uniform, 0 otherwise */
  int    npatch[5];               /* number of grid patches  */
  int    patch_npoint[5][16];     /* number of points per patch */
  int    patch_type[5][16];             
  int    log_freq;              /* log frequency */
  int    user_var;              /* number of additional user-variables being
                                   held in memory and written to disk */
  int    anl_dn;                /*  number of step increment for ANALYSIS */
  char   solv_type[64];
  char   user_var_name[128][128];
  Output output[MAX_OUTPUT_TYPES];  
  double patch_left_node[5][16];  /*  self-expl. */
  double  xbeg[3];
  double  xend[3];
  double  cfl;                    /* hyperbolic cfl number */
  double  cfl_max_var;
  double  cfl_par;           /* (STS) parabolic  cfl number */
  double  rmax_par;          /* (STS) max ratio between current time
                                step and parabolic time step */
  double  tstop;
  double  first_dt;
  double  anl_dt;          /* time step increment for ANALYSIS */
  double  aux[32];         /* we keep aux inside this structure, 
                              since in parallel execution it has
                              to be comunicated to all processors  */
} Input;

typedef struct RUNTIME{
  int nstep;
  int nfile[MAX_OUTPUT_TYPES];
  double t;
  double dt;
} Runtime;

typedef struct RGB{
  unsigned char r, g, b;
} RGB;

typedef struct IMAGE{
  int    nrow, ncol;    /* -- image rows and columns -- */
  int    slice_plane;   /* -- one of X12_PLANE, X13_PLANE, X23_PLANE -- */
  int    logscale;      /* -- YES/NO for log scale -- */
  char   *colormap;     /* -- colormap name -- */
  char   basename[32];  /* -- image base name (no extensions) -- */
  unsigned char r[256], g[256], b[256]; /* -- colortable saved here -- */
  RGB    **rgb;         /* -- rgb array containing image values -- */
  double max;           /* -- max image value -- */
  double min;           /* -- min image value -- */
  double slice_coord;   /* -- slice coord orthogonal to slice_plane -- */
} Image;

typedef struct FLOAT_VECT{
  float v1, v2, v3;
} Float_Vect;

typedef struct INDEX{
  int ntot, beg, end;
  int t1, t1_beg, t1_end;
  int t2, t2_beg, t2_end;
  char fill[28]; /* useless, just to make the structure size a power of 2 */
} Index;

/* ********************************************************************* */
typedef struct RBOX{
  int ib; 
  int ie; 
  int jb; 
  int je; 
  int kb; 
  int ke; 
  int di; 
  int dj; 
  int dk; 
  int vpos; 
} RBox;

/* ####################################################### 

            Recurrent function types 
 
   ####################################################### */

typedef void Riemann_Solver (const State_1D *, int, int, double *, Grid *);
typedef void Limiter        (double *, double *, double *, int, int, Grid *);
typedef double Reconstruct  (double *, double, int);
typedef double ****Data_Arr;

/* -- when using Finite Difference Schemes, the "Riemann Solver" 
      function computes the fluxes with high order interpolants. -- */

/* ----------------------------------------------------------
             Include module header files
   ---------------------------------------------------------- */

#if   PHYSICS == HD

 #include "HD/mod_defs.h"

#elif PHYSICS == RHD

 #include "RHD/mod_defs.h"

#elif PHYSICS == MHD

 #include "MHD/mod_defs.h"
 #ifdef SHEARINGBOX
  #include "MHD/ShearingBox/shearingbox.h"
 #endif

#elif PHYSICS == RMHD

 #include "RMHD/mod_defs.h"

#endif

#if COOLING == MINEq
 #include "Cooling/MINEq/cooling.h"
#elif COOLING == SNEq 
 #include "Cooling/SNEq/cooling.h"
#elif COOLING == POWER_LAW
 #include "Cooling/Power_Law/cooling.h"
#elif COOLING == TABULATED
 #include "Cooling/Tab/cooling.h"
#elif COOLING == H2_COOL
 #include "Cooling/H2_COOL/cooling.h"
#endif

#if THERMAL_CONDUCTION != NO
 #include "Thermal_Conduction/tc.h"
#endif

#if VISCOSITY != NO
 #include "Viscosity/viscosity.h"
#endif

#if INCLUDE_PARTICLES == YES
 #include "Particles/particles.h"
#endif

#ifdef FARGO
 #include "Fargo/fargo.h"
#endif

/* *****************************************************
   *****************************************************

                  GLOBAL VARIABLES

   *****************************************************
   ***************************************************** */

 extern int SZ;
 extern int SZ_stagx;
 extern int SZ_stagy;
 extern int SZ_stagz;
 extern int SZ_char;
 extern int SZ_float;
 extern int SZ_Float_Vect;
 extern int SZ_rgb;
 extern int SZ_short;
 extern int prank;

extern long int IBEG, IEND, JBEG, JEND, KBEG, KEND;
extern long int NX1, NX2, NX3;
extern long int NX1_TOT, NX2_TOT, NX3_TOT;

extern long int NMAX_POINT;

extern int VXn, VXt, VXb;
extern int MXn, MXt, MXb;
extern int BXn, BXt, BXb;

extern int *g_i, *g_j, *g_k;

extern int g_dir;
extern int g_maxRiemannIter;
extern long int g_usedMem;
extern long int g_stepNumber;
extern int      g_intStage;

extern double g_unitDensity, g_unitLength, g_unitVelocity;

extern double g_maxCoolingRate, g_minCoolingTemp;

extern double g_smallDensity, g_smallPressure;

extern double g_time, g_dt;
extern double g_maxMach;
#if ROTATING_FRAME
 extern double g_OmegaZ;
#endif

extern double g_domBeg[3], g_domEnd[3];

extern double g_inputParam[32];
#if EOS != ISOTHERMAL
 extern double g_gamma;
#elif EOS == ISOTHERMAL
 extern double g_isoSoundSpeed;
#endif

#ifdef CH_SPACEDIM
 extern double ch_max, ch_max_loc, coeff_dl_min;
 extern int use_glm;
 extern int LEVEL;
#endif

extern FILE *pluto_log_file;

/* ----------------------------------------------------------
    define NX_MAX, NY_MAX, NZ_MAX as the maximum 
    size on a patch.
    For the static grid of PLUTO, they simply coincide
    with NX1_TOT, NX2_TOT, NX3_TOT. Otherwise
    the grid is taken to be the maximum, according
    to maxGridSize. 
    This avoids allocating and freeing memory all the times.
   ---------------------------------------------------------- */

#ifdef CH_SPACEDIM
 #define NX_MAX NMAX_POINT
 #define NY_MAX NMAX_POINT
 #define NZ_MAX NMAX_POINT
#else
 #define NX_MAX NX1_TOT
 #define NY_MAX NX2_TOT
 #define NZ_MAX NX3_TOT
#endif 

/* ------------------------------------------------------------------
    Define total number of variables to be integrated;
    this includes:

     NFLX     = Number of equations defining the system of 
                conservation laws. For example, for the 
                HD module, it consists of density, momentum and energy.
                It is defined in the physics module header file mod_defs.h.

     NTRACER  = Number of user-defined tracers; defined in the problem
                directory header file definitions.h

     NIONS = Number of chemical species; it is defined in 
                the cooling modules cooling.h, if present.

  nv = 0...NFLX - 1; NFLX...NFLX+NIONS-1; TRC...TRC + NTRACER-1; ENTR; 
       <----------> <------------------> <------------------>  
          NFLX              NIONS               NTRACER              

                    <---------------------------------------------->
                                    NSCL

   ------------------------------------------------------------------ */

#ifndef ENTROPY_SWITCH
 #define ENTROPY_SWITCH  NO 
#endif 

#ifndef NIONS
 #define NIONS 0
#endif

#define TRC   (NFLX + NIONS)
#define TR   TRC
#define NSCL (NTRACER + NIONS + ENTROPY_SWITCH)

#if ENTROPY_SWITCH == YES
 #define ENTR  (TRC + NTRACER)
#else
 #if EOS != ISOTHERMAL && EOS != BAROTROPIC
  #define ENTR (ENG)
 #endif
#endif

#define NVAR (NFLX + NSCL)

#define KELVIN (g_unitVelocity*g_unitVelocity*CONST_amu/CONST_kB) 

#include "prototypes.h"
#endif /* PLUTO_H */