Home

User's Guide

1. CHAPTER 2 PROBLEM HEADER FILE DEFINITIONS H 28 Figure 2 3 Second order TVD lim iter functions used by PLUTO as functions of the left to right slope ra tio z AV 1 AV 1 Larger values of lim x indicate larger com pressive behavior In this sense the minmod limiter MM x and the monotonized central limiter MC x are the least and most compressive respectively where MINMOD_LIM is the most diffusive and MC_L IM is the least diffusive limiter The TVD dia eram for the various limiter functions is shown in Fig 2 3 All of the above limiters employ a 3 point stencil e CT_EMF_AVERAGE string Control how the electromotive force EMF is integrated from the face center to the edges This is discussed in more detailed in 46 2 3 3 e CT_EN_CORRECTION YES NO This option is available only in the MHD and RMHD modules The default is NO implying that energy is not corrected after the conservative update However for low beta plasma one may find useful to switch this option to YES as described in BS99 e ASSIGN VECTOR POTENTIAL YES NO When set to YES magnetic field components are initialized from the vector potential In the con strained transport algorithm CT 6 2 3 3 this guarantees that the magnetic field has zero diver gence When set to NO assignment proceeds in the usual way see 46 2 2 for more details e UPDATE_VECTOR_POTENTIAL YES NO Enable this option if you wish to evolve the vector
2. gs xl rs gs x2 rs gsxx3 rs elif GEOMETRY CYLINDRICAL g IDIR g JDIR g KDIR gsxx1 rs gSs Xx2 rs 0 0 elif GEOMETRY SPHERICAL g IDIR g IDIR g KDIR endif gs 0 0 0 0 It is also possible to prescribe the body force in terms of a vector and a potential by setting in your definitions h BODY_FORCE to VECTOR POTENTIAL Beware that non intertial effects due to a rotating frame of reference such as Coriolis and centrifugal forces should not be specified here since they are automatically handled by PLUTO by enabling the ROTATING_FRAME flag in the HD and MHD module see A word of caution about using reflective equatorial symmetric or similar boundary conditions strictly speaking gravity should be defined consistently with the antisymmmetric behavior of the veloc ity component normal to the given boundary plane More precisely the normal component of g should CHAPTER 5 INITIAL AND BOUNDARY CONDITIONS INIT C 55 be antisymmetric while the potential should be an even function about the boundary plane Consider for instance a reflective or equatorial symmetric conditions at the lower and upper boundaries 2 and ze in the z direction Then one should have O z y 2 2 y 2 2 O z y Z z P 2 y 2 z gz Y Zp 2z g2 y 2 l ABU If gravity does not satisfy this property then it must be imposed manually As an example you can look at the Test_
3. When set to MULTID a multi dimensional strategy is used by which upon shock detection i reconstruction in every direction reverts to the MINMOD limiter and ii fluxes are computed using the HLL Riemann solver The flagging strategy is set in States flag shoock c The MULTID shock flattening has proven to be an effective adaptation strategy that can noticeably increase the code robustness It is highly suggested for complex flow structures involving strong shocks CHAR_LIMITING YES NO Set it to YES to perform reconstruction on characteristic variables rather than primitive It is avail able for the HD RHD and MHD modules Although somewhat more expensive reconstruction in characteristic ariables is known to produce better quality solutions by suppressing unwanted nu merical oscillation in proximity of strong discontinuities and leading to a better entropy enforce ment LIMITER string Sets the limiter s to be applied when RECONSTRUCTION is set to LINEAR Here string can be one of FLAT_LIM set slope to zero 1 order reconstruction MINMOD_LIM use the minmod limiter most diffusive VANALBADA_LIM use the van Albada limiter function OSPRE_LIM use the OSPRE limiter VANLEER_LIM use the harmonic mean limiter of van Leer UMIST_LIM use the umist limiter MC_LIM use the monotonized central difference limiter least diffusive DEFAULT keep the default setting defined in Src States set_limiter c
4. BOX_LOOP box k J i bxsout d gt Vs BX1s k 2 JBEG j 111115 d gt Vs BX1s k i a Bxscut 1 0 Profile rpli BX endif Here STAGGERED_MHD is defined only in the MHD constrained transport and the boundary conditions are assigned on by bp only i e the orthogonal component 6 2 4 Background Field Splitting In situations where an intrinsic background magnetic field is present e g planetary magnetosphere stellar dipole fields it may be convenient to write the total magnetic field as B x t Bo x B x t where Bo is a background curl free magnetic field and B x t is a deviation The background field must satisfy the following conditions o Bo Ot In this case one can show Pow94 that the MHD equations reduce to O V Bo VxBo 0 Op ay V ou 0 O V mo BiB EE BB T Vo p V 9 O E MERO Ly Es pe pb w B1 B mg 2 TT 0 Ot CHAPTER 6 BASIC PHYSICS MODULES 64 where r p p 2i 4 Bi Bo E P pu B3 Thus the energy depends only on B a feature that turns out to be useful when dealing with low beta plasma The sets of conservative and primitive variables are the same as the original ones with B gt By E gt Es In order to enable this feature the macro BACKGROUND_FIELD must be turned to YES in your definitions h The initial and boundary conditions must be imposed on B alone while the function BackgroundField c
5. T 100 No Reaction Rate Coefficient cm s Reference d L Ht e gt Ht4 2e k 5 84x 10 T exp 157 809 0 T 1 2 Ht e 3 H hv ky 2 6x107 UT 9 1 3 H2 e gt 2H e kg 4 4 x 10710703 exp 102 000 0 T 2 4 H H gt 3H y 1 067 a exp 4 463 Toy 1 0 2472T 1 3 512 3 5 H gt 2 H gt H gt 2H ks 1 0 x 10 8exp 84 100 T 4 6 H H hH ke 3 0 x 1071 T2 1 0 0 4yTz 0 2T gt 0 08T2 5 REFERENCES 1 Eq 1e 2 Eq H17 3 AAZN97 Tab 3 Eq 13 4 WAMMO07 UMIST Database 5 Eq 3 8 9 4 Molecular Hydrogen Non Equilibrium Cooling H2_COOL This module is implemented in the Src Cooling H2_COOL directory and introduces three new variables with index X_HI X_H2 and X_HII used to label the fraction of atomic hydrogen molecular hydrogen and ionized hydrogen respectively as follows NH na NA rr LH2 LHI THIIS nH nH nH where the total hydrogen number density ny ny Np 2Np You can assign these hydrogen fractions in a similar manner like the SNEq module v X_HI v X_H2 0 4 v X_HII 1 0 v X_HI 2 0x xv X_H2 0 2 x for example in your Init function Note the value of v X_H2 should be between 0 0 and 0 5 while the remain ing two hydrogen fractions can have values ranging from 0 0 to 1 0 such that their sum is conserved The chemical evolution of molecular atomic and ionized hydrogen is governed by equations listed in
6. When set to YES initial conditions are assigned by sub sampling and averaging different values inside each cell It is useful to create smooth profiles of sharp boundaries not aligned with the grid e g a circle in Cartesian coordinates WARNING_MESSAGES YES NO Issue a warning message every time a numerical problem or inconsistency is encountered setting WARNING_MESSAGES to YES will tell PLUTO to print what when and where a numerical problem occurred PRINT TOFILE YES NO When set to YES it tells PLUTO to re direct the output log to the file pluto log If this file does not exist it will be created if the file exists but integration starts from initial conditions it will be over written Finally if you restart from a previously saved file the output will be appended INTERNAL BOUNDARY YES NO When turned to YES it allows to overwrite or change the solution array anywhere inside the computational domain This is done inside the UserDefBoundary function when side 0 see This option is particularly useful when flow variables must be kept constant in time or to assign lower upper threshold values to any physical quantity e g density or pressure SHOCK_FLATTENING NO ONED MULTID Provides additional dissipation in proximity of strong shocks When set to ONED spatial slopes are progressively reduced following following a one dimensional shock recognition pattern as in ICW84 This is done separately dimension by dimension
7. The TAUB Equation ol statel espresso oe Od Rw EES oe a a HS 71 72 OL VISCOSIl 26 25 44 942558555645 R28 n gt oe POSSESS EOE EG ea ee Pee 72 02 BoD eea a eee a A A 74 EEEE a rd EENT 75 SL Dimension rar AAA AAA AA 75 8 4 Numerical Integration of Diffusion Terms ooa 77 S41 Explicit Dime Stepping s s eea a a sk e er 77 8 4 2 Super Time Stepping STS o ee 77 9 Optically Thin Cooling 79 91 Power Law A AA IEA 80 9 2 Tabulated COLMO eros orar AAA ASA 81 9 3 Simplified Non Equilibrium Cooling SNEq 82 9 4 Molecular Hydrogen Non Equilibrium Cooling H2_COOL o o ooo o o 83 9 5 Multi lon Non Equilibrium Cooling MINEG 10 Additional Modules 86 TRE ee E 86 A 87 aora aos nerds eee eae eee 88 acosa aaa AA es 88 aa oras aorta seas 88 arar aaa AA 90 10 3 1 WENO schemest 0 00 ee ee eee 90 10 3 2 LimO3 amp MP5 002000 020 000 0000 2 ee ee ee 91 92 UA Ouiput Data Formats 425 tbs eyo bbe be Roo eae ee ee ea ob oe oe 92 11 1 1 Binary Output dbl or flt data formats o 93 AI 93 11 1 3 VTK Output vtk data format o e 94 11 1 4 ASCII Output tab Data format 4 ck neha KAR AAA 94 E a e aga 94 T TETS TE TETTERE ETETETT 95 11 2 Customizing your output s were ae eb ebb AA 96 ANITA 96 SA
8. optional flt h5 LOL i optional tab e 1 optional ppm Sal optional png nO optional log 1 analysis 1 0 1 optional Chombo HDF5 output Output_dir wif optional Checkpoint_interval 1 0 0 Plot_interval 1 0 0 Parameters SCRH 0 32 CHAPTER 4 RUNTIME INTITALIZATION FILE PLUTO INI 33 Tag labels on the left side identify appropriate field s following on the same line and must not be changed Block ordering is irrelevant Runtime parameters can be accessed anywhere in the code through the members of the Runtime structure see Doc Doxygen html structs_8h html using the func tion RuntimeGet e g cfl RuntimeGet gt cfl x Obtain the cfi number x char fname 64 sprintf fname s mydata dat RuntimeGet gt output_dir Prepend output dir to file name The quantities and related data types read from the file are now described 4 1 The Grid block In the static version of PLUTO it defines the physical domain and generates the whole computation mesh while in the AMR version is used to specify the base grid corresponding to level 0 e Xl grid integer double integer char double e X2 grid integer double integer char double e X3 grid integer double integer char double For each dimension the first integer defines the number of non overlapping adjacent one dimensional grid patches making up the computational domain Note
9. 3 4 5 YES NO int int YES NO YES NO YES NO real real YES NO YES NO YES NO YES NO real Description Sets the order of reconstruction when using PARABOLIC recon struction Allowed values are 3 uses three point stencil third order accurate 4 fourth order and 5 fifth order For more information see Mig 14 The default value is 4 as in the origi nal PPM method Used for the PVTE_LAW EOS in ionization equilibrium When set to YES replaces function evaluations of the thermal EOS p nkgT and its inverse with lookup table and bilinear interpolation This results in a considerably faster execution Default is YES Sets the number of x points used to construct the temperature table for the PVTE_LAW EOS The default value is set in Sr c EOS PVTE thermal_eos c Sets the number of y points used to construct the temperature table for the PVTE_LAW EOS The default value is set in Sr c EOS PVTE thermal_eos c Use the four velocity u yv instead of the three velocity when reconstructing left and right states in the RHD and RMHD modules Not compatible with conservative form of MUSCL Hancock scheme Default is NO If set to YES use approximate and faster expressions when computing the fast magnetosonic speed of the RMHD equa tions see Sect 3 3 of DZBLO7 Solutions of quartic equa tion is avoided and replaced with upper bounds provided by quadratic equation
10. PLOAD tries to read binary data in double precision if dbl out is present To select a different format a corresponding keyword must be supplied e g FLOAT H5 or HDF5 or a combination of them see Doc idl_tools html When PLOAD is called for the first time it initializes the following four common blocks e PLUTO_GRID contains grid information such as the number of points nx1 nx2 nx3 coordi nates x1 x2 x3 and mesh spacing dx1 dx2 dx3 e PLUTO_VAR the number NVAR and the names of variables being written for the chosen format Variable names follow the same convention adopted in PLUTO e g rho vxl vx2 bx1 bx2 prs and so on e PLUTO_RUN time stepping information such as output time t time step dt and total number of files nlast PLOAD can be used inside a normal IDL script after it has been invoked at least once or compiled with Y pload The DISPLAY Procedure DISPLAY is a general purpose visualization routine and the source code can be found in Tools IDL display pro The DISPLAY procedure can be used to display the intensity map of a 2D variable to a graphic window e g IDL gt PLOAD 3 load the third data set in double precision IDL gt DISPLAY x1 x1 x2 x2 alogl0 rho title Density vbar IDL gt DISPLAY xl x1 x2 x2 vx1 title X Velocity nwin 1 The second line displays the density logarithm and the third line displays the x component of velocity in a new window Another
11. RBox xbox int side Grid xgrid int i j k if side 0 TOT LOOP k j i k i if d gt Vc RHO d gt Vc RHO k ILLL 31 3 A more complex example consists of a time independent region of space where variables are fixed in time and should not be evolved by the algorithm If this is the case you may additionally tell PLUTO not to update the solution in the specified computational zones during the current time step by enabling the FLAG_INTERNAL_BOUNDARY flag CHAPTER 5 INITIAL AND BOUNDARY CONDITIONS INIT C 93 Example The following example taken from Test_Problems HD Stellar_Wind shows how to set up a radially sym metric spherical wind in cylindrical coordinates inside a small spherical region of radius 1 centered around the origin This is achieved by prescribing fixed inflow values for density pressure and velocity void UserDefBoundary const Data d RBox box int side Grid grid Int ip Jy Ky nv double xl Xx2 Xx3 double r r0 cs double Vwind 1 0 rho Vr grid IDIR xgc grid JDIR xgc grid KDIR xgc if side 0 EO ds cs g_inputParam CS_WIND TOT_LOOP k j 1 r sqrt xl a xt 1 x2 Jj x2 J if r lt 40 4 vr tanh r r0 0 1 Vwind rho Vwind r0x r0 vrxrx xr HO k j i rho j Vwindx x1 1 r Vwindxx2 3 r cs cs g_gammaxpow rho g_gamma FLAG INTERNAL BOUNDARY The symbol a combination of the bitwise OR ope
12. amp scrh 1 MPI_DOUBLE MP1_MAX MPI_COMM WORLD Eth m x scrh MPT Barrier MPI_COMM_WORLD endif Write ascii file averages dat to disk x if prank 0 char fname 512 static double tpos 1 0 FILE xfp sprintf fname s averages dat RuntimeGet gt output_dir if g_stepNumber 0 Open for writing only when we re starting fp fopen fname w x from beginning fprintf fp _ 37s _ 12s 12s 12s n t dt lt Ekin gt Max Eth else x Append if this is not step 0 x if tpos lt 0 0 Obtain time coordinate of to last written row x char sline 512 fp fopen fname r while fgets sline 512 fp sscanf sline slf n amp tpos tpos time of the last written row x fclose fp fp fopen fname a if g_time gt tpos x Write if current time if gt tpos x fprintf fp Sl2 6e 3l2 6e 212 66 12 66 1n q_tim e g_dt Ekin Eth max fclose fp 6 Basic Physics Modules In this chapter we describe the basic equation modules available in the PLUTO code for the solution of the fluid equations under different regimes HydroDynamics HD MagnetoHydroDynamics MHD and their relativistic extensions RHD and RMHD We remind that only first order spatial derivatives accounting for the hyperbolic part of the equations are described in this chapter whereas the reader is referred to Chap S for a comprehensive descrip
13. instead of pluto log where n 0 Np 1 and Np is the total number of processors However pout 0 also contains additional information regarding the chosen configuration Pre configured examples can be found in the Test_Problems folder e g e Configuration 04 of Test_Problems HD Mach_Reflection e Configuration 04 of Test Problems HD Stellar_Wind e Configuration 03 of Test_Problems RHD Blast e Configuration 08 of Test_Problems MHD Rayleigh_Taylor e Configuration 07 08 11 of Test_Problems MHD Rotor e Configuration 08 of Test_Problems MHD Shock_Cloud e Configuration 03 of Test_Problems RMHD KH e Configuration 01 02 of Test_Problems RMHD Shock_Cloud CHAPTER 12 ADAPTIVE MESH REFINEMENT AMR 110 12 3 Controlling Refinement Refinement is controlled by a series of runtime parameters defined in the Chombo Refinement block 34 2 of your pluto ini Zones are tagged for refinement whenever a prescribed function x U of the conserved variables and of its derivatives exceeds the threshold value x assigned to Refine_thresh in your pluto ini Generally speaking the refinement criterion may be problem dependent thus requiring the user to provide an appropriate definition of x U The default criterion is based on the second derivative error norm Loh87 and it is implemented in the function computeRefGradient in the source file Src Chombo TagCells cpp The test function adopted for this purpose is Xa Ag 1 U gt
14. lt KEND but it is considered part of the solution Similar considerations hold for by and bz components at the x2 and x3 boundaries respectively The electromotive force EMF is computed at zone edges see Fig 6 1 by a proper averaging re construction scheme set by CT_EMF_AVERAGE inside your definitions h Options are e CT_EMF_AVERAGE ARITHMETIC yields a simple arithmetic averaging BS99 of the fluxes com puted during the upwind steps In this case one has available 0 Cus Cae gt 0 Cii Es Distr ATEN i ij ki 2 2 2 during the z1 2 and x3 sweeps respectively The arithmetic average follows 1 Cajal pet 4 are EE R T dia CHAPTER 6 BASIC PHYSICS MODULES 62 Figure 6 1 Collocation points in 2 x D left and in 3D right Cell centered quantities are given as green circles face centered as red squares and edge centered as blue diamonds 1 Cee itt gel A CA a Cen ith jk T go i j kth al 1 xs ith jth h 4 CA a Coar el h a Cig ey T Eia Although being the simplest one this average procedure may suffer from insufficient dissipation in some circumstances GS05 Ld04 and does not reduce to its one dimensional equivalent algo rithm for plane parallel grid aligned flows e CT EMF AVERAGE UCT_HLL uses a two dimensional Riemann solver based on a four state HLL flux function see DBLO03 Ld04 If the fully unsplit HANCOCK or CHARACTRISTIC_TRACING scheme is used the Courant number
15. right panel collocation of physical variables on a 2D grid X and Y face centered staggered quantities are shown by squares and triangles respectively Filled symbols circles boxes and triangles are considered interior values part of the solution whereas boundary values are identified by empty symbols and must be prescribed by the user if the boundary is userdef X2_BEG X2_END boundary conditions can be assigned in the ghost zones at the beginning and end of the physical domain in the x2 direction X3_BEG X3_END boundary conditions can be assigned in the ghost zones at the beginning and end of the physical domain in the 3 direction 0 zero change control the solution inside the computational domain This feature can be used only if the macro INTERNAL_BOUNDARY has been enabled in your definitions h see 5 2 1 If say X1 beg has been tagged userdef inside your pluto ini the user has to specify the boundary values at the beginning of the x direction when side X1_BEG e xgrid a pointer to an array of Grid structures containing all the relevant grid information In this case grid IDIR is the structure relevant to the x direction ygrid JDIR to the z direction and grid KDIR pertains to the x3 direction See the code documentation for more details on the members of the Grid structure CHAPTER 5 INITIAL AND BOUNDARY CONDITIONS INIT C 51 Example 1 As a first example we show how to prescribe a fixed inflow bo
16. this has nothing to do with parallel decomposition which is separately carried out by MPI If say a uniform grid covers the whole physical domain this number should be set to 1 If two consec utive adjacent grids are used then 2 is the correct choice and so on For each patch the triplet double integer char specifies respectively the leftmost node coordinate value number of points and grid type for that patch there must be as many triplets as the number of patches Since patches do not overlap the rightmost node value of one grid defines the leftmost node value of the next adjacent one The last double specify the rightmost node coordinate value of the last segment which is also the rightmost node value in that direction If a dimension is ignored then 1 grid point only should be assigned to that grid The global domain therefore extends in each direction from the first double node coordinate to the last double node coordinate These values can be accessed from anywhere in the code using the global variables g domBeg d and g domEnd d where d IDIR JDIR KDIRis used to select the direction The grid type char entry can take the following values e vor uniform A uniform grid patch is constructed if xz and xp are the leftmost and rightmost point of the patch the grid spacing becomes LR TIL Az ON e sor stretched a stretched grid is generated Stretched grids can be used only if at least
17. tronic recombination as well as charge transfer with H and He processes see IT MMO8 Similarly the energy loss term is A nanat T X Lee Lir AT X D XpLe ne T Be 9 14 k where By is the fractional abundance of the element and j lt i is the total cooling for one ion specie that is computed and saved to external files by the tables generation program then loaded at runtime In Eq 9 14 Lp and Li r represent the energy losses in bremsstrahlung and ionization recombi nation processes respectively nat and ne are the total atom and electron number densities respectively MINEg uses a dynamically switching integration algorithm for the ion species and energy designed to maximize the accuracy while keeping the computational cost as low as possible 10 Additional Modules 10 1 The ShearingBox Module Figure 10 1 Schematic representation of the shearing boundary condition The computational domain central box is assumed to be surrounded by identical boxes slid ing with constant velocity w qQoLxz with respect to one another The shearingbox provides a local model of a differentially rotating system obtained by expanding the tidal forces in a reference frame co rotating with the disk at some fiducial radius Ry The validity of the approximation and of the module itself is restricted to a Cartesian box considered small with respect to the global flow with a steady flow consisting of a lin
18. you must invoke the python script with the with sb and the with fargo options In this case the macro FARGO_AVERAGE VELOCITY 10 2 1 is automatically turned to NO so that the background orbital velocity is prescribed analytically with the FARGO _SetVelocit y function With FARGO however the source terms in the momentum and energy equations are slightly different IMFS 12 SG10 Sm ETAK a 2 Qoprr J pz k Se pvyv2 ByBz qQo pvz 022z One can see that radial gravity disappears and therefore only the vertical component of gravity must be included in the BodyForceVector or BodyForcePotential functions The additional term in the energy equation represents the work done by Reynolds and magnetic stresses because of the radial shear This term is accounted separately for during the FARGO transport step CHAPTER 10 ADDITIONAL MODULES 88 10 2 The FARGO Module The FARGO MHD module permits larger time steps to be taken in those computations where a grid aligned supersonic or super fast dominant background orbital motion exists see MFS 12 The algorithm decomposes the total velocity into an average azimuthal contribution and a residual term v v w 10 4 where v is called the residual velocity while w is a background velocity field that must be solenoidal The MHD or HD equations are solved in two steps i a linear transport operator corresponding to the velocity w in the direction
19. 1 INTRODUCTION 14 Option Description with chombo enables support for adaptive mesh refinement AMR using the Chombo library Chapter 12 WIth E0 enables support for finite difference schemes with targo enables support for the FARGO MHD module S with sb enables support for the shearing box module 410 1 disables the curses terminal control feature of the Python script Instead a shell based setup will lt O CurSses be used This switch can be used to circumvent problems with the ncurses library present on some systems e g Snow Leopard 10 6 Table 1 2 Command line options available when running the Python setup script 1 3 Configuring PLUTO In order to configure and setup PLUTO for a particular problem four main steps have to be followed the resulting configuration will then be stored in 4 different files part of your local working directory e definitions h header file containing all problem dependent preprocessor directives required at compilation time physics module geometry dimensions etc This is the subject of Chapter e makefile needed to compile PLUTO and it depends on your system architecture This is described in Chapter f e pluto ini startup initialization file containing run time parameters grid size CFL This is the subject of Chapter 4 e init c implements initial boundary conditions etc See Chapter 5 The Python script setup py is used for the first two steps while the re
20. 10 23 gr cm 3 or 10 03 50 0 means 50 1 Km sec x 10 0 means 10 x 1 Km sec x define pressure to that sound speed 1 1 e 6 10 Km sec x v PRS v RHO cs cs g_gamma x means 100 gmm Xrho_0x v_0 2 x where CONST_PI 7 is another pre defined constant With this initialization the sound speed is exactly cs 10 Km s 5 1 2 Specifying Temperature and Gas Composition In many applications it may be more convenient to use a reference temperature to initialize pressure or velocity In PLUTO a direct relation between pressure and density in code or non dimensional units and temperature in Kelvin is provided by 2 where kg is the Boltzmann constant m is the atomic mass unit u is the mean molecular weight while p and p are in code i e non dimensional units The conversion factor K depends on UNIT_VELOCITY and it is provided by the macro KELVIN Eq can be easily used to determine pressure once the temperature is known for instance 0 5 x Density in code units 1 63 x Temperature in KELVIN x MeanMolecularWeight v RHO T KELVIN mu x Obtain pressure in code units x where y is the mean molecular weight P PMtotMu 5 3 while niot is the total number density of particle and it depends in general on the composition of the gas The mean molecular weight may be computed e by calling the MeanMolecularWeight function when the gas comp
21. Andrea PLUTO Test_Problems MHD Shock_Cloud Inffj lt lt pload m Sim time 6 000000e 02 Grid size 256 x 256 File mode single file Variables 8 1 gt rho 2 gt vxl 3 gt vx2 4 gt vx3 5 gt bx1 6 gt bx2 7 gt bx3 8 gt prs 0 0 Figure 11 5 An example of visualization of a binary datafile with mathematica PLUTO data files can be displayed with Mathematicd using a notebook interface to create an in teractive document A simple reader interface is provided by Tools Mathematica pload m and it can be launched to load and display binary datafiles written in single or double precision flt and dbl Data is stored into lists and can be handled using a variety of built in functions in Mathematica Grid and time information are also read from the out log files and stored into the variables nx and ny number of points t current time level nvar number of variables dt current time step A typical interactive session once you open an empty book is AppendTo Path ToFileName home mignone PLUTO Tools Mathematica SetDirectory home mignone PLUTO Test_Problems MHD Shock_Cloud lt lt pload m Please remember to type Shift Enter after each line to make the Wolfram Language process your input The first line simply adds the Tools Mathematica directory to the path the second line changes directory to the working location and the third invokes the reader Once executed pload m re
22. Default is NO Used in the RMHD module 46 4 to evolve the total energy minus the mass density contribution see MMO7 This is more advisable in order to avoid numerical errors in the non relativistic limit and catastrophic cancellation problems De fault is YES Shearing box module only 410 1 The value of the orbital fre quency parameter Qo see 410 1 The default value is 1 0 Shearing box module only 410 1 The value of the differential shear parameter q see The default value is 1 5 proper for a Keplerian profile Shearing box module only 10 1 Symmetrize the hydrody namical fluxes at the left and right x boundaries in order to en force conservation of hydrodynamic variables like density mo mentum and energy no magnetic field Default is YES Shearing box module only g10 1 Symmetrize the y component of the electric field at the left and right x boundaries to enforce conservation of magnetic field only in 3D see Sr c MHD Shearing_Box sb_fluxes c Default value if YES Shearing box module only 410 1 Symmetrize the z component of electric field at the left and right x boundaries to enforce conservation of magnetic field Default is YES Print during the integration log the three time steps hyper bolic parabolic and radiative from which the CFL condition is estimated Sets the value of the y parameter used to control the effi ciency of Super Time Stepping integration for parabolic
23. Figure 12 2 Density distribution with overplotted AMR levels for the disk planet interaction 12 4 3 Visualization with pyPLUTO The simulation data obtained from PLUTO Chombo is written as a HDF5 file which can now be visu alised and analysed using pyPLUTO 11 3 3 The reader for HDF5 files with AMR data in pyPLUTO has been developed by Dr Antoine Strugarek Departement of Physics University of Montreal and has same capablities as that of IDL s HDF5LOAD In order to use this reader it is required to install h5py package the Pythonic interface to HDF5 data The syntax you need is similar to that used for static grids For example in order to read data 0001 hdf5 at the equivalent resolution provided by the 4 refinement level gt ipython pylab In 1 import pyPLUTO as pp In 2 D pp pload 1 datatype hdf5 level 4 Now the pyPLUTO pload object D has all information regarding the data To visualize say the density p one can use the pyPLUTO Image class as follows In 3 I pp Image In 4 I pldisplay D D rho xl D xl x2 D x2 cbar True vertical To plot 2D R Phi data obtained from a POLAR AMR Grid In 5 1 pldisplay D D rho x1 D xl x2 D x2 cbar True vertical polar True True To plot 2D R Theta Slice from 3D POLAR AMR Data In 6 I pldisplay D D rho D n3 2 xl D x1 x2 D x2 cbar True vertical polar True False Further AMR levels in form of
24. MHD equations are discretized using the following di vergence form Op a AV pv Om Op B pui Bs At V mv BB 2 i a gt E OMe 1 Op E 109 pugur BoB pus B ae mav BaB p ij DD A OU r OMe or 1 Op B 1 Ot E a aa rsin9 E e r sin 0 75 O 5 B p V E pi p0 0 B v B oB pa 1 sin OE o 1 DES 9 Ot r sin 0 00 rsin do 7 o Bo 1 dE 1 rEg _ 6 Ot rsin0 0p r Or E OBS 1 0 rEg 10 0 Ot r Or r 00 E A 4 Note that curvature terms are present in the radial and meridional components while the azimuthal component is discretized in angular momentum conserving form The corresponding divergence oper ators are 1 O r F 1 a sin F 1 OFs ve r2 Or r sin 0 00 rsind 06 A5 UF 1 O F 1 O sin 0F 1 OF r3 Or rsin 0 00 rsin0 Od In the previous equations v v va v and B B Bo B are the velocity and magnetic field vectors Eo Eg are the components of the electromotive force v x B gis the body force vector and is the gravitational potential 116 A 2 Special Relativistic MHD Equations A 2 1 Cartesian Coordinates In Cartesian coordinates x y z the relativistic MHD equations 6 11 take the form OD dl D Di V Dv 0 OMy 2 Opt n J EV w 0 vev bsb 57 00 m Op ae V w b uyv byb By Gy Ot Oz A 6 la Dvu g
25. Ot OB Oly Oly 7 o Du o OBy Of 0 7 Ot il Oz Ox 3 oB 0 7 dE _ 9 Ot o oy where D yp is the lab density m w b v y v B b is the momentum density w is the gas enthalpy b B y v BY v vz vy vz is the velocity B B By B2 is the magnetic field in the lab frame b B y y v B v is the covariant field Ez Ey Ez are the components of the electromotive force v x B and gis the body force vector A 2 2 Polar Coordinates In polar cylindrical coordinates R z the RMHD Equations 6 11 are discretized using the following form V Dv R yy w B upv brb e Feet VE w P vow bob ra V w b v v b b a EV m Dv ld ee E dt ROH Oz OB Er OE t az OR OB 10 REs 10 ot R OR R Oo MEE Bo Bo pn Me Beato 2 PI pgz A 7 Dv 9 Note that curvature terms are present in the radial component while the azimuthal component is dis cretized in angular momentum conserving form The corresponding divergence operators are 10 RFr 10F OF R OR ROS Oz 1 O RFFR 10F oF R2 OR R 0z VF A 8 vVi F In the previous equations v vr vg vz and B Br Bg B are the velocity and magnetic field vectors Er Ep Ez are the components of the electromotive force E v x B gis the body force vector and is the gravitational potential A 2 3 Spherical Coordinates In s
26. PATH variable to the bin folder where the executable GUI_pyPLUTO py exists after the installation of source code see installation notes in Doc pyPLUTO html and then apply the following commands in the data directory gt GUI_pyPLUTO py default is for dbl files gt GUI_pyPLUTO py float for flt files gt GUI_pyPLUTO py vtk for vtk files Along with the code an example folder with some sample py files are provided for certain test prob lems The source codes from these files along with their outputs are listed in the HTML documentation It is required to first run the respective test problem and generate the data files after which the user can run the sample py files as follows from the Tools pyPLUTO examples folder gt python samplefile py where the samplefile py are listed in CHAPTER 11 OUTPUT AND VISUALIZATION Table 11 3 List of sample py files provided in the Tools pyPLUTO examples folder samplefile py sod py Rayleigh taylor py stellar wind py jet py orszag tang py Sph_disk py flow_past_cyc py Test Problem HD Sod HD Rayleigh_Taylor HD Stellar_Wind MHD Jet MHD Orszag_Tang MHD FARGO Spherical_Disk HD Viscosity Flow_Past_Cylinder 103 CHAPTER 11 OUTPUT AND VISUALIZATION 104 11 3 4 Visualization with Mathematica le A Ls DA ew ESSE e A h 3 SetDirectory C cygwin64 home Andrea PLUTO Test_Problems MHD Shock Cloud i El Untitled 1 Out 3J C cygwin64 home
27. Table 9 1 The number density of various hydrogen forms are determined by solving the chemical rate equations which have a general form as dn dt y kj NN Mi y kiat 9 10 j k J where n is the number density k is the rate of formation of i specie from all j and k species and k jis the rate of destruction of it specie due to all j species The code integrates the three hydrogen fractions defined above using the advection equation of the form OX Ot where the first term on the rhs is treated during the hydro step while only the second is integrated during the cooling step The source terms S is essentially the difference between formation and destruction rate of a particular specie see eq 9 10 Additionally the internal energy losses take into accounts various hydrogen cooling processes v VX Si 9 11 A Aci Arr Arotvib Anadiss Aerm 9 12 where Ac and Agp are losses due to collisional ionization and radiative recombination respectively The remaining terms Arotvib AHadiss and Agrain are associated with molecular hydrogen and represent CHAPTER 9 OPTICALLY THIN COOLING Araib Anodis Ag ia ARR n 10 cm Xj 0 835 Xz 0 0823 Xj 0 0004 T WY U o L Uv e lt y E lt Temperature K 84 Figure 9 1 Variation of the radia tive cooling functions A with tem perature due to various processes that can affect the total energy of a gas co
28. The RMHD module solves the following system of conservation laws D Dv dl o m wey vv bb lp al y ewe 0 6 11 B vB Bv where D is the laboratory density m is the momentum density E is the total energy including contri bution from the rest mass ns yu B D YP b B y q v Bw m wv bb wi ph B Y v BY E wey b b p B v B Pt AS a Notice that the components of the momentum tensor may also be written as BI a MY wuu bb m u m u 5 v B B Y Primitive variables are similar to the RHD module but they also contain the magnetic field V p v p B The quasi linear form of the RMHD is not available yet and algorithms using the char acteristic decomposition of the equations or the quasi linear form are not available Therefore the CHARACTERISTIC_TRACING step cannot be used and the HANCOCK scheme works by default using the conservative predictor step rather than the primitive one On the other hand Runge Kutta type integrators works well for the RMHD module The available equations of state are IDEAL and TAUB already introduced for the RHD module see for the extension of this EOS to the RMHD equations The RMHD sub menu offers some of the switches already discussed in the MHD module 6 2 or in the RHD 46 3 module Divergence control is achieved using the same algorithms introduced for MHD namely 8 wave 36 2 3 1 divergence cleaning 36 2 3 2 and
29. and their order is specified by the elements of the array input_ var until the value 1 is encountered You may provide only some of the variables used by PLUTO and not necessarily all of them The number of elements per variable should exactly match the number of grid points defined by the input grid Finally the function InputDatalnterpolate is used to map the values of the variables con tained in the input binary data on the grid employed by PLUTO using bi or tri linear interpolation at the desired coordinate location InpPutDatalnterpolate Y Xle X2p X3 J where v is the same array of primitive variables used in the Init function and x1 x2 x3 are the coordinates at which interpolates are required In the example below density and velocity components are assigned from the input binary file tmp data 0010 dbl defined on the computational domain specified in tmp grid out void Init double v double xl double x2 double x3 static int first call 1 if f1r8t_ call int k input_var 256 for k 0 k lt 256 k input_var k 1 input_var 0 RHO input_var 1 VX1 input_var 2 VX2 input var S VX3 input_var 4 InputDataSet tmp grid out input_ var InputDataRead tmp data 0010 db1 gt first_call 0 CHAPTER 5 INITIAL AND BOUNDARY CONDITIONS INIT C 48 InputDataInterpolate v xl x2 x3 Beware that interpolation is performed only on the variables specified by
30. are used third column CHAPTER 6 BASIC PHYSICS MODULES 60 6 2 3 Controlling the V B 0 Condition 6 2 3 1 Eight Wave Formulation In the eight wave formalism Pow94 PRL 99 magnetic fields have a cell centered representation Ad ditional source terms are added on the right hand side of Eqns 6 4 p 0 m B t E AA v B B V Contributions to V B are taken direction by direction Note that the 8 wave method keeps V B 0 only at the truncation level and NOT to machine accuracy More accurate treatments of the solenoidal condition can be achieved using the other two formulations The 8 wave algorithm should be used in conjunction any Riemann solver with the exception of h11d 6 2 3 2 Hyperbolic Divergence Cleaning In DKK t02 see also MT10 MTB10 for additional discussion the divergence free constraint is en forced by solving a modified system of conservation laws where the induction equation is coupled to a generalized Lagrange multiplier GLM Using the mixed GLM hyperbolic parabolic correction the induction equation and the solenoidal constraint are replaced respectively by OB 7 Op o G tV wB Bv Vy 0 a AS 6 8 where cp CFL x Almin At is maximum speed compatible with the step size cp Almincn a and Almin is the minimum cell length The free parameter a controls the rate at which monopole are damped IMT10 and its value is set by the user defined constant GLM_ALPHA default value 0 1
31. as T and vortz in the example above have not been allocated into memory In order to change the default attributes follow the example in the next subsection 11 2 1 Changing Attributes Defaults attributes which variables in which output have to be written image attributes can be easily changed through the function ChangeDumpVar located in the file Src userdef_output c To include exclude a variable from a certain output type use SetDumpVar var type YES NO Here var is a string containing the name of a variable listed in Table or an additional one de fined in your pluto ini The type argument can take any value among DBL_OUTPUT FLT_OUTPUT CHAPTER 11 OUTPUT AND VISUALIZATION 97 VTK_OUTPUT DBL_H5_OUTPUT FLT_H5 OUTPUT TAB_OUTPUT PPM_OUTPUT PNG_OUTPUT This is a sketch of how this function may be used void ChangeDumpVar Image image x a pointer to an image structure SetDumpVar bx1 FLT_OUTPUT NO SetDumpVar pes PPM QUIPUT YES SetDumpVar vortz PNG_OUTPUT YES image GetImage rho image gt slice_plane X13_PLANE image gt slice_coord 1 1 image gt max image gt min 0 0 image gt logscale 1 image gt colormap red In this example the variable bx1 is excluded from the flt output prs and vort z defined in the previous example are added to the ppm and png outputs respectively Furthermore the default image attributes of
32. case for example when the 2 5 D formalism is used where integration is performed along the first two coordinates say x y but the fluid has a non vanishing velocity component along the third 17 CHAPTER 2 PROBLEM HEADER FILE DEFINITIONS H 18 direction as well say 0v 0x 0v 0y 0 An example is an axisymmetric 2 D cylindrical problem such as a disk or a torus in the r z plane with a uniform rotation in the azimuthal direction where it is assumed 0 04 0 In all cases it is required that DIMENSIONS lt COMPONENTS 2 1 3 GEOMETRY Sets the geometry of the problem Spatial coordinates are generically labeled with x1 x2 and x3 and their physical meaning depends on the value assigned to GEOMETRY e CARTESIAN Cartesian coordinates 11 12 x3 x y z e CYLINDRICAL cylindrical axisymmetric coordinates 11 x2 r z 1 or 2 dimensions e POLAR polar cylindrical coordinates x1 12 x3 r 0 z e SPHERICAL spherical coordinates 21 12 73 r 0 Note that when DIMENSIONS 2 the third coordinate x3 is meaningless and will be set to zero similarly in 1 D x2 and x3 do not play any role Whenever present however the component of vectors both in spherical and cylindrical coordinates is integrated by discretizing the equations in angular momentum conserving form We warn that non Cartesian geometries are handled better when a multi stage unsplit integrator i e Runge Kutta is used es
33. containing d gt Vc nv k j i a four index array of primitive variables defined at the cell center The integer nv RHO VX1 NVAR 1 labels the variable see Table 5 1 while k j and i are the zone indices of the x3 x2 and x direction note the reversed order respectively d gt Vs nv k j i staggered MHD only a four index array containing the three com ponents of the staggered magnetic field BX1s BX2s BX3s if any defined at zone faces see Fig 5 1 These components only exists in the MHD or RMHD modules when using the Con strained Transport algorithm to control the V B 0 condition see 6 2 3 3 for more details Important Face centered staggered magnetic fields and cell centered fluid variables are defined on different zone stencils see Figure 5 1 The zone centering and the corresponding index range is encoded in the box structure see below All variables must be assigned at a user defined boundary with the exception of the staggered component of magnetic field normal to the interface if you are using the Constrained Transport CT method see 46 2 3 3 e xbox a pointer to a RBox structure defining the rectangular portion of the domain over which ghost zone values should be assigned Since cell centered and face centered data are defined on different box structures its usage is maily intended to discriminate between cell centered variables and face centered variables using the structur
34. diffu sion terms see Chapter 8jand If not set the default value is 0 01 122 Table B 1 Continued from previous page Name TC_SATURATED_F LUX T CUT RHOE TV_ENERGY_TABLE TV_ENERGY_TABLE_NX TV_ENERGY_TABLE_NY UNIT_DENSITY UNIT_LENGTH UNIT_VELOCITY VTK_TIME_INFO VIK VECTOR DUMP Value YES NO real YES NO int int real real real YES NO YES NO Description Include saturation effects when computing the thermal conduc tion flux Default value is YES Sets the cut off temperature in K used in the PVTE_LAW equa tion of state 87 3 Zones with temperature below T_CUT_RHOE will be reset to this value and the internal energy will be rede fined accordingly Default value is 10 K Used for the PVTE_LAW EOS in ionization equilibrium When set to YES replaces function evaluations of the caloric EOS internal energy and its inverse e T p and T e p with lookup table and bilinear interpolation This results in a consid erably faster execution Default is YES Sets the number of x points used to construct the tempera ture table for the PVTE_LAW EOS Default value is set in Sr c EOS PVTE internal_energy c Sets the number of y points used to construct the tempera ture table for the PVTE_LAW EOS Default value is set in Sr c EOS PVTE internal_energy c Sets the unit density in gr cm Default value is the proton mass per cm Sets the unit le
35. every fixed number of steps When set to NO w is computed from the FARGO SetVelocity function Default is NO or YES if FARGO is used with the shearing box module Sets how often the orbital velocity should be recomputed in the FARGO transport step Default is 10 Sets the spatial order of the reconstruction used during the lin ear transport step of the FARGO algorithm The allowed values are 3 third order PPM like reconstruction or 2 second order MUSCL HANCOCK scheme Default is 3 Sets the value of the constant a used monopole damping in the GLM formalism see 46 2 3 2 The default value is 0 1 Enable the E xtented GLM form of the MHD equations see q6 2 3 2 Default value is NO Hydrogen mass fraction X muny A p Used to compute FRAC_He and FRAC_Z in the definition of the mean molecular weight see 45 1 2 Default value is X 0 7110 Helium mass fraction Y munHeAne p Used to compute FRAC_He and FRAC_Z in the definition of the mean molecular weight see Note that the fraction of metals is always computed as Z 1 X Y The default value is Y 0 2741 121 Table B 1 Continued from previous page Name PPM_ORDER PV _ TEMPERATURE TABLE PV_TEMPERATURE_TABLE_NX PV_TEMPERATURE_TABLE_NY RECONSTRUCT_4VEL RMHD_FAST_EIGENVALUES RMHD_REDUCED_ENERGY SB_OMEGA SB_O SB_SYMMERIZE_HYDRO SB_ SYMMERTZE_EY SB_ SYMMERTZE_EZ SHOW_TIME_STEPS STS_NU Continued on next page Value
36. example shown in Fig shows how to visualize magnetic pressure and density in two different windows and overplot the velocity field For more details consult Doc idl_tools html CHAPTER 11 OUTPUT AND VISUALIZATION Magnetic Field Density a E e E E E E gt E ff F F E 5 E be E E Gi E F z i JE E H 5 E E 5 E E H 13 E E E jE le H E EL O Se SS woe oes ee N 4 D z 4 5 6 IDL gt IDL gt pload 1 silent IDL gt loadct 5 LOADCT Loading table STD GAMMA II NIDL gt display x1 x1 x2 x2 vbar title Density rho IDL gt display x1 x1 x2 x2 vbar title Magnetic Field bx172 bx2 2 nwin 1 IDL gt vecfield bx1 bx2 x1 x2 overplot col 120 IDL gt 101 Figure 11 3 An exam ple of visualiza tion in IDL using the display pro routine CHAPTER 11 OUTPUT AND VISUALIZATION 102 11 3 3 Visualization with pyPLUTO a ma X pyPLUTO Save As Quit Nstep Load data x x3 v Log v Polar v Aspect Contour v Arrows Contours Arrows Labels XRange YRange VarRange Information Dir fimport sunstar4 phybva PLUTO Test_Problem HD Disk_Planet Domain nl xn2xn3 128 x 16 x 384 A 4 t t 1 y A Time Status Nlast 10 Loaded Nstep 10 Present Time 1 0 cu y 4 4 e a A ANN xn ax a A A A A A am A A S e e we ea a e e Ms y s ne OS OSO Figure 11 4 An example of visualization with the pyPLUTO tool Binary data files dbl flt a
37. for both software packages TT hoOShD s sste asada 4 a p 322 oe BB F AAR SR Re KD Damos ania OEC JARE Aak Eea F i e AN DB data Density Appl Y Apply operators to alli 15 59644 sto Y Apply subset selection Opaci 10 El 0 683735 Figure 11 6 An example of visualization of an xmf h5 data file using VisIt left or ParaView right Visualization of HDF5 files Both Visilt and Paraview interpret the cell centered grid and data con tained in the Pixie files as node centered as a consequence the first and the last half cells in every direction are clipped from the images e g a small sector around 0 is chopped from a periodic polar plot covering the 27 angle Therefore for every h5 file PLUTO writes also a xmf text file in XDMF format that describes the content of the corresponding HDF5 file The xmf files can be directly opened by Visit and ParaView so as to provide the correct data centering and avoid the image clipping Besides we noticed that the Pixie reader crashes e g using ParaView 3 14 3 98 or incorrectly reads the h5 files e g using VisIt versions after 2 6 while all versions of both VisIt and Paraview properly open the xmf files All the variables are read as scalar quantities Visualization of VTK files PLUTO writes vtk files using a cell centered attribute rather than point centered as in p
38. gcc defs and press enter All the information relevant to the specific problem should now be stored in the four files init c assigns initial condition and user supplied boundary conditions pluto ini sets the number of grid zones Riemann solver output frequency etc definitions h specifies the geometry number of dimensions interpolation time stepping scheme and so forth and the makefile 5 Exit from the main menu Quit or press q and type PLUTO Test_Problems HD Sod gt make to compile the code 6 You can now run the code by typing CHAPTER 0 QUICK START 7 PLUTO Test_Problems HD Sod gt pluto At this point PLUTO reads the initialization file pluto ini and starts integrating The run should take a few seconds or less and the integration log should be dumped to screen Data can be displayed in a number of different ways If you have for example Gnuplot v 4 2 or higher you can display the density output from the last written file using gnuplot gt plot data 0001 db1 bin array 400 400 400 form double ind O where ind 0 1 2 may be used to select density velocity or pressure If you have IDL installed on your system you can easily plot the density by IDL gt pload 1 IDL gt plot x1 rho The IDL procedure pload is provided along with the code distribution 0 3 Running the Orszag Tang MHD vortex test 1 Change directory to PLUTO Test Problems MHD Orszag Tang 2 Choose a configura
39. identify the parameter labels i e the corresponding indices of the array g inputParam while the double values on the right are the actual user defined parameter values The number of parameters specified in this section must exactly match the number and the order given in definitions h 5 Initial and Boundary Conditions Init c The source file init c provides a set of functions that are used to define set and configure your own specific problem These include e Init sets initial conditions as functions of the spatial coordinates z1 2 3 e UserDefBoundary sets user defined boundary conditions at the physical boundary sides of your computational domain if necessary e Analysis run time data analysis and reduction e BodyForceVector BodyForcePotential defines the vector components of the acceler ation vector and or the gravitational potential e BackgroundField sets a background force free magnetic field The init c must be part of your local working directory and a template can be found in Src Templates init c Functions are described in the next sections 5 1 Inital Conditions the Init function The Init function is used to assign the initial condition as a function of the spatial coordinates Syntax vozd Init double v double xl double 22 double x2 Arguments e v a pointer to a vector of primitive quantities A particular variable is located by means of an index p v RHO
40. init c or any other source file for conditional compilation Hif SETUP_VERSION 0 x implements version 0 elif SETUP_VERSION 1 x implements version 1 x endif 2 Advanced options expert users override the default value of special constant macros used through out the code a comprehensive list of which is given in Appendix B 3 This gives the user more control on the code and it avoids copying and modifying source files in the local working directory A simple example is given by configuration 04 in the Test_Problems MHD Rayleigh_Taylor problem where the symbolic constants x Beg user defined constants do not change this line x define USE_RANDOM_PERTURBATION NO define CHOMBO_REF_VAR RHO x End user defined constants do not change this line x are be used in init c to enable disable a random perturbation and to tell PLUTO Chombo that density must be used as the refinement variable Another typical example is provided by the CGS physical units described 45 1 1 CHAPTER 2 PROBLEM HEADER FILE DEFINITIONS H 27 2 4 Frequently Used Options Besides the user defined constants discussed so far which are accessible via the Python script defini tions h contains additional switches that are frequently used and are always shown These constants cannot be modified using the Python script and must be adjusted manually by editing your definitions h INITIAL SMOOTHING YES NO
41. is invoked with the HDF5 keyword For instance in order to read data 0001 hdf5 at the equivalent resolution provided by the 4 refinement level you need IDL gt pload hdf5 2 level 4 will load data 0002 hdf5 ref level 4 Note that IDL re interpolates the required dataset to a uniform mesh with resolution given by the spec ified refinement level As an example we show how to visualize the density map for the relativistic Kelvin Helmholtz test problem as in Fig 12 1 corresponding to the output of configuration 03 of Test_Problems RMHD KH Density map Figure 12 1 Density maps of the relativistic Kelvin Helmholtz test problem at t 5 Refinement levels are displayed using the oplotbox routine The figure has been produced with the following IDL commands IDL gt PLOAD HDF5 dir DATA_03 5 lev 4 x2range 0 5 1 5 IDL gt LOADCT 6 IDL gt DISPLAY x1 x1 x2 x2 vbar rho imax 1 1 imin 0 65 title Density map IDL gt OPLOTBOX ctab 3 The last command OPLOTBOX overplots the levels of refinement utilizing the color table 3 If you are plotting a 2D map in curvilinear coordinates polar or spherical using the DISPLAY procedure setting the POLAR keyword you can use the same POLAR keyword for the OPLOTBOX procedure to correctly overplot the levels of refinement Reading Large Datasets It may occur that the dataset one wishes to load exceeds the available mem ory In that case it is useful to load only a por
42. it is possible to suppress saturation effects by turning the macro TC_SATURATED_F LUX to NO see also Appendix B 3 in this case Fe Fass The coefficients appearing in Eq 8 8 and in the definition of the saturated flux may be spec ified using the function TC_kappa in your local copy of PLUTO Src Thermal_Conduction tc_kappa c and by noting the equivalence ky gt kpar k gt knor and x phi The variable knor can be ignored in the HD case where k xy Proper setting of units and dimensions is briefly discussed in 83 1 The thermal conduction module is implemented inside Src Thermal_Conduction and works in 1 2 and 3 dimensions in all systems of coordinates Derivative terms are discretized at cell interfaces using second order accurate finite differences and assuming a uniform grid spacing Integration may proceed via standard explicit time stepping or Super Time Stepping see Note Thermal conduction behave like a purely parabolic diffusion operator in the classical limit oo and like a hyperbolic operator in the saturated limit VT gt oo Thus in the general case a mixed treatment is required where the parabolic term is discretized using standard central differences and the saturated term follows an upwind rule BTH08 MZT 12 In this case and when Super Time Stepping integration is used to evolve the equations sev eral numerical tests have shown that problem involving strong discontinuities may req
43. macro introduced in automatically implements the correct loop over the boundary ghost zones Thus at the x boundary for instance one needs to assign A at Ll Lo 4 1 L3 k for i 0 IBEG 1 3 1 j k 5 at Virta ga a ktl The component normal to the interface b in this case should NOT be assigned since it is automatically computed by PLUTO from the V B 0 condition after the tangential components have been set Example The following example prescribes user defined boundary conditions at the lower x2 boundary for a MHD jet problem in cylindrical coordinates x R x2 z Inflow conditions are given as p UR vz p Br Bz 1 0 10 1 T 0 3 for R lt 1 while a symmetric counter jet is assumed for R gt 1 1f side X2_BEG JETVAL vjet x beam jet values x R grid DIR cylindrical radius if box gt vpos CENTER select cell centered varaibles only x BOX_LOOP box k j i loop on boundary zones x for nv 0 nv lt NVAR nv vout nv d gt Vc nv k 2 xJBEG jJ 1 i vout VX2 1 0 if PHYSICS MHD vout BX1 1 0 endif Lor ny NVAR nv smooth out the two solutions x 0 nv lt d gt Vc nv k jJ i vout nv vjet nv vout nv Profile R i nv else if box gt vpos X1FACE 1 select xl staggered component x Hifdef STAGGERED_MHD Rp grid IDIR A x right interface area
44. now set the names of the 3 auxiliary parameters previously introduced To do so use the arrow keys to select each of them and explicitly write their names P_IN POUT and GAMMA and press enter to confirm 4 Finally we complete the python session by setting the architecture for the makefile In the makefile menu choose your system configuration e g Linux gcc defs for Linux Press enter to confirm MN YY You are now done with the Python script and can exit by pressing either q or selecting quit At this point you should find the following four files inside your Blastwave folder definitions h init c makefile pluto ini sysconf out Next we need to edit the two files pluto ini and init c The first one defines the computational domain and certain properties of the run i e time of integration first timestep etc The second one sets the initial conditions for the blast wave problem a circular region of high pressure in a lower pressure ambient Edit pluto ini to make the following changes e The domain should span from 1 to 1 in both dimensions with 200 points in each direction Xl grid 1 LD 200 u 1 Xx2 gGgrld 1 LO 200 u 1 20 0 e The simulation should stop when time reaches 0 04 tstop 0 04 with the first timestep being First de LL S 0 Save the files every t 0 004 in double precision and in multiple_files format dbl 0 004 1 multiple_files e At the end of the file set the numerical values for the 3 parameters
45. of orbital motion and ii a standard nonlinear solver applied to the original equations written in terms of the residual velocity v The Courant condition is then computed only from the residual velocity leading to substantially larger time steps In MFS 12 it has been shown that if the characteristic velocity of fluctuations are comparable in magnitude than the expected gain in polar coordinates is roughly 1 M 1 1 m _ Atr isk AR RAG Az At 1 1 a 10 5 e AR RAS Az where Atp and At are the FARGO time step and the standard time step respectively whereas M w A and X uv cfa is the characteristic speed in the e direction The discretization is fully conservative in both angular momentum and total energy The MHD module works only with the Constrained Transport CT method to control divergence free condition 10 2 1 Using the Module The FARGO MHD module is implemented in the directory Src Fargo and can be enabled by invoking the python script with the with fargo option It works in Cartesian polar and spherical coordinates with a dimensionally unsplit time stepping scheme i e with DIMENSIONAL_SPLITTING set to NO The background velocity can be computed by PLUTO in two different ways depending on the value of the macro FARGO_AVERAGE_VELOCITY e YES default the azimuthal velocity v or vg is averaged along the corresponding orbital di rection This operation is performed once every f
46. one uniform grid is present The stretching ratio r is computed as follows Lor z Az r r r rR TL gt r r IR TL l r Az where Az is taken from the closest uniform grid N is the number of points in the stretched grid and xz and xp are the leftmost and rightmost points of the patch e 1 a logarithmic grid is generated When 1 is invoked the mesh size is increasing with the coordinate 1 An zi lx 1 zz 10 1 A y 198 PEE a L CHAPTER 4 RUNTIME INTITALIZATION FILE PLUTO INI 34 when 1 is invoked the mesh size decreases with the coordinate 1 An zia z xp 10 1 A 77 los EAS In practice the mesh spacing in the 1 grid is obtained by reversing the spacing in the 1 grid Note The interval should not include the origin when using a logarithmic grid In CYLINDRICAL or SPHERICAL coordinates a radial logarithmic grid has the advantage of pre serving the cell aspect ratio at any distance from the origin In addition the condition to obtain approximately squared cells aspect ratio 1 is Ar r Ag where Ar rz e 1 is the radial spacing of the first active computational zone This condition can be be used to determine either the number of points in the radial direction or the endpoint 2 Ad 2 Ad log ee N log rL Beware that non uniform grids may introduce extra dissipation in the algorithm Changes in the grid spacing are cor
47. selecting Setup problem see Fig 2 1 If you do not have an existing definitions h a new one will be created for you otherwise the Python script will try to read your current setup from the existing one eoo gt gt Setup problem lt lt PHYSICS DIMENSIONS 2 COMPONENTS 3 GEOMETRY CARTESIAN BODY_FORCE NO COOLING NO RECONSTRUCTION LINEAR TIME_STEPPING RK2 DIMENSTONAL_SPLITTING NO NTRACER o USER_DEF_PARAMETERS 3 Figure 2 1 The Setup prob lem menu needed for your def initions h and makefile creation by moving the arrow keys you should be able to browse through different options The header file definitions h also contains other more advanced switches that are not accessible via the Python script 42 4 and should be changed manually We now describe the options accessible through the Python script 2 1 1 PHYSICS Specifies the fluid equations to be solved The available options are e HD classical hydrodynamics described by the Euler equations e MAD single fluid ideal resistive magnetohydrodynamics e RHD special relativistic hydrodynamics e RMHD special relativistic magnetohydrodynamics 2 1 2 DIMENSIONS amp COMPONENTS DIMENSIONS sets the number of spatial dimensions of your problem whereas COMPONENTS sets the number of vector components such as velocity and magnetic field present in the integration Usually DIMENSIONS COMPONENTS but one can also have more COMPONENTS than DIMENSIONS This is the
48. standard Given suitable hardware this allows the transfer of data out in the user s buffer to proceed concurrently with computation A separate request complete call is needed to complete the I O request i e to confirm that the data has been read or written and that it is safe for the user to reuse the buffer Overall this results in an improved performance for intensive I O computations More details may be found in http Awww prace project eu IMG pdf petascale_astrophysical_simulations_with_pluto pdf Note This is an experimental feature that can be used in the current version of the code only for flt or dbl binary files for saving cell centered data 3 2 HDF5 Library Support If your system is already configured with serial or parallel HDF5 libraries you may enable support for HDF5 I O in the static grid version of PLUTO by turning the makefile variable USE_HDF5 to TRUE inside your defs file If you do not have HDF5 installed you may follow the installation guidelines given in Note that the same HDF5 library can be used for both the static and AMR versions of PLUTO although support for HDF5 in the AMR version is enabled differently see Once USE_HDF5 has been set to TRUE add the HDF5 library path to the list of directories to be searched for header files as well as the corresponding linker option for HDE5 library files Note that different pathnames should be given if you are building PLUTO in serial or parallel mode These info
49. the array input_var The remaining variables if any must still be set inside Init Note When the input geometry differs from the one used by PLUTO vector components are not automatically transformed to the current geometry A configuration example may be found in the Test_Problems HD Blast directory where the initial condition sets an isothermal blast wave propagating in a non uniform density medium The inital den sity distribution is created by the separate file Turbulence c in the same directory and interpolated at runtime by PLUTO using the method outlined above Note Staggered magnetic fields may not be assigned in this way since the divergence free condi tion is not necessarily maintained Using the vector potential components is more advisable CHAPTER 5 INITIAL AND BOUNDARY CONDITIONS INIT C 49 5 2 User defined Boundary Conditions The UserDefBoundary function is used to assign user defined boundary conditions to a particular physical boundary see Fig 5 1 if this has been tagged userdef inside your pluto ini Alternatively this function may also be used to control the solution array at the beginning of every time step inside the computational domain set floor values override the solution etc by first enabling INTERNAL_BOUNDARY to YES inside definitions h see Syntax void UserDefBoundary const Data d RBox xbox int side Grid xgrid Arguments e xd a pointer to the PLUTO data structure
50. the constrained transport 36 2 3 3 Computation of the fast characteristic speeds can be perfomed by replacing the numerical solution of a quartic equation see MMO07 with the analytical solution of an approximate quadratic equtions thus making computation faster This is achieved by setting RMHD_FAST_EIGENVALUES to YES as in DZBLO07 see the Appendix B 3 7 Equation of State In the current implementation PLUTO describes a thermally ideal gas obeying the thermal Equation of State EOS p Mu p nkgT kel 7 1 where p is the pressure n is the total particle number density kg is the Boltzmann constant T is the tem perature p is the density m is the atomic mass unit and y is the mean molecular weight The thermal EOS describes the thermodynamic state of a plasma in terms of its pressure p density p temperature T and chemical composition y Eq is written in CGS physical units Using code units for p and p while leaving temperature in Kelvin the thermal EOS is conveniently re expressed as T ps ga Ki 7 2 Ku p where K is the KELVIN macro which depends explictly on the value of UNIT_VELOCITY Another fundamental quantity is the specific internal energy e whose rate of change under a phys ical process is regulated by the first law of thermodynamics de dQ pd 7 3 where Q represents the heat absorbed or released The internal energy is a state function of the system and can also be related to tempe
51. the rest mass density p three velocity v Vz1 Uz2 z3 and pressure p The relation between U and V is more complicated and is expressed by D py M phYv phyu E ph p where h is the specific enthalpy see Chapter 7 for available equation of states In order to express the equations in primitive quasi linear form one assumes dp c2de where cs is the adiabatic speed of sound Op 1 1 Op Ov 1 v Op 1 ET Tp a oh Ot a 6 10 Op 1 phV v 1 c2 v Vp 0 Ot l 1 vee For more detailed expressions and the characteristic decomposition see MPBO05 Spatial reconstruction may be performed on the four velocity rather than on the three velocity by enabling the macro RECONSTRUCT_4VEL to YES manually in your definitions h see also Appendix B 3 Using the four velocity in place of the three velocity offers in some circumstances the advantage that the total velocity v u V1 u is always less than 1 by construction for any 0 lt u lt oo This is not always the case when the three velocity is used and precautionary measures are used to ensure that v lt 1 CHAPTER 6 BASIC PHYSICS MODULES 66 6 4 The RMHD Module The RMHD module implements the equations of special relativistic magnetohydrodynamics in 1 2 or 3 dimensions Velocities are always assumed to be expressed in units of the speed of light Source and definition files are located inside the Src RMHD directory 6 4 1 Equations
52. to TRUE the same variable name is passed to PLUTO as a define directive with value 1 As an example if USE_HDF 5 is set to TRUE inside a defs file then any C source file containing instruc tions inside a preprocessor directive ifdef USE_HDF5 tendif statement will be compiled Note These switches are effective only in the static grid version of the code and have no effect when creating a PLUTO Chombo makefile 912 2 3 1 MPI Library Parallel Support Parallel executables for the static grid version of PLUTO are created by specifying the name of a MPI C compiler e g mpicc and by setting the makefile variable PARALLEL to TRUE in your defs file mpicc or similar PARALLEL Ht HE HE HE TE FE FE FE HE HE HE FEFE AE FE HE HE FE FE FE FE HE HE E E H H H H H MPI additional spefications HA HH ifeg strip PARALLEL TRUE USE_ASYNC_IO if TRUE enable Asynchronous binary I O endif In this case you may also modify existing variables or add new ones inside the conditional statement beginning with i feq When parallel mode is enabled C source code sections that are specific to MPI should be enclosed inside ifdef PARALLEL endif statements 29 CHAPTER 3 MAKEFILE SELECTION MAKEFILE 30 3 1 1 Asynchrounous I O By enabling the USE_ASYNC_IO to TRUE inside the chosen defs PLUTO allows to replace the standard blocking calls of MPI with non blocking and split collective calls available in the MPI 2 I O
53. vz v VX1 vy v VX2 and so on Although VX1 VX2 and VX3 should be used in any coordinate system in order to avoid confusion an alternative set may be adopted if the geometry is not Cartesian see columns 2 4 in Tablel5 1 Temperature in Kelvin can also be used to initialize density or pressure see 45 1 2 Note PLUTO 3 Users The old array indexing style using DN VX VY and VZ etc used in previous versions of the code is no longer supported e x1 x2 x3 coordinates x1 2 x3 of the computational cell where v is initialized Example 1 The following code sets a disk with radius 1 centered around the origin in a 2D Cartesian domain The flow is stationary and the disk has higher density and pressure p 10 p 30 with respect to the background state p 1 p 1 42 CHAPTER 5 INITIAL AND BOUNDARY CONDITIONS INIT C 43 void Init double xv double x1 double x2 double x3 double r r sqrt xlaxl x2 x2 v VX1 v VX2 0 0 if lt 1 014 v RHO 10 0 v PRS 30 07 jelse v RHO v PRS Tla LD With a small modification the same initial condition can be written in a dimension independent way e g a line in 1D a disk in 2D and a sphere in 3D void Init double xv double xl double x2 double x3 double r2 r Xx2xx2 Xx3 x3 The macro D EXPAND a b c is used to write dimension independent code by conditionally compiling one two or three comma separat
54. 04 0 1 236510 6 step 3 y where step gives the current integration step t is the current integration time dt is the current time step n is the percentage of integration The two numbers in square brackets are respectively the max imum Mach number and maximum number of iterations required by the Riemann solver if iterative e g two shock during the previous step For non iterative Riemann solvers the last number will al ways display 0 The maximum Mach number is a very sensitive function of the numerical method it may be used as a robustness indicator Very large Mach numbers or rapid variations usually indicate problems and or fixes during the computation 1 4 1 Command line options When running PLUTO a number of command line switches can be given to enable or disable certain features at run time Some of them are available only in the static grid version see Table for a description of the available flags CHAPTER 1 INTRODUCTION Option dec n1 n2 i fname hb5restart n makegrid maxsteps n no write no xlpar NO KZpar no xspar restart n show dec x13et x2jet x33Jet xres nl As Description Enable user defined parallel decomposition mode The integers n1 n2 and n3 specify the number of processors along the x1 x2 and x3 directions There must be as many integers as the number of dimen sions and their product must equal the total number of p
55. 1 Ad a 1 RK NO At Ca lt ae oe 2 si 2 o Naim 1 in 2D ENCK CATE NO At max 37 max a 0n ijk Alq ijk NAlG 1 2 in 3D Table 2 1 CFL conditions used by PLUTO for different explicit time stepping methods For a given direction d Alq represents the cell physical length in that direction Ag provides the largest signal speed while Tq accounts for diffusion processes Here HNCK and ChTr stand for HANCOCK and CHARACTERISTIC_TRACING respectively These limits are based on a stability analysis on the constant coefficient advection diffusion equation by by Beckers 1992 Bec92 2 1 8 DIMENSIONAL_SPLITTING Set this feature to YES if you intend to use Strang operator splitting to solve multidimensional equations by a sequence of 1D problem If DIMENSIONAL_SPLITTING is set to NO flux contributions are evaluated from all directions simultaneously Dimensionally unsplit schemes avoid the errors due to operator splitting and are generally preferred Table 2 2 gives a brief description of commonly used setups CHAPTER 2 PROBLEM HEADER FILE DEFINITIONS H 21 RECONST LINEAR PARABOLIC WENO3 LimO3 LINEAR LINEAR PARABOLIC LINEAR TIME ST RK2 RK3 HANCOCK COTE Chir Ghir HANCOCK DIM oP Lh YES NO YES NO YES YES YES NO Cost 2Ndim 3IN4im Naim Nata Naim 2N dim Comments Default setup Compatible with almost every algorithms of th
56. 15 lt lt Working dir Users mignone PLUTO Test Problems MHD Shock Cloud PLUTO dir Users mianone PLUTO roblem Change makefile Auto update Save Setup Quit Figure 1 2 Python script main menu 1 4 Compiling amp Running the Code After the four basic configuration files init c definitions h makefile and pluto ini have been created PLUTO can be compiled from your local working directory by typing MyWorkDir gt make gmake is also fine It is important to remember that the makefile created by Python Chapter 3 guarantees that your work ing directory is always searched before PLUTO Src This turns out to be useful when modifying PLUTO source files 1 5 If compilation is successful type MyWorkDir gt pluto flags for a single processor run or MyWorkDir gt mpirun pluto args for a parallel run are options given to MPI such as number of processors etc while args are command line options specific to PLUTO see Table 1 3 For example MyWorkDir gt pluto restart 5 maxsteps 840 will restart from the 5 th double precision output file and stop computation after 840 steps During execution the integration log will look something like step 0 t 0 0000e 00 dt 1 0000e 04 O 0 000000 0 step 1 t 1 0000e 04 dt 1 0000e 04 O 1 236510 10 step 2 t 2 0000e 04 dt 1 1000e 04 0 1 236510 7 t 3 1000e 04 dt 1 1000e
57. A number of tests suggests that the optimal range can be found for 0 05 S a 0 3 In the mixed formulation divergence errors are transported to the domain boundaries with the maximal admissible speed and are damped at the same time By default y is set to zero in the initial and boundary conditions but the user is free to change it at a user defined boundary by prescribing d gt Vc PSI_GLM k j i inside UserDefBoundary which has the usual cell centered representation The scalar multiplier is not written to disk except for the double format needed for restart The advantage of this formulation GLM MHD is that the equations retain a conservative form no source terms are introduced all variables including magnetic fields retain a cell centered representa tion and standard 7 wave Riemann solvers with a single value of the normal component of magnetic field may be used A slightly different version called the extended GLM formulation that breaks momentum and en ergy conservation has been found to be more robust in problems involving strongly magnetized media see for example configuration 11 in Test_Problems MHD Blast The extended form of the equations can be enabled by adding the user defined constant GLM_EXTENDED to definitions h and setting its value to YES from the Python script For a complete description of the GLM and Extended GLM MHD formulation and its implementation in PLUTO refer to MT10 MTB10 CHAPTER 6 BAS
58. Aga U x U ea se 12 1 2 Da lAa 130 Aa U Usos where U U is a conserved variable A 1 U are the undivided forward and backward differences in the direction d e g A 1U Uj 1 U see also section 4 1 in IMZT 12 The last term appearing in the denominator U4 ref prevents regions of small ripples from being refined and it is defined by Us ref Uia 2 0 Usa 12 2 with e 0 01 similar expressions hold when d x2 or d 23 By default U is the total energy density if the EOS is not isothermal or the mass density other wise see CHOMBO_REF VAR in Appendix B3 A different variable q q U e g q m2 2p can be used to replace U in Eq by copying the source file Src Chombo TagCells cpp in your local work ing area setting CHOMBO_REF_VAR to 1 and defining the appropriate expression through the function computeRefVar Beware however that x U may become ill defined if Uz ref changes sign This occurs for example when U is a vector component e g momentum or magnetic field and a better solution would be to replace Uy ref with a constant reference value 12 4 Output and Visualization PLUTO Chombo output employs the HDF5 format and the frequency of output is controlled at runtime by specifying the relevant parameters described in HDF5 is a data model library and file format for storing and managing large amounts of data It supports an unlimited variety of datatypes and is designed for flex
59. Block ooo oo ee ee ee 35 4 CLL TPR e eaa a aaa a ea ea aT 36 roo osas aaa 37 45 ThelBounaaryl DIOCK gt neu kee ee Gate eee G44 bao GE Ge eRe ee Se es 38 4 6 e state Gra OnipuUl DIOCK exter cda seras E EEE OEE HA SBE 39 Pe Be See Peewee oe eee eee ee eee a eee 40 4 8 eL amero DIOS 5 6 86 oa oo ds OO Oe we oe ee Se ey eA See 41 42 42 A 44 Serer ee 45 e Ese 47 Sueca enue see se44 eee ne eae ue oo ee 49 5 2 1 Internal Boundary cord ooo ea Shoo dG S S564 8S oR OSS 52 oa aa peas oa sarro 53 OM IE AIR 55 6 Basic Physics Modules 57 6 1 TheHDModule 0 0 0 0 ee a 57 Ook BOAON 2 256248 e 506255 5 8e8 BEES HE oS REE eee 57 62 The MHD Module 0 a 58 oael MOUGHOUG tees es Gee eee so See Peewee AAA 58 eee ee eee eee ee oes ta oe 59 oh aes ere eee eee cae oe ows 60 6 2 4 Background Field Splitting o o e 0 000 000 63 euros ia eres arar e 65 Gus OQUAUOUS 9er AAA 65 poderoso os Ge eue ee oe eae ees oe os 66 Ol TGUGUONS cig2 228 3 Bae eee eee Bee sa oe ee ee SS r 66 67 weer Tee ee eee ee ee eee eee 67 7 2 Me DEAL Equation Of Stale 4s 46 6 04 bho wo ew SES OS eGo EE SE a G eH wm 67 7 3 The PVTE_LAW Equation of State 2 ee 0 68 7 3 1 Example EOS for a Partially Ionized Hydrogen Gas in LTE 69 7 3 2 Analytic vs tabulated approach o o oo 70 74
60. D linear reconstruction is applied to primitive variables It is 2 order accurate in space Stencil is 3 point wide e WENO3 provides 3 order weighted essentially non oscillatory reconstruction YC09 inside a cell using is 3 point stencil e LimO3 provides 3 order limiter function CT09 based on a 3 point stencil e PARABOLIC piecewise parabolic method PPM as implemented by Mig14 The stencil requires zones The default is LINEAR Both WENO3 and LimO3 employ a local three point stencil to achieve piecewise quadratic reconstruction for smooth data and preserves their accuracy at local extrema thus avoiding clipping of classical second order TVD limiters and PPM Non uniform grid spacing and curvilinear coordinates are handled more correctly with LINEAR and PARABOLIC using the approach presented in Mig14 Note that although 3 order reconstructions are available the finite volume version of the code retains a global 2 4 order accuracy as fluxes are computed at the interface midpoint On the contrary genuine 3 and 5 order accurate schemes can be employed using the conservative finite difference framework 10 3 2 1 7 TIME STEPPING PLUTO has several time marching algorithms which can be used in either a spatially split or unsplit fashion If At t 1 t is the time increment between two consecutive steps and denotes the discretized spatial operator on the right hand side of Eq 1 1 the possible o
61. DERS myheader h This will instruct make that PLUTO has to be compiled and linked together with the user supplied file myfile c which depends on myheader h This is particularly useful when the user wants to compile and link the code together with supplementary routines contained in external files 4 Runtime initialization file pluto ini At start up the code checks for the pluto ini input file or a different one if the i command flag is given that contains all the run time information necessary for integration A template for this file can be found in the Src Templates directory The initialization file is divided into eight different blocks enclosed by a pair of square brackets Each block contains a set of labels and associated mandatory or optional fields Grid Xl grid 1 0 0 100 u Tr O x2 grid 1 0 0 100 u 10 x3 grid il 0 0 aL u Le Chombo Refinement Levels 4 Ref ratio 2 2222 Regrid_interval 2 2 2 2 Refine_thresh 0 3 Tag_buffer_size 3 Block tTactor 4 Max_grid_size 32 Fill_ratio dis TO Time CFL 0 4 CFL max_var etl CFL_par MS optional rmax_par 40 0 optional tstop 1 0 first_dt 1 e 4 Solver Solver tvdlf Boundary Xl beg outflow Xl end outflow X2 beg outflow X2 end outflow X3 beg outflow X3 end outflow Static Grid Output uservar 0 output_dir optional dbl RaO T single_file FLE 105 si single_file vtk Ls sit single_file optional dbl h5 1 0 2 40h
62. Doxygen documentation system and an on line documentation browser can be found in Doc test_problems html Test problem documentation is still being added and more examples will be available in future releases CHAPTER 0 QUICK START 10 0 6 Migrating from PLUTO 4 1 to PLUTO 4 2 PLUTO 4 2 provides several bug fixes and does not introduce major changes with respect to version 4 1 in the syntax of the basic functions defined in init c Some optimizations and improvements have been performed in the source distribution and a few minor changes have been introduced mainly for uniformity and efficiency reasons The file CHANGES lists the most relevant ones Structural Changes User defined constants have been removed from the python script and can now be added manually by editing a dedicated section in definitions h The old array indexing style using a two letter index e g DN PR VX is no longer supported Use RHO PRS VX1 instead The macro STS_nu has been changed to STS_NU it has been removed from the python script and it can be inserted manually in the user defined constant section of definitions h ENTROPY_SWITCH does not take the value YES NO anymore but it has more options see the documentation The Shearing Box module has been changed o It works with FARGO IDEAL Eos o The functions BodyForceVector and BodyForcePotential are slightly different from before check Test_Problems MHD Shearing_Box in
63. E 7 oE Ot Oz Ox OB Ey 7 OE Ot Ox Oy 0 _ ae oe _ ae o 9 57 Oz A 1 pu g 0 0 0 where v Ur Yy vz and B Bz By Bz are the velocity and magnetic field vectors Ez Ey Ez are the components of the electromotive force E v x B gis the body force vector and Y is the gravitational potential 114 A 1 2 Polar Coordinates In polar cylindrical coordinates R z the conservative ideal MHD Equations are discretized using the following divergence form Op Ml aT V pv eS OMR Op o puz Ba J V mpv BrB JR Pp ox 5 R OMe R 1 Op 1 0 BB Om Op OD zU B B a ae O Z E 08 V E p po v B v B pvu g OBR n 10 OEg 9 Ot R 09 Oz OB Er OE _ 9 Ot Oz OR o OB 1O REs 1 OER o ot R R R Oo Note that curvature terms are present in the radial component while the azimuthal component is dis cretized in angular momentum conserving form The corresponding divergence operators are 10 RFr 10F OF R OR ROS Oz 1 O RFFR 10F OF R2 OR R06 0z V F A 3 vep In the previous equations v vR vg vz and B Br Bg B are the velocity and magnetic field vectors Er Ep Ez are the components of the electromotive force E v x B gis the body force vector and is the gravitational potential 115 A 1 3 Spherical Coordinates In spherical coordinates r 0 p the ideal
64. I 98 11 3 1 Visualization with Gnuplot oo o 98 11 3 2 Visualization WihIDL coords rosada eho ED Eo 100 oe See Ge OE ae Gs oe a Ae we a eee ee 102 poe etaaneee tans eoeke E a 104 jo Sede ere oe nee oe ese eee 105 12 Adaptive Mesh Refinement AMR 106 Aes eee tee Ghee ee eRe Hee eee eee eee eee ae ey ee 106 12 LL nstaline FIDE oo oe bok eon AAA AAA ee 107 See nee Ge De eee ae eae es eee ee ee Se 107 Dee eee eee eee eee eee eee eee 109 ro ore a aaa er sa 109 125 Controlling Kennementt sese arca amp oe pe AAA Boe oe 110 E E ae ee ee 110 12 4 1 Visualization withIDO 00 0 0 0000000000008 110 12 4 2 Visualization with VisIt 0 0 a eee eee eee 112 12 4 3 Visualization with pyPLUTO o ooo oo o 113 114 a a ES 114 A ee E er ee 114 ar sao osa a 115 A 1 3 Spherical Coordinates oo essa esas a ee aw as 116 A 2 Special Relativistic MHD Equations o o 117 A 2 1 Cartesian Coordinates ee 117 A 2 2 Polar Coordinates 6 0 o a a 117 B A 2 3 Spherical Coordinates Predefined Constants and Macros B 1 Predefined Physical Constants B 3 Advanced Options B 2 Predefined Function Like Macros 0 Quick Start 0 1 Downloading and unpacking PLUTO PLUTO can be downloaded from http plutocode ph unito it Once downloaded extract all the files fro
65. IC PHYSICS MODULES 61 6 2 3 3 Constrained Transport CT In this formulation BS99 Ld04 GS05 two sets of magnetic fields are used e face centered magnetic field b hereafter e cell centered magnetic field B hereafter The primary set is the first one where the three components of the field are located at different spatial points in the control volume that is IR Oro ij 1 k y Oro ij del see Fig In Cartesian coordinates for instance b is located at X faces whereas by lives at Y faces etc see the boxes and triangles in Fig This feature must be used only in conjunction with an unsplit integrator With CT the solenoidal condition is maintained at machine accuracy as long as field initialization is done using the vector potential 46 2 2 The staggered magnetic field is treated as an area weighted average on the zone face and Stoke s theorem is used to update it Ob dbz 1 e Grvxe 4se 0 E fed 0 6 9 Please note that the staggered components are initialized and integrated also on the boundary interfaces in the corresponding staggered direction In other words the interior values are IBEG 1 lt 1i lt IEND b 1 2 JBEG lt j lt JEND 21 45 j k i KBEG lt k lt KEND IBEG lt i lt IEND cia a eo D KBEG lt k lt KEND IBEG lt 1 lt IEND Desig kth JBEG lt y lt JEND KBEG 1 lt k lt KEND Thus boii 45 is NOT a boundary value for i IBEG 1 JBEG lt 7 lt JEND KBEG lt k
66. LUTO INI 38 4 5 The Boundary Block Specifies the boundary conditions to be applied in the ghost zones of the computational domain X1 beg string X1 end string X2 beg string X2 end string X3 beg string X3 end string Assuming that q is a scalar quantity and n is the coordinate direction orthogonal to the boundary plane st ring can be one of the following outflow B a y v_y 0B_ a oe CU Sets zero gradient across the boundary LE LeCrtive Reflective rigid walls boundary conditions Variables are symmetrized across the boundary and normal components of vector fields flip signs Un gt Un Vt gt Ut qq Bn gt Bn B gt B where n t is normal tangential to the interface axisymmetric Same as reflective except for the angular component of vg or By which also changes sign Un gt Un Ud U Vagis 7 Varis q m q J Da gt Dy Bg Bg ais r Daris where axis is r 0 z or r 0 0 in cylindrical or spherical coordinates egtsymmetric Sets equatorial symmetry with respect to a given plane It is similar to reflective but with reversed sign for the magnetic field Un gt Un Ui gt Ut qq Ban gt Bn B gt B periodic Sets periodic boundary conditions on both sides of the computational domain shearingbox Shearingbox boundary conditions are similar to periodic except that they are sheared in one di rection only X1 beg and X1 end suppor
67. Me No MeNe Y Mp Np No p 2N Np mMp kpT p Weng tro kel 14 a ng kl E p where u 1 1 zx is the mean molecular weight x np Nnp no is the degree of ionization computed from Saha equation 2 3 2 Xv _ QamekpT e x kaT 7 10 l z h np no l where x 13 6 eV and n no p Mp The specific internal energy includes two contributions 3kgT 3 e 12 ER 7 11 2 ump Mp 2p Mp where the first one represents the standard kinetic energy while the second one corresponds to the ionization energy neutral atoms have a potential energy that is lower than that of ions Note that the latter introduces a temperature or equivalently a velocity scale in the problem so that computations are no longer scale invariant but depend on the value of UNIT_VELOCITY used to obtain the temperature in Kelvin and UNIT_DENSITY used in Saha equation that must be defined in your Init function see Fig 7 1 shows the classical Sod shock tube solution at 0 2 obtained with the IDEAL equation of state and the PVTE_LAW with UNIT_DENSITY 10 m and UNIT_VELOCITY 10 Km s The equivalent I defined as p Peg 1 7 12 pe is no longer a constant but a function of the temperature see Fig The implementation of this particular EOS can be found in Src EOS PVTE pvte_law_H c simply copy it to your working directory as pvte_law c Eq 7 11 is implemented by the InternalEnergyFunc function while the mean molecular weigh
68. O_1 will contain the actual value of the first user defined parameter g inputParam FOO_2 the second one and so forth The maximum number is 31 and parameter names should be chosen with care in order to avoid overlapping conflicts with names that are already defined in the code Although there are no strict rules we strongly advise to use capital letters to avoid short labels such as VO or VX and to choose a more representative name that explains the use of the variable on its own e g PAR_INFLOW_VEL Parameter names and values are automatically inserted inside pluto ini in the correct order after the execution of the python script However if you use a different initialization file than pluto ini you may have to set the parameter names together with their values manually CHAPTER 2 PROBLEM HEADER FILE DEFINITIONS H 23 2 2 Physics Dependent Options After the physics module has been selected further options become available in a following menu Depending on the selected physics module different options may appear They are described in the following 2 2 1 BACKGROUND FIELD If set to YES the magnetic field is split into a static curl free contribution and a time dependent deviation This option is only available for the MHD module and it is described in detail in 46 2 4 2 2 2 DIVB CONTROL Numerical methods do not naturally preserve the condition V B 0 For this reason this option allows the user to select a contro
69. P_IN the high pressure of a region yet to be specified P OUT the ambient pressure and GAMMA polytropic index CHAPTER 0 QUICK START 9 P_IN 8 e2 POUT 8 0 GAMMA 1 6060666066666566 71 Save and exit the editor Next you need to edit init c e Define inside the function Init the radius r a floating point value which we will be used to set a circular region of high pressure double r e Set the global variable g_gamma polytropic index and the radius r Define the initial ambient pressure P_OUT and put an IF statement to specify the high pressure region inside a circle of r 0 3 P_IN g_gamma g_inputParam GAMMA calls the auxiliary parameter GAMMAx xlexl x2 x2 1 0 x initial density array x 0 0 f initial Vx array 0 0 je initial Vy array 0 0 gi je initial Vz array nputParam P_OUT calls the auxiliary parameter P_OUT x 0 3 v PRS g_inputParam P_IN calls the input parameter P_IN x Save and exit the editor Compile the code and run PLUTO with a the following set of commands Blastwave gt make Blastwave gt pluto In order to visualize the results follow the instructions described in the two previous sections 0 5 Supplied test problems The official distribution of PLUTO comes with several examples and test problems that can be found under the Test_Problems folder Documentation is extracted from comments at the beginning of init c sources files using the
70. Problems MHD Rayleigh_ Taylor or Test_Problems MHD Shearing Box setups where reflective and equatorial symmetric boundaries are used in the y and z directions Note Relativistic flows Body forces are marginally compatible with the relativistic modules Only VECTOR may be used 5 4 The Analysis function The Analysis function can be used to perform run time data analysis reduction in order to save intensive I O for data post processing This function call frequency is set in pluto ini see Syntax void Analysis const Data xd Grid grid Arguments e da pointer to the PLUTO data structure as in e xgrid a pointer to an array of Grid structures containing all the relevant grid information Example In the next example we show how to compute at run time the total integrated kinetic energy and the maximum internal energy ol l 4 1 1 9 B D Fin Ay J zPY dx dy dz AV 2 5 Ping kV AVi j k PC rias mak 75 i j k 1 where AV is the total volume of the computational domain AV jx Azr Ay Az is the volume of a single cell i j k extend over the entire computational domain a Cartesian 3D domain is used for simplicity The output file name is averages dat and it is written as a 4 column tabulated ascii file con taining the current integration time the time step the volume integrated kinetic energy and maximum internal energy for the required time level The example works also for parallel computat
71. RMAL or IDEAL equations of state Initial conditions are specified as usual in the Init function where the orbital speed must be set to vy gQox A simple example corresponding to p 1 p c2p and plasma 10 is given below Isothermal sound speed x 1 0 LL 0 0 SB_Ox SB_OMEGAxx1 Orbital velocity x 0 0 if EOS IDEAL v PRS v RHO csx cs elif EOS ISOTHERMAL g_isoSoundSpeed cs endif Hif PHYSICS MHD beta l es v BX1 0 0 v BX2 0 0 v BX3 csx sqrt 2 0 beta Vertical field net flux x endif The numerical value of q and 29 is prescribed starting with PLUTO 4 2 using the macros SB_Q and SB_OMEGA which by default take the value of 3 2 and 1 respectively Should you change the default value add them as user defined constants as explained in Only gravitational forces must be given through the BodyForceVector function since Coriolis term are separately included by PLUTO An example containing several configurations can be found in the Test_Problems MHD Shearing Box directory Boundary conditions must be prescribed as shearingbox at the X1_BEG and X1_END boundaries periodic in the azimuthal y direction but can be freely assigned in the vertical direction z Compatibility with the FARGO module The shearingbox module is fully compatible with the FARGO algorithm and a significant gain may be obtained for boxes with large aspect ratio L gt L To enable both modules
72. UNIVERSITA DEGLI STUDI 24 A DI TORINO fish fees ALMA UNIVERSITAS E TAURINENSIS y SCAT SuperComputing Applications and Innovation PLUTO v 4 2 Aug 2015 User s Guide http plutocode ph unito it e mignone ph unito it Developer A Mignone Contributors C Zanni AMR zanni oato inaf it B Vaidya EoS Cooling pyPLUTO bhargav vaidya unito it T Matsakos Resistivity Thermal Conduction STS G Muscianisi Parallelization I O P Tzeferacos Viscosity MHD STS Finite Difference O Tesileanu Cooling 1 Dipartimento di Fisica Turin University Via P Giuria 1 10125 Torino TO Italy 2 INAF Osservatorio Astronomico di Torino Via Osservatorio 20 10025 Pino Torinese TO Italy 3 Dept of Astronomy amp Astrophysics University of Chicago 5640 S Ellis Ave Chicago IL 60637 USA 4 Consorzio Interuniversitario CINECA via Magnanelli 6 3 40033 Casalecchio di Reno Bologna Italy 5 FLASH Center University of Chicago USA 6 Department of Physics University of Bucharest Str Atomistilor nr 405 RO 077125 Magurele Ilfov Romania Terms amp Conditions of Use PLUTO is distributed freely under the GNU general public license Code s development and support requires a great deal of work and for this reason we expect PLUTO to be referenced and acknowledged by authors who use it for their publications Co authorship may be solicited for those publications demanding considerable
73. _DENSITY 1 67e 23 define UNIT_LENGTH 3 1e18 define UNIT_VELOCITY 1 65 x End user defined constants do not change this line x Note that when the relativistic modules are used vg must be the speed of light Output files can be directly saved in cgs units using the flt or vtk data format see 94 6 Example 2 Consider a flow with typical number densities of the order of n 10 cm temperatures T 10 K corresponding to typical sound speeds of cso 10 Km s Suppose also that the flow propagates with uniform speed v 50 Km s and the typical scale size of the problem is L 1 pc 3 1 10 cm Then one may choose Po Nom 1 67 107 gr cm vo 1Km s 10 cm s Lo 3 1 101 cm From the python script this is done by including the following line in definitions h x Beg user defined constants do not change this line x define UNIT_DENSITY 10 0 CONST_mp define UNIT _LENGTH CONST_pc define UNIT VELOCITY 1 e5 x End user defined constants do not change this line where CONST_mp and CONST_pc are pre defined symbolic constants proton mass and parsec in c g s units and are defined together with several other constants in Appendix B 1 Please remem ber using parenthesis around a macro expression to avoid incorrect expansion With this choice of units the piece of code describing the initial condition becomes CHAPTER 5 INITIAL AND BOUNDARY CONDITIONS INIT C 45 LOs means 1 1 67
74. _MHD BOX LOOP box k 7 2 d gt Vs BX2s k jJ i lt 30 inputParam B PRE endif Jelse if box gt vpos X3FACE x z staggered field x ifdef STAGGERED_MHD BOX_LOOP box k j 1 d gt Vs BX3s k Jj 1 g_inputParam B_PRE endif As in the previous example the macro BOX_LOOP is interchangable for cell centered data box gt vpos CENTER with the macro X1_END_LOOP k j i but not rigorously for staggered magnetic fields which are defined on a larger stencil Function like macros are described in the code documentation Doc Doxygen html macros_8h html 5 2 1 Internal Boundary When UserDefBoundary is called with side 0 and INTERNAL_BOUNDARY has been turned to YES inside your definitions h the user has full control over the solution array This feature can be used to adjust the value of selected cell centered primitive variables inside a specific region of the computa tional domain rather than at boundaries In this case the TOT_LOOP macro should be employed to loop over the local computational domain and a user defined criterion typically spatially or variable dependent is used to modify the solution array in the selected zones A typical example may occur when a lower or upper threshold value should be imposed on phys ical variables such as density pressure or temperature For instance the following piece of code sets a floor value of 107 on density void UserDefBoundary const Data d
75. a loop over the bottom boundary zones and for cell centered data it is equivalent to the macro X2_BEG_LOOP k j i Similar macros may be used at any of the other boundaries X1_BEG X1_END X2_END X3_BEG X3_END although the BOX_LOOP macro has the advantage of being more general since it automatically embeds the stencil index range for the corresponding variable position i e centered or staggered CHAPTER 5 INITIAL AND BOUNDARY CONDITIONS INIT C 32 Example 2 As a second example we discuss the user defined boundary condition employed in the shock cloud problem Test_Problems MHD Shock_Cloud Here we want to prescribe at the X1 end boundary con stant pre shock values on both cell centered quantities and staggered magnetic fields The variable box gt vpos is used to select the desired data set void UserDefBoundary const Data d RBox box int side Grid xgrid ine i jki select the boundary side x select the variable position x Loop over boundary zones x if side X1_END if box gt vpos CENTER BOX_LOOP box k J3 1 d gt Vc RHO k J EXPAND d gt Vc d gt Vc d gt Vc d gt Vc PRS k EXPAND d gt Vc d gt V6 d gt Vc ES e A L 12536 j OAs P 0 0 0 0 g_ puras PRE g_inputParam B_PRE H H H O H H H O bo e e Ne LH e LI Ne Jelse if box gt vpos X2FACE x y staggered field x ifdef STAGGERED
76. a two shock Riemann solver Roe or hlld alternatively Setting DIMENSIONAL_SPLITTING NO yields the spatially unsplit fully corner coupled method of Col90 MPBO05 This scheme is stable under the condition CFL lt 1 in 2D and CFL lt 1 2 in 3D and it is slightly more expensive than RK2 Time Step Determination The time step At is computed using the information available from the previous integration step and it can be controlled by the Courant Friedrichs Lewy CFL number C a within the limits suggested in Table 2 1 see Bec92 Thus one immediately sees that if Al is the cell physical length the time step roughly scales as Al for hyperbolic problems and as Al when 18 4 1 On the contrary when parabolic terms are included via Super 18 4 2 the time step can be much larger being computed solely from the advection time scale i e Tg 0 in the table below Multi step algorithms RK2 RK3 work in all system of coordinates and are the default choice Single step schemes HANCOCK CHARACTERISTIC_TRACING are more sophisticated have less dissipation and have been tested mainly on Cartesian and cylindrical grids Have a look at Table 2 2 for a compari son between different suggested integration schemes commonly adopted in testing the code SCHEME DT Ms GPLEIT CFL Limit Ad 274 RK YES At l Ca lt 1 Ea max i xa gt Ad aTa ANCK CATE YES At max max J Ca lt 1 d ijk NAla Als
77. active processing of large amounts of data The PLUTO code distribution comes with a number of useful routines written in the IDL programming language to read and visual ized data written with PLUTO Several functions and procedures for data visualization and analysis can be found in the Tools IDL directory which we strongly advise adding to the IDL search path IDL function documentation can be found in Doc idl_tools html The PLOAD Procedure The PLOAD procedure can be considered the main driver for reading data stored in one of the following formats dbl flt vtk dbl h5 flt h5 or hdf5 PLOAD requires the infor mation stored in the corresponding data log file e g grid out dbl out flt out etc to initialize common block variables and grid information shared by other functions and procedures in the Tools IDL subdi rectory Because of this reason it should always be called at the very beginning of an IDL session For example IDL gt Read all variables from the the third output in double precision IDL gt gt and store it into memory rho vxl vx2 prs zes IDL gt PLOAD 3 IDL gt Read only pressure from the fifth output in single precision IDL gt and store it into memory IDL gt PLOAD 5 FLOAT VAR prs IDL gt Read output 9 from the directory Large_Data do not store it IDL gt into memory but use IDL file association preferred for large datasets IDL gt PLOAD dir Large_Data 9 ASSOC By default
78. additional support and or changes to the code Contents 0 Quick Start 0 3 Running the Orszag Tang MHD vortex test 0 6 Migrating from PLUTO 4 1 to PLUTO 4 2 1 Introduction paa oe Gee obese Benson Esa PA o aos wad asa a se ao poesias a 2 Problem Header File definitions h 2 1 1 PHYSICS 2 1 2 DIMENSIONS amp COMPONENTS GEOMETRY BODY FORCE COOLING RECONSTRUCTION 2 1 7 TIME STEPPING 2 1 8 DIMENSIONAL SPLITTING NTRACER 2 1 10 USER DEF_ PARAMETERS 2 2 Physics Dependent Options 2 2 1 BACKGROUND_FIELD Choe es p bY on bots be bee ae eee ong adae en cane DOS HON a shades Bede h nae weae Shawn eRA ooh ode chee obekeas 2 24 ENTROPY SWITCH 229 RESTOTIVITY 2 2 6 ROTATING_FRAME 2 2 7 THERMAL CONDUCTION 220 a L ad rA ee na AAA 2 3 User defined Constants 2 4 Frequently Used Options ooa sas ESA RAS 3 Makefile Selection makefile 29 3 1 MPI Library Parallel Support 2 824 445s 6 aaa aa a 29 Ol Asynehrounousl Ol 26 25 a5 6S ceda A 30 32 IDES o es o gt tobe eke heer Stee erreen ees eee Ss 30 J9 PNG Library SUPPO 4 resis BGS OS bd Se SS EERE TELE EP ETS 31 3 4 Including Additional Files locall_Make oo o 31 4 Runtime initialization file pluto ini 32 4 1 The Grid block o 33 4 2 The Chombo Refinement
79. ader prompts for the output number single or double precision and then the variable to display The out put of this session is shown in Fig for the MHD Shock Cloud interaction test conf 01 in Test_Problems MHD Shock_Cloud The function ArrayPlot is used to display 2D datasets and is directly included in the interface reader For instance to change colormap and visualize pressure in log scale use ArrayPlot Log data 8 DataReversed gt True ColorFunction gt RoseColors DataRange gt xmin xmax ymin ymax For different colormaps consult http reference wolfram com language guide ColorSchemes html A 1D cut of pressure for instance through the x direction at constant y be plotted using ListPlot data 8 ny 2 The current directory can be displayed by typing Directory while the home directory can shown by typing SHomeDirectory 8http www wolfram com CHAPTER 11 OUTPUT AND VISUALIZATION 105 11 3 5 Visualization with VisIt or ParaView PLUTO data written using VTK or HDF5 both h5 and hdf5 files formats can be easily visualized using either VisIt or Para View available at https wci lInl gov codes visit home html and http www paraview org respectively VisIt is an open source interactive parallel visualization and graphical analysis tool for viewing scientific data ParaView is an open source mutiple platform application for interactive scien tific visualization An example is shown in Fig
80. aim is the number of spatial dimensions while for dimensionally split methods one has CFL S 1 Certain combinations of algorithms may have more stringent limitations a second order Runge Kutta algorithm with parabolic reconstruction for example requires CFL lt 0 4 for stability reasons CFL_max_var double Maximum value allowed for At At maximum time step growth between two consecutive steps CFLpar double optional When parabolic terms are integrated using operator splitting with Super Time Stepping 48 4 2 it controls the diffusion Courant number The default value is 0 8 Naim This parameter has no effect when parabolic terms are included via standard explicit integration rmax_par double optional When parabolic terms are integrated using operator splitting with STS it sets the maximum ratio between the actual time step and the explicit parabolic time step i e At At The default value is 100 This parameter has no effect when parabolic terms are included via standard explicit integration tstop double Integration ends when t t stop unless a maximum number of steps 1 4 is given tstop has to be gt 0 0 first det double The initial time step A typical value is 107 CHAPTER 4 RUNTIME INTITALIZATION FILE PLUTO INI 37 two_shock roe ausm hlld hllc hll tvdlf HD E vo Vv vv MHD a ANY RHD Vo ee y RMHD ee y Table 4 1 Available Riemann solvers for the different phys
81. al ysis 5 1968 506 517 John Lighton Synge The relativistic gas North Holland Publishing Company Interscience Publishers Amsterdam New York 1957 O Tesileanu A Mignone and S Massaglia Simulating radiative astrophysical flows with the PLUTO code a non equilibrium multi species cooling function A amp A488 2008 429 440 Eleuterio F Toro Riemann solvers and numerical methods for fluid dynamics 2 ed Springer Verlag Berlin Heidelberg 1997 B van Leer Towards the ultimate conservative difference scheme V A second order sequel to Godunov s method Journal of Computational Physics 32 1979 101 136 B Vaidya A Mignone G Bodo and S Massaglia Astrophysical fluid simulations of thermally ideal gases with non constant adiabatic index numerical implementation A amp A580 2015 A110 J Woodall M Ag ndez A J Markwick Kemper and T J Millar The UMIST database for astrochemistry 2006 A amp A466 2007 1197 1204 N K Yamaleev and M H Carpenter A systematic methodology for constructing high order energy stable WENO schemes Journal of Computational Physics 228 2009 4248 4272 126
82. an be added to your init c to assign Bo void BackgroundField double xl double x2 double x3 double BO Note that when writing output datafiles to disk only the deviation B is written Examples can be found in the 4 configuration of Test _Problems MHD Rotor and in the 4 or 5 configurations of Test_ Problems MHD Blast Note Background field splitting works at present with the CT and GLM divergence cleaning techniques with most Riemann solvers but only with RK type integrators CHAPTER 6 BASIC PHYSICS MODULES 65 6 3 The RHD Module The RHD module implements the equations of special relativistic fluid dynamics in 1 2 or 3 dimensions Velocities are always assumed to be expressed in units of the speed of light The special relativistic module comes with 2 different equations of state and it also works in curvilinear coordinates Gravity in Newtonian approximation can also be incorporated The relevant source files and definitions for this module can be found in the Src RHD directory 6 3 1 Equations The special relativistic module evolves the conservative set U of state variables T U D m1 Ma MZ E where D is the laboratory density M 1 2 13 are the momentum components E is the total energy in cluding contribution from the rest mass The evolutionary conservative equations are D Dv m V mv ol 0 Ot dd E E m where v is the velocity p is the thermal pressure Primitive variables V always include
83. an in code units vtk double integer string string cgs Like 1t but for VTK legacy file format see 11 1 3 dolne double integer string Like db1 but for haf 5 double precision format 411 1 2 This format can also be used for restarting the code by supplying the h5restart n command line argument cdo ais double integer string Like db1 but for ha 5 single precision format 411 1 2 tab double integer string Sets the time and the number of steps interval for tabulated ascii format 11 1 4 ppm double integer string Sets time and the number of step intervals for 2D color images in PPM format 411 1 5 CHAPTER 4 RUNTIME INTITALIZATION FILE PLUTO INI 40 e png double integer string Sets time and the number of step intervals for 2D color images in PNG format 411 1 5 e log integer Sets the interval in number of steps of the integration log to screen or pluto log e analysis double integer Sets time and number of steps interval between consecutive calls to the function Analysis see 4 7 The Chombo HDF5 output Block Relevant only for AMR Pluto with the Chombo library 12 this block controls how often restart and plot files are dumped to disk in the AMR version of the code The relevant options are e Output dir string Sets the output directory name for writing and reading simulation data Both plotfiles and check point files will be written to and read from this directory If Output _d
84. an unsplit fashion F Fip t Fia 8 10 where hyp and par are respectively the hyperbolic and parabolic fluxes see also 3 1 of MZT 12 Such methods are however subject to a rather restrictive stability condition since in the diffusion dominated limit At Al n where n is the maximum diffusion coefficient see Table 2 1 for the exact limiting factor Clearly high resolution and large diffusion coefficients may lead to drastic reduction of the time step thus making the computation almost impracticable 8 4 2 Super Time Stepping STS STS AAG96 is a technique that considerably accelerates the standard explicit treatment of parabolic terms In this case parabolic terms are treated in a separate step using operator splitting and the solution vector is evolved over a super time step equal to the advective one The superstep consists of N sub steps properly chosen for optimization and stability depending on the diffusion coefficient the grid size and the free parameter y lt 1 STS_NU 2N _ 1 _ 5 2N At Pe eee Osos era with Atpar T Da ea a eo ee 2 Da Cp 8 11 Here Atpar is the explicit parabolic time step computed in terms of the diffusion coefficient D and phys ical grid size Al The previous equation is solved to find N for given values of At Atpar and v For v 0 SIS is asymptotically N times faster that the standard explicit scheme However very low val ues of y may result in a
85. ar problems involving strong discontinuities to limit the scale disparity between advection and diffusion time scales by restricting the time step At to be at most Tmax Atpar With Atpar defined by Eq 8 11 and rmax a free parameter see 4 3 In this cases max may be lowered by lowering rmax par in pluto ini from its default value 100 to 40 or even less Note that although this method is in many cases considerably more efficient than the explicit one it is found to be slightly less accurate due to operator splitting The method is by definition first order accurate in time although different values of the y parameter are found to affect the accuracy On the other hand STS bypasses the severe time constraint posed by second derivative operators in high resolution simulations During the STS step momentum magnetic field or total energy are evolved in time even if the ENTROP Y_SWITCH has been enabled 9 Optically Thin Cooling PLUTO can include time dependent optically thin radiative losses in a fractional step formalism in which the hydrodynamical evolution and the source step are solved separately using operator splitting This preserves 2 order accuracy in time if both the advection and source steps are at least 2 order accurate During the cooling source step specifically PLUTO solves the internal energy and chemical reaction network equations 7 A n T X XX gy ny 9 1 Ot a i where A is a cooling or heating te
86. are allowed The grid type u s or 1 is ignored and a uniform grid is always assumed unless the CHOMBO_LOGR switch is enabled to generate a logarithmic radial grid see Appendix B 3 Cells must not necessarily have the same physical length in each direction e g squares in 2D cubes in 3D and can have an aspect ratio different from 1 The refinement options are set in the Chombo Refinement section CHAPTER 4 RUNTIME INTITALIZATION FILE PLUTO INI 35 4 2 NX1_TOT Figure 4 2 Computational grid in 2 dimensions with NX1 NX2 4 and 1 ghost zone Internal zones solid boxed are spanned by IBEG lt i lt IEND JBEG lt j lt JEND Dashed boxes represent boundary ghost zones The Chombo Refinement Block Controls the grid refinement if PLUTO has been compiled with the Chombo library It is ignored otherwise The relevant parameters for refinement are Levels integer Sets the finest allowable refinement level starting from the base grid level 0 defined by the Grid block 0 means there will be no refinement Ref ratio integer integer Sets the refinement ratios between all levels First number is ratio between levels 0 and 1 second is between levels 1 and 2 etc There must be at least Levels 1 elements or an error will result Regrid_interval integer integer Sets the number of timesteps to compute between regridding A negative value means there will be no regridding There must be at least Le
87. ation SIAM Journal on Numerical Analysis 29 1992 no 3 701 713 D S Balsara and D S Spicer A Staggered Mesh Algorithm Using High Order Godunov Fluxes to Ensure Solenoidal Magnetic Fields in Magnetohydrodynamic Simulations Journal of Computational Physics 149 1999 270 292 T J M Boyd and J J Sanderson The Physics of Plasmas January 2003 D S Balsara D A Tilley and J C Howk Simulating anisotropic thermal conduction in supernova remnants I Numerical methods MNRAS386 2008 627 641 L L Cowie and C F McKee The evaporation of spherical clouds in a hot gas I Classical and saturated mass loss rates ApJ211 1977 135 146 Phillip Colella A direct eulerian muscl scheme for gas dynamics SIAM Journal on Scientific and Statistical Computing 6 1985 no 1 104 117 P Colella Multidimensional upwind methods for hyperbolic conservation laws Journal of Computational Physics 87 1990 171 200 M Cada and M Torrilhon Compact third order limiter functions for finite volume methods Journal of Computational Physics 228 2009 4118 4145 P Colella and P R Woodward The Piecewise Parabolic Method PPM for Gas Dynamical Simulations Journal of Computational Physics 54 1984 174 201 L Del Zanna N Bucciantini and P Londrillo An efficient shock capturing central type scheme for multi dimensional relativistic flows II Magnetohydrodynamics A amp A400 2003 397 413 A Dedner F Kemm D K
88. ault The default is to write ALL fields in dbl format whereas to exclude staggered magnetic field components if any from the flt format 11 1 2 HDF5 Output dbl h5 or flt h5 data formats HDF5 output format can be used in the static grid version if PLUTO has been succesfully compiled with the serial or parallel version of the HDF5 library see The file extension is h5 and not hdf5 as used by PLUTO Chombo data files and output files are written in Pixie format a single block rectilinear mesh file using HDF5 may be related to the Polar lonospheric X Ray Imaging Experiment The conventions used in writing dbl h5 or flt h5 files are the same ones adopted for the dbl and flt data formats However with HDF5 all variables are written to a single file Pixie files can be opened and visualized directly by different softwares like VisIt and Paraview Since we have found compatibility issues with some versions of these visualization softwares each file comes along with a supplementary xmf text file in XDMF format that describes the content of the correspond ing HDF5 file and can be opened correctly by VisIt and Paraview see Restart can be performed from double precision HDE5 data files by invoking PLUTO with the h5restart n command line option 1 4 1 where n is the output file number from which to restart In this case an additional file restart out will be dumped to disk Default The default is to write ALL fields in dbl h5 format
89. beg gt lt xl_end gt lt nx1 gt point s lt ngh gt ghosts X2 lt x2_beg gt lt x2_end gt lt nx2 gt point s lt ngh gt ghosts X3 lt x3_beg gt lt x3_end gt lt nx3 gt point s lt ngh gt ghosts KKEKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKK The rest of the file is made up of 3 sections one for each dimension giving the interior number of point followed by a tabulated multi column list containing from left to right the point number left and right cell interfaces nxl lt point number gt lt cell left edge gt lt cell right edge gt and similarly for the x2 and x3 directions CHAPTER 11 OUTPUT AND VISUALIZATION 96 11 2 Customizing your output Output can be customized by editing two functions in the source file Src userdef_output c in the PLUTO distribution We recommend to copy this file into your working directory and modify the default set tings if necessary Changes can be made by i introducing new additional variables and ii altering the default attributes 11 2 0 1 Writing Supplementary Variables New variables can be written to disk in any of the available formats previously described The number and names of these extra variables is set in your pluto ini initialization file under the label uservar The function ComputeUserVar located inside Src userdef_output c tells PLUTO how these variables are computed As an example suppose we want to compute and write te
90. bling this selective update neither the total energy nor the entropy will generally be conserved at the numerical level eo ALWAYS the entropy equation is used everywhere in the computational domain to update the solution array i e F a 1 always in Eq 2 5 This choice is consistent only with smooth flows e CHOMBO_REGRID the energy equation is used everywhere in the computational domain that is F o 0 in Eq 2 5 However pressure is still computed from entropy after projection coarse to fine prolongation and restriction operations AMR only This option violates energy conservation but has the advantage of preserving entropy and pressure positivity in those situations where kinetic and or magnetic energies are the dominant contributions to the total energy density The ENTROPY_SWITCH is also compatible with super time stepping although it will only be used during the hydro steps Also beware that in the current code release the ENTROPY_SWITCH may not compatible with all divergence control methods in resistive MHD CHAPTER 2 PROBLEM HEADER FILE DEFINITIONS H 25 2 2 5 RESISTIVITY Include resistive terms in the MHD equations see The available options are e NO resistivity is not included e EXPLICIT resistivity is included explicitly e SUPER_TIME_STEPPING resistivity is treated using super time stepping 48 4 2 2 2 6 ROTATING FRAME When set to YES it solves the equations in a frame of reference rotating with con
91. boxes can be overplotted using the oplotbox routine Here we plot the boxes for all 4 refine levels in addition to the base coarse grid In 7 I oplotbox D AMRLevel lrange 0 4 cval r g k c m geom D geometry The figure shows the total magnetic pressure obtained for the MHD Rotor problem in 2D at the equivalent resolution provided by the 4 refinement level also overplotted are the AMR levels in different colored lines for all of these 4 levels MagneticPressure Figure 12 3 AMR Data visualisation using pyPLUTO Note The HDF5 reader is not yet integrated into the pyPLUTO s graphical user interface 113 A Equations in Different Geometries In this section we give the explicit form of the MHD and RMHD equations written in different systems of coordinates Non ideal terms such as viscosity resistivity and thermal conduction are not included here The discretizations used in the Src MHD rhs c and Src RMHD rhs c strictly follow these form Equations for the non magnetized version HD and RHD are obtained by setting the magnetic field vector B 0 A 1 MHD Equations A 1 1 Cartesian Coordinates In Cartesian coordinates x y z the conservative ideal MHD Equations are discretized using the following divergence form Op aT V pv OM OP OMy Opt BB PE Om Opt O gE pb V gt E pe p v B v B OB J oE 7 Ey Ot Oy Oz OBy 4 O
92. braries are automatically named by Chombo after the chosen configuration The default configuration can be set by editing manually Chombo 3 2 lib mk Make defs local where depending on your local system and configuration you need to set make variables To this end gt 2d CROMDO F27 14 6 gt make setup create Make defs local from template a imi The command make setup will create this file from a template that contains instructions for set ting make variables that Chombo uses These variables specify the default configuration to build what compiler to use together with its flags where the HDF library can be found and so on At this point you should edit Make defs local The normal procedure is to define a default configura tion e g 2D serial Configuration variables DIM 2 DEBUG FALSE OPT TRUE PRECISION DOUBLE PROFILE FALSE CXX g EC gfortran MPI FALSE Note don t set the MPICXX variable if you don t have MPI installed MPICXX mpic OBJMODEL XTRACONFIG Optional features USE_64 USE_COMP LEX USE_EB USE_CCSE SE_HDF TRUE DF INCFLAGS I usr local lib HDF5 serial include DFLIBFLAGS L ust local lib ADF5 serial lib 1hdf5 1z Note don t set the HDFMPIx variables if you don t have parallel HDF installed DFMPIINCFLAGS I usr local lib HDF5 parallel include DEMPILIBFLAGS L usr local lib HDF5 parallel lib lhdf5 1z oot oo GH AE Se Defaults are used for the remaini
93. c gnuplot html Ascii Data Files If you enabled the tab output format in pluto ini you can plot 1D data from e g the first output file by typing gnuplot gt plot data 0001 tab 3 title Density gnuplot gt replot data 0001 tab 9 u Ls u Ls title Pressure overplot Here the first column corresponds to the x coordinate the second column to the y coordinate and flow data values start from the third column Fig shows the density and pressure profiles for the Sod shock tube problem conf 03 in Test_Problems HD Sod using the previous commands Two dimensional ascii datafiles can also be visualized using the splot command Fig shows a simple contour drawing of the final solution of the Mach reflection test problem remember to enable tab output using gnuplot gt set contour gnuplot gt set cntrparam level incremental 0 1 0 2 20 Uniform levels from 0 1 to 20 gnuplot gt set view map gnuplot gt unset surface gnuplot gt unset key gnuplot gt splot data 0001 tab u 1 2 3 w lines 2 http www exelisvis com 3https wci lInl gov codes visithome html http www paraview org gt http www wolfram com Shttp www gnuplot info Version 4 2 or higher is recommended CHAPTER 11 OUTPUT AND VISUALIZATION 99 Density Pressure X Figure 11 1 Density and pressure plots for the Sod shock tube using Gnuplot Figure 11 2 Density and pressure plots for the Sod shock tube using Gnuplot Binary Data Fi
94. cal composition X Also e 1 must be monotonically increasing CHAPTER 7 EQUATION OF STATE 69 needed to evaluate the sound speed cs yI p p Note that T has the upper bound of 5 3 and may not be straightforward to compute Fortunately its value is only needed to estimate the wave propagation speed during the Riemann solver and an approximate value should suffice Two different implementations are provided with the current distribution pvte_law_H c is suitable for a partially hydrogen gas in LTE described in the next section while pvte_law_dAngelo c can be used for molecular and atomic hydrogen cooling as in D Angelo G et al ApJ 2013 778 More technical details can be found under the Src EOS PVTE folder in the API reference guide or following this link Note The PVTE_LAW EOS is not compatible with algorithms requiring characteristic decomposi tion and cannot be used with the ENTROPY_SWITCH We suggest to use RK time stepping and the tvdlf hll or hllc Riemann solvers This EOS is at present available for the HD and MHD modules only 7 3 1 Example EOS for a Partially Ionized Hydrogen Gas in LTE As a simple non trivial example consider a partially ionized hydrogen gas in Local Thermodynamic Equilibrium LTE no cooling see also 32 4 of VMBMI15 Let the particle number densities be No neutrals Np Ne charge neutrality gt n le 28 Density and pressure can then be written as P MpNp Mp
95. e code and work in all system of coordinates and physics mod ules The dimensionally unsplit version is stable up to CFL lt 1 Naim Where Naim is the number of dimensions Similar to the previous setup but it has better stability proper ties for higher than and order interpolants The dimensionally unsplit version is stable up to CFL lt 1 Naim MUSCL Hancock second order scheme of van79 Tor97 Com putationally more efficient than RK integrators it is prob ably the fastest 2 order algorithm Works well for the HD RHD modules and the MHD RMHD modules with DIVERGENCE_CLEANING or EIGHT_WAVES particularly on Cartesian 1 2 3 dimensions or cylindrical geometries More sophisticated than the previous one it yields the Piecewise Linear Method of Col85 Gives the original Piecewise Parabolic Method PPM of CW84 Suggested for HD RHD or MHD with DIVERGENCE_CLEANING or EIGHT_WAVES in Cartesian 1 2 3 dimensions or cylindrical geometries It is stable up to CFL lt 1 and it has small dissipa tion Yields the Corner Transport Upwind method of and for the HD RHD or MHD RMHD modules the RMHD version works only with HANCOCK It is fully un split and stable up to 1 in 2 D and 0 5 in 3D It is one of the most sophisticated algorithms available It is suitable for com putations in Cartesian and cylindrical grids in the HD RHD and MHD module Table 2 2 Suggested algorithm configurations The cost 4th column is
96. e member box gt vpos which specifies the location of the variable inside the cell CENTER X1FACE X2FACE X3FACE provide an efficient way of looping through the ghost boundary zones using the macro BOX_LOOP box k j i which automatically takes care of the index range of definition Note Using the box structure is not strictly mandatory and the usual macros X1_BEG_LOOP X3_END_LOOP may still be employed without any modifications However these macros perform loops over cell centered data stencils and staggered field are not completely defined since the loops do not include one row of zones at the furthest left edges of the boundary zones On the contrary the BOX_LOOP macro takes into account the full range of definition of the variable and should be used whenever possible e side an input integer label specifying on which side of the physical domain user defined values should be prescribed It can take on the following values X1_BEG X1_END boundary conditions can be assigned in the ghost zones at the beginning and end of the physical domain in the x direction CHAPTER 5 INITIAL AND BOUNDARY CONDITIONS INIT C 50 Interior values m X Face interior values AY Face interior values Boundary values O X Face boundary values AY Face boundary values l l l l l l 1 Cell Centered Variables Face Staaaered Variables Figure 5 1 Schematic representation of cell centered left panel and face centered
97. e CONST_au 49597892e13 xx lt Astronomical unit x define CONST_c 99792458e10 lt Speed of Light x define CONST_eV 602176463158e 12 x lt Electron Volt in erg define CONST_G 6726e 8 xx lt Gravitational Constant x define CONST_h 62606876e 27 x lt Planck Constant as define CONST_kB 3806505e 16 lt Boltzmann constant kf define CONST_ly 9461e18 xx lt Light year x define CONST_mp 67262171e 24 e Proton mass x7 define CONST_mn 67492728e 24 xx lt Neutron mass x define CONST_me 1093826e 28 x lt Electron mass x define CONST_mH 6733e 24 xx lt Hydrogen atom mass x define CONST_Msun e33 x lt Solar Mass x define CONST_Mearth 9736e27 xx lt Earth Mass x define CONST_NA 0221367e23 lt Avogadro Contant ef define CONST_pc 0856775807e18 x x lt Parsec xy define CONST_PI 14159265358979 fe XES ipi Aris T define CONST_Rearth 378136e8 lt Earth Radius x define CONST_Rsun 96e10 lt Solar Radius x define CONST_sigma 67051e 5 lt Stephan Boltmann constant define CONST_sigmaT 6524e 25 fex lt Thomson Cross section x DONA WWAUONFOFF DEFRA VDEFENEF EF B 2 Predefined Function Like Macros PLUTO comes with a number of pre defined function like macros to implement simple arithmetic op erations such as maximum MAX minimum MIN or looping over a specific portion of the computa tional domain e g DOM_LOOP or TOT_LOOP Please re
98. e divergence free condition in Godunov type schemes for ideal magne tohydrodynamics the upwind constrained transport method Journal of Computational Physics 195 2004 17 48 M S Liou A Sequel to AUSM AUSM Journal of Computational Physics 129 1996 364 382 L D Landau and E M Lifshitz Fluid mechanics 2 ed Pergamon Press Oxford 1987 X D Liu S Osher and T Chan Weighted Essentially Non oscillatory Schemes Journal of Computational Physics 115 1994 200 212 R Lohner An adaptive finite element scheme for transient problems in CFD Computer Methods in Applied Mechanics and Engineering 61 1987 323 338 A Mignone and G Bodo An HLLC Riemann solver for relativistic flows I Hydrodynamics MNRAS364 2005 126 136 A Mignone G Bodo S Massaglia T Matsakos O Tesileanu C Zanni and A Ferrari PLUTO A Numerical Code for Computational Astrophysics ApJS170 2007 228 242 A Mignone M Flock M Stute 5 M Kolb and G Muscianisi A conservative orbital advection scheme for simulations of magnetized shear flows with the PLUTO code A amp A545 2012 A152 A Mignone A simple and accurate Riemann solver for isothermal MHD Journal of Computational Physics 225 2007 1427 1441 High order conservative reconstruction schemes for finite volume methods in cylindrical and spherical coordinates Journal of Computational Physics 270 2014 784 814 T Miyoshi and K Kusano A multi state HLL appr
99. e expression above holds for an isotropic viscous stress and the resulting tensor is symmetric The actual diffusion terms will then be given by V M and V v M on the right hand side of the mo mentum and the energy equation respectively These fluxes are added to the hyperbolic momentum and energy fluxes computed with a Riemann solver in a fully conservative and unsplit fashion if EXPLICIT is chosen In curvilinear geometries additional geometrical source terms coming from the tensor s di vergence are added to the right hand side of the equations On the other hand if VISCOSITY is set to STS advection and diffusion terms are treated separately using operator splitting The implementation of the previous expressions together with the equation module can be found under the directory Src Vis cosity Derivative terms are discretized at cell interfaces using second order accurate finite differences and assuming a uniform grid spacing Note that when using FARGO MHD this module can operate only with STS The viscous transport coefficients v shear and v2 bulk are defined in the function Visc_nu in the source file PLUTO Src Viscosity visc_nu c This file should be copied from its original folder to the actual working directory before doing any modification 72 CHAPTER 8 DISSIPATIVE EFFECTS 73 The Visc_nu function has the following syntax Syntax void visc nu double v double xl double x2 double x3 double nul double szn
100. ear shear velocity dlog Q R q th A A 10 1 O a ER a 10 1 where 29 is the local constant angular velocity and q is a local measure of the differential rotation q 3 2 for a Keplerian profile The module solves the HD or MHD equations in a non inertial frame so that the momentum and energy equations become a V puv BB Vp pg 2Q02 x pu Ot 95 10 2 3 FW E p v w B B pvu 9s where g Mi 2qxt 22 is the tidal expansion of the effective gravity while the second term in Eq 10 2 represents the Coriolis force The continuity and induction equations retain the same form as the original system While the computational box should be periodic in the azimuthal y direction radial x boundary conditions are determined by image boxes sliding with relative velocity w q 0L relative to the computational domain Fig 10 1 In other words the boundary conditions at the left right x boundaries are dut dge ELTU l 403 OLD U2 Opec leg Uo ew where q is any other flow quantities except vy The ShearingBox module is implemented inside Src MHD ShearingBox and works at present with the HD equations or with CONSTRAINED_TRANSPORT MHD Parallelization can be performed in all three spatial dimensions 86 CHAPTER 10 ADDITIONAL MODULES 87 10 1 1 Using the module The shearingbox module is enabled by invoking the Python setup script with the with sb option It is compatible with the ISOTHE
101. ed arguments on the value taken by DIMENSIONS Likewise the macro EXPAND a b c is used to write component independent code depending on the value taken by COMPONENTS Function like macro are documented in the file macro h of the API reference guide see Doc Doxygen html macros_8h html Index Cylindrical Polar Spherical Quantity Physics Module RHO rest mass density ALL VX1 iVR iVR iVR x velocity ALL VX2 iVZ iVPHI iVTH x2 velocity ALL VX3 iVPHI iVZ iVPHI x3 velocity ALL PRS thermal pressure ALL BX1 BR BR BR x1 cell centered magnetic field MHD RMHD BX2 iBZ iBPHI iBTH x2 cell centered magnetic field MHD RMHD BX3 iBPHI iBZ iBPHI x3 cell centered magnetic field MHD RMHD BX1s iBRs iBRs iBRs x staggered magnetic field MHD RMHD BX2s iBZs iBPHIs iBTHs x2 staggered magnetic field MHD RMHD BX3s iBPHIs iBZs iBPHIs x3 staggered magnetic field MHD RMHD TRC tracer passive scalar Q ALL hey S Table 5 1 Array indices used for labeling primitive variables Staggered components fields in the boundary conditions see 46 2 3 3 suffix are used only for magnetic CHAPTER 5 INITIAL AND BOUNDARY CONDITIONS INIT C 44 5 1 1 Units and Dimensions In general PLUTO works with non dimensional or code units so that flow quantities can be properly scaled to reasonable numbers Although it is possible in principle to work directly in c g s units i e cm sec and gr we strongly recommend to
102. esolu tion in the other coordinate directions accordingly 16 work w AMR Yes Yes No Yes Yes Yes No Yes Table 1 3 Command line options available when running PLUTO Compatibility with AMR version is given in the last column y on parallel architectures only 1 5 Modifying the Distribution Source Files PLUTO source files are compiled directly from the PLUTO Src directory Should you need to modify a C source file other than your init c we strongly advise to copy the file to your local working directory and then edit it since the latter is always searched before PLUTO Src during the compilation phase In other words if you want to modify say boundary c copy the file to your working area and introduce the appropriate changes When make is invoked your local copy of boundary c is compiled since it has priority over PLUTO Src boundary c which is actually ignored In such a way you can keep track of the problem dependent modification without affecting the original distribution Note Header files h or H do not follow the same convention and must not be copied to the local working directory Modifications to header files must therefore be done in the original directory 2 Problem Header File definitions h This chapter explains how to create the configuration header file definitions h for a specific problem 2 1 Basic Options The header file definitions h is created by the Python script setup py by
103. fault are set to be equal to solar abundances X 0 711 Y 0 2741 see also Appendix B 3 e by calling the GetMu function when the PVTE_LAW equation of state is adopted with equili brum ionization see In this case temperature in Kelvin and density in code units must be supplied as input arguments In Kelvin In code units Means rho UNIT_DENSITY in cgs units mu Example 3 Consider a flow with typical number densities of the order of n 4 x 10 cm temperature T 2 5 x 10 K and Mach number M v c 15 while the typical length scale is Ly 10 AU Suppose also that a magnetic field with strength of 10 uG is also present Units can be set in definitions h x Beg user defined constants do not change this line x define UNIT_DENSITY 1 e3 CONST_mp define UNIT LENGTH 10 0 CONST_au define UNIT VELOCITY 1 65 x End user defined constants do not change this line x The sound speed c is defined for an adiabatic equation of state by the relation Pp TT Cs AA p Ku The initial condition is then implemented as follows v RHO 4 0 means 4 1073 1 67e 24 qgr cm 3 T 2 5e3 Temperature in Kelvin Hif COOLING SNEq CompEquil n T v Compute ionization fraction at equilibrium endif mu MeanMolecularWeight v v PRS v RHO T KELVIN xmu Pressure in units of rho0xv0 2 x Define sound speed and velocity lt c
104. fer to the Doc Doxygen html macros_8h html page in the API refernece guide Doc Doxygen html index html 119 B 3 Advanced Options PLUTO allows a number of switches to be fine tuned directly from your definitions h in the user defined constant section see 42 3 These advanced options are described in Table B 1 Table B 1 PLUTO advanced options Name Value Description The amount of artificial Lapidus type viscosity y added to the two shock Riemann solver flux only Fr gt Fu max Uun L Un R ONU Up B 1 where vn is the velocity normal to the interface This term intro duces an extra amount of numerical dissipation CW84 useful to reduce small amplitude oscillations occurring when a char acteristic speed associatedwith a strong shock measured rela tive to the grid vanishes Typical values are around 0 1 By default it is not used Note this constant has no effect with other Riemann solvers ARTIFICIAL_VISC real In curvilinear coordinates set this switch to YES to enforce an gular momentum conservation during the prolongation and re striction operation with PLUTO Chombo Default value is YES when CHOMBO_EN_SWITCH is set to YES Enable this switch if you want to produce an equally spaced logarithmic grid in the radial direction in POLAR or SPHERICAL CHOMBO_LOGR YES NO coordinates when using PLUTO Chombo A logarithmic grid has the advantage of preserving the cell aspect ratio both close to and far a
105. ff temperature e gmaxCoolingRate limit the time step so that the maximum fractional thermal losses cannot exceed yg maxCoolingRate In general 0 lt g_maxCoolingRate lt 1 the default is 0 1 e gminCoolingTemp sets the cut off temperature below which cooling is artificially set to 0 79 CHAPTER 9 OPTICALLY THIN COOLING 80 9 1 Power Law Cooling Power law cooling is the most simple form of cooling where the loss term in the internal energy equation becomes A ap T 9 2 There are no new species when this form of cooling is selected When an ideal equation of state is used the source step becomes dp 2 a a de L aa Kn dt and since density is not affected during this step integration is done analytically 1 pn SAS a 7 for aX pS 9 3 p exp CAt for a 1 where C T 1 a p Ky is a constant The default power law accounts for bremsstrahlung cooling by solving dPcgs T 1 a Des T K dp o 2 9 4 dt oes U M dt p with p t and p given in code units and Frei poLo kgumg v CG ti where po vo and Lo are the reference density velocity and length defined in 45 1 1Jand ap 2 107 in expressed in c g s units The implementation of this cooling step with a 1 2 can be found under Src Cooling Power_Law cooling c CHAPTER 9 OPTICALLY THIN COOLING 81 9 2 Tabulated Cooling The tabulated cooling module provides a way to solve the internal energy equatio
106. file is then automatically created by the Python script by dumping Chombo makefile variables into the file make vars part of your local working direc tory Although system dependencies have already been created during the Chombo compilation stage the Change makefile option from the Python menu is still used to specify the name and flags of the C compiler used to compile PLUTO source files This step is achieved as usual by selecting a suitable defs file from the Config directory see Chapter 3 Beware that during this step additional variables such as PARALLEL USE_HDF5 etc normally used in the static grid version have no effect since Chombo has its own independent parallelization strategy and I O Fortran and C compilers are the same ones used to build the library Initial and boundary conditions are coded in the usual way but definitions h and pluto ini may contain some AMR specific directives 12 2 1 Running PLUTO Chombo Once PLUTO Chombo has been compiled and the executable pluto has been created PLUTO runs in the same way i e MyWorkDir gt pluto flags where the supported command line options are given in Table 1 3 in Note that restart must be followed by the restart checkpoint file number An error will occur otherwise Parallel runs proceeds in the usual way e g MyWorkDir gt mpirun np 8 pluto flags Note that when running in parallel each processor redirects the output on a separate file pout n
107. formation stored in out files We thus strongly recommend to backup these files together with the actual datafiles Note Restart is possible only using the dbl or dbl h5 data formats 11 1 1 Binary Output dbl or flt data formats Binary data can be dumped to disk at a given time step as i one single file containing all variables by selecting single_file in pluto ini or ii as a set of separate files for each variable multiple_files We recommend the second option for large data sets The base name is set to data for a single data file containing all of the fields or takes the name of the corresponding variable if multiple sets are preferred see Table 11 1 Restart can be performed from double precision binary data files by invoking PLUTO with the restart n command line option where n is the output file number from which to restart In this case an additional file restart out will be dumped to disk The corresponding log file dbl out or flt out is a multi column ascii files of the form nout E dt nstep single_file little varli var2 where nout t dt and nstep are respectively the file number time time step and integration step at the time of writing The next column single_file multiple_files tells whether a single file or multiple files are expected The following one little big gives the endianity of the architecture whereas the remaining columns list the variable names and their order in this particular format Def
108. fusion terms Some of the C functions normally used in the static grid version of PLUTO have been replaced by C codes in order to interface the structure of PLUTO with the Chombo library For instance the main function main c has been replaced by amrPluto cpp 12 1 Installation In order to properly install PLUTO Chombo you will need check also Table 1 1 e C C and Fortran compilers e the MPI library for parallel runs e GNU make e the Chombo library version 3 2 is recommended e the HDF5 library version lt 1 8 14 is recommended e the chombo 3 2 patch tar gz provided with the PLUTO distribution which replaces some of the library source files The following sections give a quick headstart on how these libraries can be built for being used by PLUTO Please consult the libraries respective documentation for additional information lhttps commons Ibl gov display chombo 2http www hafgroup org HDF5 106 CHAPTER 12 ADAPTIVE MESH REFINEMENT AMR 107 12 1 1 Installing HDF5 The HDF5 library can be downloaded from http www hdfgroup org HDF5 and it can be used for the static grid version 411 1 2 while it is mandatory for the AMR version and it must be installed before compiling Chombo Note Both PLUTO static and PLUTO Chombo AMR have been succesffuly tested for serial and parallel computation using with HDF5 v lt 1 8 14 while parallel I O problems were found on Ubuntu systems using newer ver
109. given in terms of number of Riemann problems per cell per step Naim is the number of spatial dimensions ChTr stands for CHARACTERISTIC_TRACING CHAPTER 2 PROBLEM HEADER FILE DEFINITIONS H 22 2 1 9 NTRACER The number of passive scalars or colors denoted with Q obeying simple advection equations of the form 90 IQ OX k l AP Ot WU VO 0 At The second form gives the conservative equation and it is the one actually being discretized by the code The array index used to access tracer variables 45 1 45 2 is TRC for the first tracer TRC 1 for the second one and so on The maximum number is 4 V pQzv 0 2 2 2 1 10 USER DEF PARAMETERS 600 gt gt User defined Parameters lt lt FOO_1 FOO_2 NAME or VALUE gt F00 38 Figure 2 2 User defined pa rameter names are chosen in this sub menu Sets the number of user defined parameters that can be changed at runtime and accessed from any where in the code The explicit numerical value is read at runtime from pluto ini and can be changed before execution without re compiling the code The parameters are identified by means of a label corresponding to an index of the global array g_inputParam visible anywhere in the program If for instance USER_DEF_PARAMETERS has been set equal to 3 you will be prompted to define 3 differ ent labels say FOO_1 FOO_2 and FOO_3 as in Fig These names are the integer indexes of the g inputParam array g_inputParam FO
110. he internal energy equation modifies to Op 2 ap 1 U Vet pc V v P 1D 9 3 d 8 4 The resistive tensor is assumed to be diagonal with components 1 diag Mais 1x2 x3 8 5 The module is implemented in the Src MHD Resistivity directory and the functional form of 7 can be specified by editing your local copy of PLUTO Src MHD Resistivity res_eta c which includes the function Resistive_eta Syntax void Resistive _eta double xv double xl double x2 double x3 double xJ double xeta Arguments e v a pointer to a vector of primitive variables e x1 x2 x3 local spatial coordinates e J a pointer to the electric current vector e eta a pointer to an array containing the three components of the resistive diagonal tensor The resistive module works in 1 2 and 3 dimensions in all systems of coordinates on both uniform and non uniform grid although higher accuracy can be achieved on uniform grid spacing Both cell centered and staggered MHD are supported using either EXPLICIT or SUPER_TIME_STEPPING inte gration The choice of the V B 0 control strategy determines two different ways to compute V x B For cell centered MHD the three components of the current are calculated at the zone interfaces normal to the sweep integration direction For staggered MHD current components are computed just once at the beginning of the step at cell edges using the staggered field and the three components of J have different spatial locat
111. hould contain the full path name to your HDF5 library e g usr local lib HDF5 serial lib in the example above Please make sure to add for example gt setenv LD LIBRARY PATH usr local lib HDF5 serial lib SLD LIBRARY PATH to your tcshrc if you re using the tcsh shell or gt export LD_LIBRARY_PATH usr local 1lib HDF5 serial lib SLD_LIBRARY_PATH if you re using bash If you do not want shared libraries then add disable sharedto the configure command 12 1 2 Installing and Configuring Chombo Chombo 3 2 can be downloaded by direct access to the SVN server repository after free registration The Chombo source code distribution should be unpacked under PLUTO Lib and some of the library source files must be replaced using the chombo 3 2 patch tar gz patch archive provided with the PLUTO distribution A typical session is gt get the 3 2 release of Chombo gt svn username username co https anag repo lbl gov svn Chombo release 3 2 Chombo 3 2 gt tar xzvf chombo 3 2 patch tar gz C Chombo 3 2 apply PLUTO Patch CHAPTER 12 ADAPTIVE MESH REFINEMENT AMR 108 In order to use Chombo you may have to build different libraries depending on the chosen compiler serial parallel build number of dimensions optimizations etc If you intend to run PLUTO Chombo for serial or parallel computations in one two or three dimensions we suggest to compile all possible configurations that is 1 2 and 3D serial or 1 2 and 3D parallel Li
112. hs and so on Important for creating the makefile e Doc documentation directory e Lib repository for additional libraries e Src main repository for all c source files with the exception of the init c file which is left to the user The physics module source files are located in their respective sub directories HD classi cal hydrodynamics RHD special relativistic hydrodynamics MHD magnetohydrodynamics RMHD relativistic magnetohydrodynamics Cooling viscosity thermal conduction and addi tional physics models are located under the folders with similar names e g Cooling Viscosity Thermal_Conduction The Templates directory contains templates for the user dependent files such as init c pluto ini makefile and definitions h e Tools Collection of useful tools such as Python scripts IDL visualization routines pyPLUTO etc e Test_Problems a directory containing several test problems used for code verification PLUTO should be compiled and executed in a separate working directory which may be anywhere on your local hard drive Although most of the current algorithms can be considered in their final stable version the code is under constant development and updates are released once or twice per year When upgrading to a newer version of PLUTO it is recommended that the entire PLUTO directory tree be deleted Syntax changes are usually listed in the file CHANGES in the PLUTO root directory CHAPTER
113. ible and efficient I O HDF5 data can be visualized by a number of commercial or open source packages and at present Chombo data files have been successfully opened and visualized with ID VisIf ParaView and pyPLUTO Examples are provided in the following A comprehensive list of application software using HDF5 may be found at http www hdfgroup org toolsSapp html A set of utilities for manipulat ing visualizing and converting HDF5 data files is provided by Hbutils a set of utilities available at http www hdfgroup org products hdf5_tools H5utils offers a simple tool for batch visualization as PNG images and also includes programs to convert HDF5 datasets into the formats required by other free visualization software e g plain text Vis5d and VTK In what follows we describe some of the routines provided with PLUTO Chombo for viewing and analyzing HDF5 data in the IDL programming language 12 4 1 Visualization with IDL PLUTO Chombo comes with a set of visualization routines for the IDL programming language For more information consult idl_tools html 3http www exelisvis com https wci In gov codes visit home html gt http www paraview org CHAPTER 12 ADAPTIVE MESH REFINEMENT AMR 111 The procedure HDF5LOAD located in Tools IDL hdf5load pro can read a HDE5 data file and store its content on the usual set of variables used during a typical IDL session HDF5LOAD is directly called from PLOAD 411 3 2 when the latter
114. ic field in the usual way by directly prescribing the values for us BX1 us BX2 and us BX3 2 When ASSIGN_VECTOR_POTENTIAL is set to YES the vector potential A is used instead and the magnetic field is recovered from B V x A This option guarantees that the initial field has zero divergence in the discretization which is more appropriate for the underlying formulation i e cell or face centered fields 46 2 3 Note In 2D only the third component of A that is us AX3 should be used Likewise the third component of the magnetic field B cannot be assigned through the vector potential and must be prescribed in the standard way see the third example in Table 6 1 Table 6 1 shows some examples of magnetic field initializations with and without using the vector potential Magnetic Field Standard Using Vector Potential us BX1 0 0 us AX1 0 0 A us BX2 5 0 us AxX217 0 0 era us BX3 0 0 us AX3 x1x5 0 B 0 5 0 us BX1 0 0 us AX1 0 0 Cvlindrical 2D us BX2 5 0 us AX2 0 0 y l wera 0 05 asla 0 5 15 0 BX1 sin x2 us AX1 0 0 B sin y sin 2x 2 es B da us AX2 0 0 Cartesian 2 5D a T SED NS EEAST Fs TAX3 cos x2 0 5 cos 2 0 x1 us BX3 2 0 us BX3 2 0 Table 6 1 Examples of how the magnetic field may be initialized Direct initialization standard is possible when ASSIGN_VECTOR_POTENTIAL is set to NO Otherwise the components of the vector potential
115. ics module 4 4 The Solver Block e Solver string Sets the Riemann solver for flux computation From the most accurate i e least diffusive to the least accurate i e most diffusive two_shock The Riemann problem is solved exactly or approximately depending on the particular solver implemented for a given physics module at every interface this is usually more accurate but computationally intensive See CW84 for the HD module and MB05 for the relativistic hydro equations roe Linearized Roe Riemann solver based on characteristic decomposition of the Roe matrix Roe8 1 ausm Advection Upstream Splitting Method of Lio96 only for the HD module hild The hlld approximate Riemann solver of UKO05 for the adiabatic case Mig07 for the isothermal case and MUBO9 for the relativistic MHD equations hiic Harten Lax Van Leer approximate Riemann Solver that restores with the middle con tact discontinuity h11 Harten Lax Van Leer approximate Riemann Solver tvdlf A simple Lax Friedrichs scheme is used Note that not all solvers are available for a given physics module see Table We warn the user that under some circumstances high Mach number flows low density plasmas more dif fusive solvers such as HLL or TVDLF turn out to be more robust than accurate solvers However hybrid adaptive strategies can be turned on when SHOCK FLATTENING is set to MULTID CHAPTER 4 RUNTIME INTITALIZATION FILE P
116. ing TABULATED The HD and MHD modules share the same RightHandSide function inside MHD rhs c Likewise source terms have been separately implemented in MHD rhs_source c Runtime information can now be accessed from anywhere in the code using the RuntimeGet function which gives access to a copy of the runtime structure e g RuntimeGet gt cfl RMHD module allows faster computation of eigenvalues using approximate expressions RMHD_FAST_EIGENVALUES Added the macro RECONSTRUCT _4VEL YES NO which allows to reconstruct the 4 velocity rather than the 3 velocity It works for for RHD and RMHD Fixed Bugs IDL pload function with var now works also for multiple files Fixed a number of minor bugs in H2_COOL module and improved CompEquil function PPM method works correctly for non uniform grid in any coordinate system FlagShock now called when using chombo MULTID shock flattening and ENTROPY SWLECH PLUTO Chombo now writes plotfiles correctly when restarting from a checkpoint file changed AMR cpp in library GLM complained with Background field splitting print removed Resistivity background can be used removed QUIT statement from MHDY bckornd fi le E PathCTU 1D has been removed and all 1D AMR computations are done With PatchCTU cpp or Patchhuler cpp VTK files can now be visualized with PARAVIEW PVTE_LAW now uses cubic spline interpolation rather than linea
117. ing on the specific configuration you intend to use The minimal set to get PLUTO running on a workstation with a static grid no AMR requires Python a C compiler and the make utility These are usually installed by default on most Linux Unix platforms A comprehensive list is shown in Table 1 1 Adaptive Grid serial parallel parallel Python gt 2 0 yes yes yes es Ces compiler yes yes gt Fortran compiler yes yes GNU make yes yes yes yes PMPIMibrary ys yes Chombo library 432 yes yes A PNG library opt opt Table 1 1 Software requirements for different applications of PLUTO Here opt stands for optional serial refers to single q pp pP p 8 processor runs and parallel to multiple processor architectures Starting with PLUTO 4 parallelization is handled internally and ArrayLib used in previous versions of the code is no longer necessary The Chombo library is required for computations making use of Adaptive Mesh Refinement Chapter 12 while the PNG library should be installed only if PNG output is desired The HDF5 library is required for I O with the Chombo library and may also be used with the static grid version of the code 1 2 Directory Structure Once unpacked your PLUTO root directory should contain the following folders e Config contains machine architecture dependent files such as information about C compiler flags library pat
118. ion Note The resistive module is only partially compatible with the entropy switch 92 2 4 CHAPTER 8 DISSIPATIVE EFFECTS 75 8 3 Thermal Conduction Thermal conduction can be included for the hydro HD or MHD equations by introducing an additional divergence term in the energy equation OE a TV E pi v B Ww B V Fe 8 6 where F is a flux limited expression that smoothly varies between the classical and saturated thermal conduction regimes F class and Fiat respectively F sat Pec gt Fiat F F ass F sizes 8 7 see Spi62 OBR 08 In the MHD case thermal conductivity is highly anisotropic being largely suppressed in the direction transverse to the magnetic field Denoting with b B B the unit vector in the direction of magnetic field the classical thermal conduction flux may be written as Bal86 Pas y b b vr ite VI b 6 vr 8 8 where the subscripts and denote respectively the parallel and normal components to the magnetic field T is the temperature and k are the thermal conduction coefficients along and across the field In the purely hydrodynamical limit no magnetic field Eq reduces to Fe Ky VT Saturated effects are accounted for by making the flux independent of VT for very large tempera ture gradients Spi62 CM 77 In this limit the flux magnitude approaches Fsat 5 pc where is the isothermal speed of sound and lt 1 is a free parameter Note however that
119. ion YE Hell 0 23 x for example x in your Init function The fractions of all ion species can also be automatically set for equilibrium conditions using the CompEquil function in Src Cooling MINEq comp_edquil c double CompEquil double N double T double xv where N and T are the plasma number density and temperature respectively and v is a vector of primitive variables The function will return the electron density as output and v will contain the computed ionization fractions the other variables are not affected The routine solves the system of equations for abundances in equilibrium Note The number of ions for C N O Ne and S may be reduced from 5 to a lower number gt 1 by editing Src Cooling MINEq cooling h This may reduce computational time if the expected temperatures are not large enough to produce high ionization stages e g IV or Vif T lt 10 K The current default value is 3 The elements abundances are set in radiat c from the Src Cooling MINEq folder When using the MINEq module the cooling coefficients tables are generated at the beginning of the simulation by the routines in Src Cooling MINEq make_tables c Update or customization of the atomic data can be done by editing this file The ion fractions are integrated through advection equations of the form OX ay v VX Si 9 13 where the source term S is computed taking into account collisional ionization radiative and dielec
120. ions and can be safely used at restart since the last position of the file is automatically searched for and subsequent writing is appended starting from the correct row CHAPTER 5 INITIAL AND BOUNDARY CONDITIONS INIT C void Analysis const Data xd Grid grid int ip J K double dV vol scrh double Ekin Eth_max vx2 vy2 vz2 double x xdx xdy dz fe Bet pointer shortcuts dx grid IDIR dx dy grid JDIR dx dz grid KRDIR ax x Main loop x Ekin Eth_max 0 0 DOM_LOOP k j i dav dx i dy j dz k Cell volume Cartesian coordinates x vx2 d gt Vc VX1 k jJ i d gt Vc VX1 kK j i x velocity squared x vy2 d gt Vc VX2 k 3111 d YelVvxX2 11K1191 1141 y velocity squared vz2 d gt Vc VX3 k jJ i d gt Vc VX3 kK 91 11 z velocity squared scrh 0 5xd gt Vc RHO k J 11 vXx2 vy2 vz2 cell kinetic energy f Ekin scrhxdvV scrh d gt Vc PRS k 3 i1 g_gamma 1 0 cell internal energy Eth_max MAX Eth_max scrh vol g_domEnd IDIR g_domBeg IDIR Compute total domain volume x vol g_domEnd JDIR g_domBeg JDIR vol g_domEnd KDIR g_domBeg KDIR Ekin vol Compute kinetic energy average Parallel data reduction Hifdef PARALLEL MPI Allreduce amp Ekin amp scrh 1 MPI_DOUBLE MP1_SUM MPI_COMM WORLD Ekin scrh MPT Allreduce amp Eth_max
121. ir is omitted the current directory is used Not to be confused with output dir in e Checkpoint_interval double integer string Assigns the output interval s for checkpoint restart files The first field double can be used to set the time interval or period in code units The second field can be either an integer ora string An integer value determines the number of steps between consecutive outputs Alternatively a string value can be used to set the wall clock time between successive out puts a value of 3 55h for instance means that a checkpoint file is saved every 3 hours and 55 minutes Checkpoints files contain conservative variables e Plot_interval double integer Sets the output frequency for plot data files The meaning of the fields is the same used for Checkpoint_interval except that no wall clock interval is permitted Output files are stored using the HDF5 file format and numbered as data nnnn hdf5 where n is a zero padded sequentially increasing integer Data files contain primitive variables Note that a negative number means that checkpoint or plot files are never written 0 means that checkpoint files are written before the initial timestep and after the final one CHAPTER 4 RUNTIME INTITALIZATION FILE PLUTO INI 41 4 8 The Parameters Block e PAR _NAME_1 double e PAR_NAME_n double User defined parameter values are read at runtime in this section The labels on the left
122. it c o The shear parameter and orbital frequency are no longer assigned by the global variables sb_q and qb_Omega but using the macros SB_Q and SB_OMEGA which by default take the values 1 5 and 1 0 Userguide has been re organized Added the MeanMolecularWeight function and gas composition see userguide The global variable g_stepNumber now works with PLUTO Chombo Chombo HDF5 files can be written by specifying the output interval in clock time hours minutes seconds as well as the output directory Added GLM_ ALPHA to the list of user defined constants Runtime structure renamed to Restart Input structure renamed to Runtime and setup c renamed to runtime_setup c Suppressed replaced constant names ARTIFICIAL VISCOSITY gt ARTIFICIAL VISC see appendix tn the userguide o CHOMBO_EN SWITCH suppressed replaced by the new ENTROPY_SWITCH O CH TRACING REP STATE gt GHTR REF STATE o DN PR are no longer supported gt use RHO and PRS instead o EGLM gt GLM_EXTENDED o INTERPOLATION gt RECONSTRUCTION o MHD_FORMULATION gt DIVE CONTROL o RESISTIVE_MHD gt gt RESLTS Liv ry o STS_nu gt STS_NU o USE_FOUR_VELOCITY gt RECONSTRUCT_4VEL works for both RHD and RMHD RHD primitive variables always contain the 3 velocity although 4 vel can CHAPTER 0 QUICK START 11 be reconstructed by enabling the RECONSTRUCT_4VEL switch Src Cooling Tab changed to Src Cool
123. ity or the x component of velocity are written with scalar attribute as they are However by setting VIK_VECTOR_DUMP to YES vector fields such as velocity and magnetic field can be saved with the vector attribute and their components are automatically transformed to Cartesian See also Table B 1 If a VTK file is written to disk the log file vtk out is updated in the same manner as dbl out or flt out Default By default all variables except staggered magnetic field components if any are written 11 1 4 ASCII Output tab Data format The tab format may be used for one dimensional data or relatively small two dimensional arrays in serial mode only We warn that this output is not supported in parallel mode The output consists in multi column ascii files named data nnnn tab of the form x 1 y 3 varl i J var2 i Jj vars ipj sss where the index j changes faster and a blank records separates blocks with different index Default By default all variables except staggered magnetic field components if any are written 11 1 5 Graphic Output ppm and png data formats PLUTO allows to take two dimensional slices in the 1112 1173 Or 212123 planes and save the results as color ppm or png images The Portable Pixmap ppm format is quite inefficient and redundant al though easy to write on any platform since it does not require additional libraries The Portable Network Graphics png is a bitmap image format that employ
124. ive or diffusion terms which at present include e Viscosity HD MHD described in e Resistivity MHD described in e Thermal conduction HD MHD described in Each modules can be individually turned on from the physics sub menus accessible via the Python script Numerical integration of diffusion processes viscosity resistivity and thermal conduction requires the solution of mixed hyperbolic parabolic partial differential equations which can be carried out using either a standard explicit time stepping scheme or the Super Time Stepping STS technique see Depending on the time step restriction you may include diffusion processes by setting the correspond ing sub menu choice s to EXPLICIT or to SUPER_TIME_STEPPING respectively 8 1 Viscosity The viscous stresses enter the HD and MHD equations with two parabolic diffusion terms in the mo mentum and the energy conservation Adding the viscous stress tensor to the original conservation law Eq 1 1 we obtain a mixed hyperbolic parabolic system which in compact form may be expressed by the following oU ae e 8 1 where I represents the viscous stress tensor whose components are given by a Aa U 2 m be gt va qn V VOsj 8 2 Coefficients v and 1 are the first shear and second bulk parameter of viscosity respectively v and v denote the covariant derivatives whereas h hj are the geometrical elements of the respective direc tion Th
125. ive description In the adaptive grid version the code relies on the Chombo library for adaptive mesh refinement AMR written in C and Fortran Chapter 12 A detailed description of the AMR implementation is given in MZ T 121 Doxygen is used as the standard documentation system and the Application Programming Interface API reference guide can be found in Doc API ReferenceGuide html PLUTO has been successfully ported to several parallel platforms including Linux Windows Cyg win Mac OS X Beowulf clusters IBM power4 power power6 SGI Irix IBM BluGene P and several others Figure 1 1 shows the strong scaling on a BlueGene P machine up to 32 768 processors on a pe riodic domain with 512 computational grid zones Time to Solution Speedup Time to Sol Speed l i i me deal F deal ZL Figure 1 1 Strong scaling of PLUTO on a periodic domain problem with 512 grid zones Left panel average execution time in seconds per step vs number of processors Right panel speedup factor computed as 11 Ty where T is the inferred execution time of the sequential algorithm and Ty is the execution time achieved with N processors Code execution time is given by black circles dotted line while the solid line shows the ideal scaling 12 CHAPTER 1 INTRODUCTION 15 1 1 System Requirements PLUTO can run on most platforms but some software prerequisites must be met depend
126. ixed number of time steps set by the macro FARGO_NSTEP_AVERAGE default is 10 e NO the velocity is prescribed analytically using the FARGO_SetVelocity function that can be implemented in your init c This must be the default if FARGO is used together with the shearing box module Initial and boundary conditions are assigned as usual by prescribing the total velocity and not the residual Likewise output files are always written using the total velocity and not the residual The order of reconstruction used during the linear transport step is set by the constant FARGO_ORDER which by default is 3 third order PPM The default value of the three switches FARGO_ORDER FARGO_NSTEP_AVERAGE and FARGO_AVERAGE_VELOCITY can be changed inside your definitions h see 23 EN FARGO MHD is typically used to model supersonic accretion disks and test problems can be found in the directory Test Problems HD Disk Planet configurations 2 4 and 6 as well as in Test Problems MHD FARGO Spherical Disk For more information see the test problem documentation at Doc test_problems html 10 2 2 A Note on Parallelization The FARGO MHD algorithm is fully parallelized in all coordinate directions with the requirement that the number of zones per processor in the orbital direction must be larger than the expected transport shift denoted with m CHAPTER 10 ADDITIONAL MODULES 89 With a large number of processors 2048 the resulting auto decomposition m
127. ject to some restrictions e The allowed modules are HD and MHD special relativistic counterparts are not yet implemented e Inthe case of the MHD module only cell centered magnetic fields are supported i e DIV_CLEANING Temporal integration can be performed only with RK3 split or unsplit e Only Cartesian coordinates are supported in any number of dimensions FD schemes are based on a global Lax Friedrichs flux splitting and the reconstruction step is per formed for robustness issues on the local characteristic fields computed by suitable projection of the positive and negative part of the flux onto the left conservative eigenvectors For this reason these schemes are more CPU intensive than traditional FV schemes approximately a factor 2 to 3 5 although can achieve the same accuracy with much fewer points Unlike the FV schemes currently present in PLUTO possessing an overall 2 order accuracy schemes provided by the conservative FD module are genuinely third or fifth order accurate The latter in particular have shown to outperform traditional second order TVD schemes in terms of reduced numerical dissipation and faster convergence rates for problem involving smooth flows Figure shows as a qualitative example a comparison between traditional FV methods such as Muscl Hancock or PPM and some FD methods on a problem involving circularly polarized Alfven waves see Test_Problems MHD CP_Alfven Although FD schemes can correc
128. l strategy to enforce the V B 0 constraint Possible values are e NO divergence constraint is not enforced Recommended for one dimensional problems or 2D config urations with purely azimuthal fields eo EIGHT WAVES magnetic fields retain a cell average representation and the eight wave formulation introduced by Powell Pow94 is used see 6 2 3 1 eo DIV CLEANING magnetic fields retain a cell average representation and the mixed hyperbolic parabolic diver gence cleaning technique of DKK 02 is used see 46 2 3 2 A new scalar variable the generalized Lagrange multiplier Y PSI_GLM is introduced e CONSTRAINED_TRANSPORT the magnetic field has a staggered representation and the constrained transport is used see 46 2 3 3 2 2 3 EOS Select the equation of state EOS e IDEAL use the perfect gas law with constant specific heat ratio This EOS is available for all physics module and it is described in e ISOTHERMAL use an isothermal equation of state 87 1 In this case the energy equation is not included eo PVTE_LAW allows the user to employ implement more complex equations of state by specifying the caloric EOS as e e p T This EoS works for the HD and MHD module and is described in e TAUB use the Taub Matthews equation of state only for relativistic modules see Please refer to Chapter 7 for a detailed description CHAPTER 2 PROBLEM HEADER FILE DEFINITIONS H 24 2 24 ENTROPY_SWITCH By enabling this swi
129. le extension listed above By default data files are written in the local working directory unless a different location has been specified using output dir in your pluto ini see There s no distinction between serial or parallel mode _Basename Variable Single record size rho Density Ni x No X N3 prs Pressure Nix Na X N3 vxl x velocity Ni x No x N3 vx2 x2 velocity Ni x No x N3 vx3 x3 velocity Nix Na x N3 bx1 x1 mag field Nix Na x N3 bx2 x2 mag field Nix Na x N3 bx3 x3 mag field Ni x Na x N3 bx1s x stag mag field Ni 1 x No x N3 bx2s x2 stag mag field Ny x N2 1 x N3 bx3s x3 stag mag field Ny x No x N3 1 trl first tracer Ni xX No x N3 Table 11 1 Base prefix for multiple data set The size is in units of 4 for the f1t format or 8 for the db1 format bytes Bitmap image format that employs lossless data compression 92 CHAPTER 11 OUTPUT AND VISUALIZATION 93 For each format it is possible to dump all or just some of the variables Additional user defined variables may be written as well The default setting is described separately for each output in the next subsections and may be changed if necessary see Each format has an independent output frequency and an associated log file i e dbl out flt out vtk out and so forth keeping track of the dump history Two additional files grid out and sysconf out contain grid and system related information respectively Important some visualization packages need the in
130. le x1 if GEOMETRY CARTESIAN return 1 0 sqrt xlxx1 x2xx2 x3x x3 elif GEOMETRY CYLINDRICAL double x2 double x3 return 1 0 sqrt xls xl x2 x2 elif GEOMETRY SPHERICAL return 1 0 x1 endif for the three coordinate systems e The acceleration vector can be employed when the BODY_FORCE flag is set to VECTOR and the three components of g are prescribed using the function BodyForceVector void BodyForceVector double xv double xl double xg double x2 double x3 where xv a pointer to a vector of primitive quantities e g v RHO v VX1 etc xg a three component array g IDIR g JDIR g KDIR specifying the gravity vector g components along the coordinate directions x1 x2 x3 local zone coordinates As an example let s consider again a point mass source located at the origin of coordinates Then one needs to define depending on the geometry CARTESIAN CYLINDRICAL or SPHERICAL void BodyForceVector double v double xl double x2 double x3 double xg double gs rs Hif GEOMETRY CARTESIAN rs sqrt lex x2xx2 Xx3 x3 elif GEOMETRY CYLINDRICAL rs sqrt xl xx1 x2 x2 elif GEOMETRY SPHERICAL rs x1 endif spherical radius in cart coords spherical radius in cyl coords Ta coords spherical radius in sph 1 0 rs rs spherical gravity if GEOMETRY CARTESIAN g IDIR g JDIR g KDIR
131. leration vector g 45 3 Note that the induction equation may equivalently be written as OB VxE 0 6 6 V xE 0 6 6 where v x B n J is the electric field and 77 is the resistivity tensor The right hand side of the system of Eqns is implemented in the RightHandSide function inside Src MHD rhs c employing a conservative discretization that closely follows the expression given in A 1 2Jand A 1 3 for Cartesian polar and spherical geometries The sets of conservative and primitive variables U and V are given by T T U p m E B y V p V P B The maps U V and its inverse are provided by the functions PrimToCons and ConsToPrim The primitive form of the equations is P ey VOR po 0 O 1 1 Tio Vo Bx VxB Vp ot p p 2 B V v B V v v V B v V B O ay TU Vet pavo 0 VO g 6 7 where the V B on the right hand side of the third equation is kept for reasons of convenience although zero at the continuous level LA factor of 1 47 has been absorbed in the definition of magnetic field CHAPTER 6 BASIC PHYSICS MODULES 59 6 2 2 Assigning Magnetic Field Components Magnetic field components are initialized in your Init function just like any other flow quantity Depending on the value of ASSIGN_VECTOR_POTENTIAL in your definitions h two alternative initializa tions are possible 1 By setting ASSIGN_VECTOR_POTENTIAL to NO default you can assign the component of magnet
132. les Starting with Gnuplot 4 2 raw binary files are also supported In this case grid information being stored in separate files must be supplied explicitly through appropriate keywords making the syntax a little awkward To ease up this task one can take advantage of the scripts provided with the code distribution in Tools Gnuplot For this we recommend to define the GNUPLOT_LIB environment variable in your shell which will be appended to the loadpath of Gnuplot gt export GNUPLOT_LIB PLUTO_DIR Tools Gnuplot use setenv for tcsh users You can also define the loadpath directly from Gnuplot gnuplot gt set loadpath lt pluto full path gt Tools Gnuplot A typical gnuplot session can be started as follows gnuplot gt load grid gp read and store grid information gnuplot gt dsize 8 load macros gp gnuplot gt load pm3D_setting gp set the display canvas for pm3d plot style The first line invokes the grid gp script which is used to read grid information the second script sets the data file size 8 or 4 for double or single precision while the last one initializes a default environment for viewing binary data files using the pm3d style of Gnuplot For additional documentation and examples please refer to Doc gnuplot html CHAPTER 11 OUTPUT AND VISUALIZATION 100 11 3 2 Visualization with IDL IDL Interactive Data Language is a vectorized programming language commonly used in the astro nomical community for inter
133. m the archive gt gunzip plu uto xx tar gz gt tar xvf pluto xx tar this will create the folder PLUTO in your home directory At this point we advise to set the environment variable PLUTO_DIR to point to your code directory Depending on your shell e g tcsh or bash use either one of gt export PLUTO_DIR home user PLUTO If you re using the bash shell gt setenv PLUTO_DIR home user PLUTO If you re using the tcsh shell 0 2 Running a simple shock tube problem PLUTO can be quickly configured to run one of the several test problems provided with the distribution Assuming that your system satisfies all the requirements described in the next chapter i e C compiler Python etc you can quickly setup PLUTO in the following way 1 Change directory to any of the test problems under PLUTO Test_Problems e g gt cd SPLUTO DIR Test Problems HD Sod 2 Copy the header and initialization files from a configuration of our choice e g 01 PLUTO Test_Problems HD Sod gt cp definitions_01 h definitions h PLUTO Test_Problems HD Sod gt cp pluto_01 ini pluto ini 3 Run the Python script using PLUTO Test_Problems HD Sod gt python PLUTO_DIR setup py and select Setup problem from the main menu see Fig 1 2 You can choose by pressing Enter or modify the default setting using the arrow keys 4 Once you return to the main menu select Change makefile choose a suitable makefile e g Linux
134. maining files pluto ini and init c should be appropriately edited by the user Templates for all four files can be found in the Src Templates directory Several examples are located in the test directories under the Test_Problems directory In order to run the Python script anywhere from your hard disk we recommend to set the shell variable PLUTO_DIR to point to your PLUTO distribution Depending on your environment shell use either one of gt setenv PLUTO_DIR home user PLUTO if you re using tcsh shell gt export PLUTO_DIR home user PLUTO if you re using bash shell The setup py script can now be invoked with MyWorkDir gt python SPLUTO_DIR setup py options Command line options are listed in Table or can be briefly described by invoking setup py with help By default the Python script uses the ncurses library for enhanced terminal control However this option may be turned off by invoking the setup script with the no curses switch You should ther see the menu shown in Fig Additional menus depending on the physics module will display later Python will first create an architecture dependent file named sysconf out containing system related information this file does not have any specific purpose but may be helpful for the user Whenever an internet connection is available Python will also notify if new versions of the code are available CHAPTER 1 INTRODUCTION 15 600 xterm gt gt Python setup Aug 20
135. mperature T p p and the z component of vorticity w 0 v y O vz Then one has to set uservar 2 T vortz in your pluto ini under the Static Grid Output block This informs PLUTO that 2 additional variables named T and vortz have to be saved They are computed at each output by editing the function ComputeUserVar void ComputeUserVar const Data xd Grid int 1 17 double T VOrtz double xxxp x xrho VX x VY double xdx xdy GetUserVar T GetUserVar vortz d gt Vc d gt Vc d gt Vc dqd gt Vc RHO PRS VX1 VX2 x pointer shortcut to density pointer shortcut to pressure polnter shortcut to x velocity pointer shortcut to y velocity l grid IDIR dx 4 shortcut to dx Grid UDIR dx shortcut toa de 4 7 DOM_LOOP k 3 1 T k j i vortz k j i PLUTO automatically allocates static memory area for the new variables T and vortz when calling the GetUserVar function The DOM_LOOP macro performs a loop on the whole computational domain boundary excluded in order to compute T k 3 i and vortz k 3 i Once PLUTO runs these two variables will automatically be written in all selected formats except for the ppm and png formats by default Beware that if the number of uservar is reset to zero but the previous function is still executed a segmentation fault error will occur since user defined variables such
136. mprising of atomic molecu lar and ionized hydrogen Here the total number density n is set to be 10 cm 3 while the fractions of dif ferent hydrogen species are X_HI 0 835 X H2 0 0823 and X_HII 0 0004 losses due to rotational vibrational cooling dissociation and gas grain processes Their variation with temperature for a particular set ofm hydrogen fractions is shown in fig 9 1 Depending on the require ment the user can add more components to the cooling function for e g cooling due to fixed fractions of standard molecules like CO OH H20 etc or contributions from collisional excitation of lines as indi cated in the SNEg module CHAPTER 9 OPTICALLY THIN COOLING 85 9 5 Multi lon Non Equilibrium Cooling MINEq This module computes the dynamical evolution and ionization state of the plasma using the multi ion model of including with 28 ion species namely HI Hel Hell and the first five ionization stages of C N O Ne and S For each ion PLUTO introduces an additional variable the fractional abundance of the ion with respect to the element it belongs ida lion Melem The names of the additional variables for the corresponding species are X HI Hel HeI1 X_CI X_CII CIT ACT XOV NT NEL XNIITXNIV XN XOI XOTI XOIII XOTV XOV NET X_NeII X_NeITI X_NelV X_NeV X SI X_SII X_SITI X_SIV X_SV Ionized hydrogen is simply 1 v X_HT You can assign the fraction of any ion specie by setting in the usual fash
137. must be CFL lt 0 7 in 2D and CFL lt 0 35 in 3D e CT_EMF_AVERAGE UCTO or CT_EMF_AVERAGE UCT_CONTACT employs the face to edge inte gration procedures proposed by GS05 where electromotive force derivatives are averaged from neighbor zones UCTO0 or selected according to the sign of the contact mode UCT_CONTACT The former has reduced dissipation and is preferably used with linear interpolants and RK integrators while the latter shows better dissipation properties All algorithms with the exception of the arithmetic averaging reduce to the corresponding one dimen sional scheme for grid aligned flows However in our experience UCT_HLL and UCT_CONTACT show the best dissipation and stability properties The CT formulation works with any of the Riemann solvers Assigning Boundary Conditions Within the CT framework user defined boundary conditions b c must be assigned on the staggered components as well This is done in your UserDefBoundary function using the d gt Vs nv k j i array where nv gives the staggered component BX1s BX2s or BX3s Note In PLUTO we follow the convention that the cell center owns its right interface e g 1 means i 5 Thus baitt jk IVs BX18 k j i dno ijt k d gt Vs BX2s k 131 i bes ijk i IVs BX3s k j i CHAPTER 6 BASIC PHYSICS MODULES 63 Beware that the three staggered components have different spatial locations and the BOX_LOOP
138. n A r2A T with n 9 5 HMu when the cooling heating function A T is not known analytically but rather is available as a table sampled at discrete not necessarily equidistant points i e A A Z In order to use this module the user must provide a two column ascii files in the working directory named cooltable dat of the form T J Lambda 3 with the temperature expressed in Kelvin and the cooling heating function A in ergs cm s An example of such fil can be found in Src Cooling Tab cooltable dat As usual the dimensionalization is done automatically by the cooling module once UNIT_DENSITY UNIT_LENGTH and UNIT_VELOCITY have been defined in Init Alternatively the TABULATED cooling module can be used to provide a user defined cooling func tion A A V 9 6 where V is a vector primitive variables The explicit dependence of A can be given by i copying Src Cooling Tab radiat c into your local working directory and ii make the appropriate changes Generated with Cloudy 90 01 for an optically thin plasma and solar abundances thanks to T Plewa CHAPTER 9 OPTICALLY THIN COOLING 82 9 3 Simplified Non Equilibrium Cooling SNEq This module is implemented in the Src Cooling SNEq directory and introduces a new variable with index X_HT used to label the fraction of neutrals z pz n LHI 9 7 nH You can assign the fraction of neutrals by setting in the usual fashion v X_HI 0 2 x for exam
139. n unstable integration whereas values close to 1 can decrease STS s efficiency By default y 0 01 a value which in many cases retains stability whereas giving substantial gain see Fig To change the default value of y STS_NU redefine it in the user defined symbolic constant section of definitions h see Figure 8 1 Length of a super step in units of the ex plicit one AT Atpar as function of the number of sub steps N using different values of y 107 green plus sign v 107 red asterisk default v 107 pur ple square The upper dotted lines gives the y 0 limit AT N whereas the lower one represents the explicit limit AT N If AT Atpar 100 for ex ample explicit integration would require 100 steps while super time stepping only 21 for v 1072 or 11 for y 1073 steps Stability analysis for the constant coefficient diffusion equation Bec92 indicates that the value of Cp parabolic Courant number should be lt 1 Naim Naim is the number of spatial dimensions and it may be used to adjust the size of the spectral radius for strongly nonlinear problems A reduction of Cp CHAPTER 8 DISSIPATIVE EFFECTS 78 will results in increased stability at the cost of more substeps N The default value is C 0 8 Naim but it may be changed in your pluto ini through CFL_par see 44 3 Since STS treats parabolic equations in an operator split formalism it may be advisable for highly nonline
140. nd VTK files vtk can be visualized using the pyPLUTO code This tool is included in the current code distribution in the directory Tools pyPLUTO and provides python modules Python version gt 2 7 is recommended to load visualize and analyse data Additionally for the purpose of a quick check Graphic User Interface GUI is provided requires Python Tkinter Details of the Installation and Getting Started can be found in the Doc pyPLUTO html On successful installation the user can load data in the following manner gt 1python pylab In 1 import pyPLUTO as pp for loading data 0010 dbl In 2 D pp pload 10 w_dir lt path to data dir gt for loading data 0010 f1t In 3 D pp pload 10 w_dir lt path to data dir gt datatype float Here D is a pload object that has all the information regarding the variable names and their values which are stored as arrays It also has the respective grid and time information For example D x1 is the numpy x array D rho is the numpy density array D vx1 is the numpy vxl array and so on These numpy arrays can be easily visualized using matplotlib python s plotting library The pyPLUTO s also provides two classes Image and Tools They have some frequently used functions for analysis and data plotting Details about these classes along with their usage can be found in HTML document referred above In order to use the GUI version for visualizing the data append
141. ng c as explained in section 3 3 of aa The allowed values are 1 2 or 3 1 use the cell centered value wf wio This choice is slightly more diffusive but has found to work well for flows containing strong discontinuities 2 No reference state is introduced and the interpolated states at the base time level are used w wi This is found to be a good choice in presence of smooth flow and equlibrium configurations 3 reference states are constructed as in the original PPM al gorithm CW84 Col90 to minimize the size of the term susceptible to characteristic limiting see Eq 29 and 30 of IMZT 12 The default value is 3 except for PARABOLIC WENO reconstruc tion in characteristic variable for which 2 is used Sets the maximum shock strength above which fluid variables inside a given computational zone can be safely updated us ing the entropy equation see and the source file Sr c flag_shock c It has effect only when ENTROPY_SWITCH has been enabled A lower value will trigger the flattening proce dure in more zones Default is 0 05 Sets the minimum shock strength above which the MULTID shock flattening algorithm flags a zone to be inside a shock see Src flag shock c It has effect only when SHOCK_FLATTENING is set to MULTID A lower value will trigger the flattening pro cedure in more zones Default is 5 0 Set this to YES if the FARGO orbital velocity w should be computed by averaging the azimuthal velocity
142. ng field beginning with a FP Libraries can now be built under Chombo 3 2 lib with gt make lib Do not try make all since it won t work after the Chombo patch file has been unpacked Alternative configurations can be made from the default one by specifying the configuration vari ables explicitly on the make command line For example gt make DIM 3 MPI TRUE lib will build the parallel version of the 3D library Additional information may be found in the Chom bo README file and by consulting the library documentation CHAPTER 12 ADAPTIVE MESH REFINEMENT AMR 109 12 2 Configuring and running PLUTO Chombo In order to configure PLUTO with Chombo you must start the Python script with the with chombo option Python assumes that Chombo libraries have been built under PLUTO Lib Chombo 3 2 work gt python PLUTO_DIR setup py with chombo This will use the default library configuration 2D serial in the example above To use a configuration different from the default one enter the make configuration variables em ployed when building the library e g work gt python SPLUTO_DIR setup py with chombo MPI TRUE do not use spaces in MPI TRUE Note that the number of dimensions DIM is specified during the Python script and should NOT be given as a command line argument The setup proceeds normally as in the static grid case by choosing Setup problem from the Python script to change configure your test problem The make
143. ngth in cm Default value is 1 astronomical unit Sets the unit velocity in cm sec Default value is 1 Km sec Enable writing of time information to vtk output files Notice that this information is useful only when reading data files with VisIt and may give problems with other visualisation softwares 411 1 3 Default value is NO Enable writing of vector fields velocity and magnetic field during VTK output 11 1 3 Default value is NO all variables are written with the scalar attribute 123 Bibliography AAG96 AAZN97 Bal86 BCCD08 Bec92 BS99 BS03 BTHO8 CM77 Col85 Col90 CT09 CW84 DBL03 DKK 02 DZBLO07 GP98 GS05 HEOC87 HM79 1596 Vasilios Alexiades Genevieve Amiez and Pierre Alain Gremaud Super time stepping acceleration of explicit schemes for parabolic problems Communications in Numerical Methods in Engineering 12 1996 no 1 31 42 T Abel P Anninos Y Zhang and M L Norman Modeling primordial gas in numerical cosmology New Astronomy2 1997 181 207 S A Balbus Magnetized thermal conduction fronts ApJ304 1986 787 798 R Borges M Carmona B Costa and W S Don An improved weighted essentially non oscillatory scheme for hyperbolic conservation laws Journal of Computational Physics 227 2008 3191 3211 J M Beckers Analytical linear numerical stability conditions for an anisotropic three dimensional advection diffusion equ
144. ode may result in sub domains that violate this condition and an error message is issued To avoid this problem you can specify the parallel decomposition manually using the dec nl n2 n3 command line argument 1 4 1 and ensure that not too many processors are used along the direction As an example suppose you wish to use 4096 processors but only 8 along the orbital direction x2 You may specify the domain decomposition by giving say 32 8 and 16 in the three directions with MNpilrun ap 4096 pluto dec 32 8 16 CHAPTER 10 ADDITIONAL MODULES 90 10 3 High order Finite Difference Schemes An alternative to the Finite Volume FV methodology presented in the previous chapters and to the reconstruction algorithms described in Chapter 2 is the employment of conservative high order Finite Difference FD schemes 3 and 5 order accurate in space interpolation can be used in PLUTO invoking setup py with the following extension MyWorkDir gt python SPLUTO_DIR setup py with fd The available options in RECONSTRUCTION will now be e LIMO3_FD third order reconstruction of CTO9 e WENO3_FD an improved version of the classical third order WENO scheme of JS96 based on new weight functions designed to improve accuracy near critical points YC09 e WENOZ_FD improved WENO5 scheme proposed by e MP5_FD the monotonicity preserving scheme of based on a fifth order interface value The use of high order FD schemes is sub
145. odic domain with 120 x 20 zones The different curves plot the max imum value of Bz as a function of time and thus give a measure of the intrinsic nu WENO3_FD GLM merical dissipation Selected finite volume WENOZ_FD GLM schemes employing constrained transport MPO_FD GLM CT are MUSCL HANCOCK MH CT Runge Kutta 2 RK2 CT and PPM CT Fi nite difference schemes employ the GLM formultation and are respectively given by WENO3 WENOZ and MP5 10 3 2 LimO3 amp MP5 As an alternative to the previously described WENO schemes LimO3 and MP5 interpolations are also available The former is a new and efficient third order limiter function proposed by CTO9 Utilizing a three point stencil to achieve piecewise parabolic reconstruction for smooth data LimO3 preserves its accuracy at local extrema avoiding the well known clipping of classical second order TVD limiters Note that this reconstruction is also available in the finite volume version of the code PLUTO s MP5 originates from the monotonicity preserving MP schemes of SH97 which achieve high order interface reconstruction by first providing an accurate polynomial interpolation and then by limiting the resulting value in order to preserve monotonicity near discontinuities and accuracy in smooth regions The MP algorithm is better sought on stencils with five or more points in order to distinguish between local extrema and a genuine O 1 discontinuities For an inter scheme compari
146. ons define PV_TEMPERATURE_TABLE YES define TV_ENERGY_TABLE NO tell PLUTO to replace the thermal EOS with a temperature table T T p pj while still using the CHAPTER 7 EQUATION OF STATE 71 analytical approach i e direct function evaluation or root finder for the caloric EOS The tables T p pj and pe T p are initialized at runtime and used throughout the integration The number of points needed to construct such tables is fixed by the constants PV_TEMPERATURE_TABLE_NX PV_TEMPERATURE_TABLE NY for the first table and TV_ENERGY_TABLE_NX TV_ENERGY_TABLE NY for the second one To avoid the occurrence of spurious waves in the solution of the Riemann problem a monotone cubic spline is always used in the temperature grid see Appendix C of VMBM15 7 4 The TAUB Equation of state The Taub Matthews IM equation of state is available to describe a relativistic perfect gas for which the adiabatic exponent is a function of the temperature The actual expression for the Synge gas Syn57 is rather complex and PLUTO employs a quadratic approximation to the theoretical relativistic perfect gas EOS T gt 5 3 in the low temperature limit and T gt 4 3 in the high temperature limit see MM07 1 a 2 1 42 j 7 13 where h is the specific enthalpy related to the internal energy and pressure through h 1 e p 8 Dissipative Effects In this chapter we give an overview of the code capabilities for treating dissipat
147. osition does not explicitly depend on temperature and density e g mu MeanMolecularWeight v where v is an array of primitive variables of which only ion fractions need to be defined den sity and pressure are ignored The mean molecular weight is implemented in the source file mean_mol_weight c for a variety of different cooling reaction network Ch 9 and also when no cooling is used For the SNEg cooling module for instance the mean molecular weight is computed as _ Ap Anefue Azfz 2 Jar Jue 21z 64 where each metal contributes for one electron In the previous equation Ay CONST_AH and Ape CONST_He are the atomic mass numbers of hydrogen and helium whereas fur fue and fz are the number fractions of neutral hydrogen helium and metals with respect to hydrogen Nutr Ny Nye Y An O Nz Z Ay xe Mr A Tie Ny Ay X fai where X Muny Ay p Y myunyeAnue p while Z 1 X Y Notice that while fye and fz are fixed fyz is a time dependent quantity that evolves with flow variables CHAPTER 5 INITIAL AND BOUNDARY CONDITIONS INIT C 46 Without any cooling network the mean molecular weight is computed from the given mass fractions assuming a fully ionized gas fH 0 qa An Aefue Azfz 2 fue fz Az 2 where Az Nz 2 is the number of electrons due to metals The value of X and Y can be specified through the user defined constants HMASS_FRAC and He_MASS_FRAC which by de
148. owell Approximate Riemann solver for magnetohydrodynamics that works in more than one dimen sion Tech report March 1994 K G Powell P L Roe T J Linde T I Gombosi and D L De Zeeuw A Solution Adaptive Upwind Scheme for Ideal Magnetohydrodynamics Journal of Computational Physics 154 1999 284 309 P Rossi G Bodo S Massaglia and A Ferrari Evolution of Kelvin Helmholtz instabilities in radiative jets II Shock structure and entrainment properties A amp A321 1997 672 684 P L Roe Approximate Riemann Solvers Parameter Vectors and Difference Schemes Journal of Computa tional Physics 43 1981 357 372 J Saltzman An Unsplit 3D Upwind Method for Hyperbolic Conservation Laws Journal of Computational Physics 115 1994 153 168 J M Stone and T A Gardiner Implementation of the Shearing Box Approximation in Athena ApJS189 2010 142 155 A Suresh and H T Huynh Accurate Monotonicity Preserving Schemes with Runge Kutta Time Stepping Journal of Computational Physics 136 1997 83 99 125 SO89 Spi62 Str68 Syn57 TMMO08 Tor97 van79 VMBM15 WAMMO07 YCO9 C W Shu and S Osher Efficient Implementation of Essentially Non oscillatory Shock Capturing Schemes II Journal of Computational Physics 83 1989 32 78 L Spitzer Physics of Fully Ionized Gases 1962 G Strang On the Construction and Comparison of Difference Schemes SLAM Journal on Numerical An
149. oximate Riemann solver for ideal magnetohydrodynamics Journal of Computational Physics 208 2005 315 344 J M Mart and E M ller Extension of the Piecewise Parabolic Method to One Dimensional Relativistic Hydrodynamics Journal of Computational Physics 123 1996 1 14 A Mignone and J C McKinney Equation of state in relativistic magnetohydrodynamics variable versus constant adiabatic index MNRAS378 2007 1118 1130 A Mignone T Plewa and G Bodo The Piecewise Parabolic Method for Multidimensional Relativistic Fluid Dynamics ApJS160 2005 199 219 A Mignone and P Tzeferacos A second order unsplit Godunov scheme for cell centered MHD The CTU GLM scheme Journal of Computational Physics 229 2010 2117 2138 A Mignone P Tzeferacos and G Bodo High order conservative finite difference GLM MHD schemes for cell centered MHD Journal of Computational Physics 229 2010 5896 5920 A Mignone M Ugliano and G Bodo A five wave Harten Lax van Leer Riemann solver for relativistic magnetohydrodynamics MNRAS393 2009 1141 1156 A Mignone C Zanni P Tzeferacos B van Straalen P Colella and G Bodo The PLUTO Code for Adaptive Mesh Computations in Astrophysical Fluid Dynamics ApJS198 2012 7 S Orlando F Bocchino F Reale G Peres and P Pagano The Importance of Magnetic Field Oriented Thermal Conduction in the Interaction of SNR Shocks with Interstellar Clouds ApJ678 2008 274 286 K G P
150. pecially if angular coordinates are present and or steady state solutions are sought 2 1 4 BODY FORCE Include a body force in the momentum and energy equations Possible values are e POTENTIAL body force is derived from a scalar potential pa pV O e VECTOR body force is expressed as a three component vector pa pg e VECTOR POTENTIAL body force is prescribed using both pa p V9 g More details can be found in 2 1 5 COOLING Optically thin thermal losses can be included by appropriately setting this flag to one of the following e POWER_LAW radiative losses are proportional to p T 49 1 e TABULATED radiative losses are computed as n A T where A T is a user supplied tabulated function of temperature see Alternatively this module can be used to provide user defined cooling functions e SNEq simplified non equilibrium cooling function for atomic hydrogen See 9 3 for more details e H2_COOL optically thin cooling function for molecular and atomic hydrogen See e MINEq multi ion non equilibrium cooling model It evolves the standard equations augmented with a chemical network of 29 ions see q9 5and the work by TMMO08 CHAPTER 2 PROBLEM HEADER FILE DEFINITIONS H 19 2 1 6 RECONSTRUCTION Sets the spatial order of integration In the standard finite volume version of the code the following options are available e FLAT first order reconstruction The stencil is 1 point e LINEAR piecewise TV
151. pherical coordinates r 0 6 the RMHD equations 6 11 are discretized using the following diver gence form OD OE V Dv 0 OM 9 OD Mava MgV gt V w b u v bb Sy P9r 7 Bo Bo By By S w Bo ge S Buo Pa Ome 9 1 Op MgUr Cot IMgUE AY V w b ugv bob D9 pgo i Bo B Bg B mn B SLE e B E v jw 7 cot 0 E v Wwe OMe m 9 1 Op E OE V w b JU U byb us rsin Od Pgo DA Dv g Ot OB 1 O sin eee B 1 DES 0 Ot r sin 0 00 rsin0 do 7 3 Bo 1 O 7 1 rEg _ 9 Ot rsin0 dp r Or 7 OBS 1 0 rEp _ 10 Ot r or r 00 A 9 Note that curvature terms are present in the radial and meridional components while the azimuthal component is discretized in angular momentum conserving form The corresponding divergence oper ators are 1 OTE 1 O sin OF 1 OF 7 ua ara 24 Y X r2 Or rsinOd 00 rsin 09 A 10 oo d O r3 F yl O sin 9 Fa 1 OF r Or r sin 0 00 rsin Od 118 B Predefined Constants and Macros B 1 Predefined Physical Constants PLUTO has several predefined physical and astronomical constants in c g s units which may be used anywhere in the code see macro h define CONST_AH xx lt Atomic weight of Hydrogen x define CONST_AHe i xx lt Atomic weight of Helium x define CONST_AZ xx lt Mean atomic weight of heavy elements x define CONST_amu 66053886e 24 fea Atomic mass unit ul defin
152. ple lt in your Init function The fraction of neutrals is treated as a passive scalar during the hydro step while it is governed by the following ODE during the cooling step O HT o HL S ne er ci fa cr 9 8 together with the energy equation t A NenNH gt Jk Wi r 9 9 where the summation over k accounts for 16 different line emissions coming from some of the most common elements k Ly a H a Hel 584 623 CI 9850 9823 CII 1561 CII 2325A NI 5200 A NII 6584 6548 A OI 63u OI 6300 6363 A OI 3727 MglI 2800 Sill 351 SII 6717 6727 Fell 2511 Fell 1 61 The coefficient 7 in has dimensions of erg sec cm and is computed from hi 2r hvk Jk fka VkBMeme 1 nelq21 A21 where k is the index of a particular transition fk nk ny is the abundance for that particular species Here 86 10 Qi hvy 8 6 107 Qa q12 IT 7 kT q21 IT as where 212 221 is the collision strength and is tabulated In Eq w r represents the thermal energy lost by ionization and recombination T Wir Ci X 13 6 x 1 6 107 fn cr x 0 67 x 1 6107 1 fn a where c and c are the hydrogen ionization and recombination rate coefficients 26 104 Po 1 08 10 VT i VT 7 13 6 VT CHAPTER 9 OPTICALLY THIN COOLING 83 Table 9 1 Summary of the chemistry reaction set T is the temperature in Kelvin Tey is the temperature in electron volts and T
153. potential in time and save it to disk Note that ASSIGN_VECTOR_POTENTIAL must be enabled 3 Makefile Selection makefile The makefile contains instructions to compile and link C source code files and produce the executable pluto The Python script creates a new makefile every time you choose Change makefile from the menu otherwise it automatically updates the existing one after you have finished the problem setup If you choose to create a new makefile Python will ask you to select an appropriate defs file con taining architecture dependent flags from the Config directory The template Config Template defs can be used to create a new configuration from scratch The simplest example is a definition file for a single processor without any additional library In this case it suffices to set GE cc CFLAGS E 0 LDFLAGS 1m PARALLEL FALSE TRUE FALSE enable disable parallel mode USE _HDF5 FALSE TRUE FALSE enable disable support for HDF5 library USE_PNG FALSE TRUE FLASE enable disable support ofr PNG library where CC is the name of your C compiler cc gcc mpicc etc CFLAGS are command line options such as optimization search path etc and LDFLAGS contains options to be passed to the linker The variables PARALLEL USE_HDF5 and USE_PNG can be set to either TRUE or FALSE to enable or disable parallel mode support for HDF5 library and support for PNG library respectively in the static grid version of PLUTO When set
154. ptions are e EULER 1 explicit Euler algorithm is used to evolve from U to U U SUFAN e RK2 RK3 2 or 3 order TVD Runge Kutta is used to advance the solution to the next time level RK2 RK3 U U At L U Ae UH L U U AL UH 1 U 2U 2AP L When DIMENSIONAL SPLITTING YES the operator in Eq is one dimensional Setting DIMENSIONAL_SPLITTING to NO makes the scheme dimensionally unsplit and the right hand side include contributions from all directions simultaneously Unsplit implementation of the Runge Kutta algorithms usually requires a more restrictive CFL condition see Table 2 1 e CHARACTERISTIC_TRACING HANCOCK they evolve U according to Unt UP Arey a where V is computed by suitable Taylor expansion Although the final step is in diver gence form these methods require the primitive formulation of the equations not yet available CHAPTER 2 PROBLEM HEADER FILE DEFINITIONS H 20 for all modules They are 2 order accurate in space and time and less dissipative than the previous multi step algorithms HANCOCK should be combined with linear reconstruction while CHARACTERISTIC_TRACING which does a more sophisticated characteristic limiting can be com bined with all reconstruction algorithms The original PPM scheme of is available for the HD MHD and RHD modules by selecting TIME_STEPPING to CHARACTERISTIC_TRACING together with RECONSTRUCTION to PARABOLIC and
155. r mation are supplied using the INCLUDE_DIRS and LDFLAGS variables respectively USE_HDF5 TRUE HEPES EE HEHEHE HE HHH EH HE HH HDF5 library options HARTA AAA ifeg strip USE_HDF5 TRUE HDF5_LIB usr local lib HDF5 1 xx INCLUDE_DIRS I HDF5_LIB include LDFLAGS LS HDF5_LIB lib 1hdf5 1z endif Note PLUTO uses the HDF5 1 6 API although it may be linked with HDF5 1 8 x without any problem since the H5_USE_16_API macro defined in hdf5_io c forces the library to use HDE5 1 6 macro definitions 1 Contrary to a blocking call which will not return until the I O request is completed a non blocking call initiates an I O operation but does not wait for it completion CHAPTER 3 MAKEFILE SELECTION MAKEFILE 31 3 3 PNG Library Support Whenever the portable network graphics PNG library is installed on your system you may enable support for 2D output using PNG color images To do so simply set to TRUE the corresponding USE_PNG variable inside your defs file and add the linker option to the LDFLAGS variable USB PNG TRUE HARE EEE HEHEHE HE EHH HH HE EF PNG library options PERRA HEHEHE HEE EHH HH HE EHH HH HE HEF ifeq USE_PNG TRUE LDFLAGS lpng endif 3 4 Including Additional Files local_make Additional e g user defined files may be added to the standard object list created by Python in your makefile To this end create a new file named local_make like OBJ myfile o HEA
156. r along the temperature axis Chombo LevelPluto step alpha lt 0 crash has been partially solved Tabulated cooling takes into account the mean molecular weight 1 Introduction PLUTO is a finite volume finite difference shock capturing code designed to integrate a system of conservation laws OU Ot where U represents a set of conservative quantities T U is the flux tensor and S U defines the source terms MBM 07 MZT 12 An equivalent set of primitive variables V is more conveniently used for assigning initial and boundary conditions The explicit form of U V T U and S U depends on the particular physics module selected V T U S U 1 1 e HD Newtonian classical hydrodynamics e MHD ideal resistive magnetohydrodynamics e RHD special relativistic hydrodynamics e RMHD special ideal relativistic magnetohydrodynamics PLUTO adopts a structured mesh approach for the solution of the system of conservation laws 1 1 Flow quantities are discretized on a logically rectangular computational grid enclosed by a boundary and augmented with guard cells or ghost points in order to implement boundary conditions on a given computational stencil Computations are done using double precision arithmetic The grid can be either static or dynamically adaptive as the flow evolves In the static grid version PLUTO comes as a stand alone package entirely written in the C programming language see MBM 07 for a comprehens
157. rator followed by the equal sign turns the FLAG_INTERNAL_BOUNDARY bit on in the 3D array d gt flag This is used by the code to reset the right hand side of the conservative equations in the selected zones to zero These computational cells are thus not evolved in time by PLUTO Note The box structure should not be used here and staggered magnetic field variables should not be altered 5 3 Body Forces Body forces are introduced by enabling the BODY_FORCE flag in your definitions h The force is computed in terms of the acceleration vector a a V0 9 5 5 where is the scalar potential and g 91 92 93 is a three component acceleration vector e The scalar potential can be employed when the BODY_FORCE flag is set to POTENTIAL in your definitions h In this case g 0 and the function BodyForcePotential should be used to prescribe the analytical form of 21 2 3 double BodyForcePotential double x1 double x2 double x3 where x1 x2 x3 are the local zone coordinates and the return value of the function gives the potential In this way PLUTO employs a conservative discretization that conserves total en ergy gravitational energy see Eq and Eq 6 4 The gravitational potential however must not change in time As an example a spherically symmetric point mass potential 1 r can be defined using CHAPTER 5 INITIAL AND BOUNDARY CONDITIONS INIT C 54 double BodyForcePotential doub
158. rature and density via the caloric equation of state Tor97 e e T p 7 4 The thermal and caloric equations of state given by Eq and constitutes the basis for the consideration discussed in the next sections 7 1 The ISOTHERMAL Equation of State In an isothermal gas the temperature is constant and the pressure is readily obtained as p pce 7 5 where Ciso the isothermal sound speed can be either a constant value or a spatially varying quantity This EOS is available only in the HD and MHD modules No energy equation is present and the labels ENG and PRS are undefined The value of ciso can be set using the global variable g isoSoundSpeed in your init c e g g_isoSoundSpeed 2 0 sets the sound speed to be 2 x If not set the default value is cis 1 In order to have a space dependent isothermal speed of sound one has to copy the source file Sr c HD eos c to your local working directory and make the appropriate modification 7 2 The IDEAL Equation of State For a calorically ideal gas the ratio of specific heats I is constant and the internal energy can be written pe PEA 7 6 CHAPTER 7 EQUATION OF STATE 68 The value of I is stored in the global variable g gamma and can be modified in your Init function default value 5 3 For a relativistic flow the constant EOS is more conveniently expressed through the specific en thalpy T ph p 5 3P 7 7 The ideal EOS is compatible with all physics mod
159. re e NO viscous terms are not included e EXPLICIT viscosity is treated explicitly e SUPER_TIME_STEPPING viscosity is treated using super time stepping See 8 1 for details CHAPTER 2 PROBLEM HEADER FILE DEFINITIONS H 26 2 3 User defined Constants In addition to the options described so far the user can insert an arbitrary number of user defined symbolic constants macros in the header file definitions h This should be done manually in the section delimited by x Beg user defined constants do not change this line x End user defined constants do not change this line of this file Only lines beginning with define should appear in this section as they will not be changed by the python script The value of a symbolic constant can be either a number or another symbolic constant previously defined by the code e g YES or NO and cannot be changed at runtime User defined symbolic constants are useful in the following circumstances 1 Conditional compilation useful when your initial configuration contains computationally expen sive code blocks that should be compiled separately As an example define in your definitions h the symbolic constant name as SETUP_VERSION and give it the value of 0 or 1 x Beg user defined constants do not change this line define SETUP_VERSION 1 x End user defined constants do not change this line x This symbolic macro name can then be used inside
160. rectly accounted for when RECONSTRUCTION is set to either LINEAR PARABOLIC or WENO3 Example 1 A simple uniform grid extending from zz 0 0 to zr 10 0 with 128 zones can be specified using xXi grid I Q 128 u 10 0 Example 2 consider a one dimensional physical domain extending from 0 0 to 10 0 with a total of 18 zones but a finer grid is required for 0 lt x lt 3 Then one might specify Xi grid 2 00 TZ gt 0 6 5 10 0 which generates a uniform grid with 12 zones for 0 lt x lt 3 and a stretched grid with 6 zones for 3 lt x lt 10 see Fig 4 1 12 zones u 6 zones s Figure 4 1 Example of one dimensional grid with a uni form left and stretched seg ment right in red covering the 18 zones interval 0 10 When the computational grid is generated each processor owns a domain portion defined by the global integer variables IBEG lt i lt IEND JBEG lt j lt JEND and KBEG lt k lt KEND Ghost cells are added outside the local computational domain to complete the stencil at the boundaries see Fig The global variables NX1 NX2 and NX3 define the total number of points boundaries excluded such that IEND IBEG 1 NX1 JEND JBEG 1 NX2 KEND KBEG 1 NX3 The total number of zones for a given processor boundaries included is given by the global variables NX1_TOT NX2_TOT and NX3_TOT see Fig However the cell aspect ratio can be different from unity that is rectangular cells
161. revious versions Although this has not been found to be a problem for VisIt many filters in ParaView such as streamlines may require to apply a Cell Data to Point Data filter 12 Adaptive Mesh Refinement AMR PLUTO provides adaptive mesh refinement AMR functionality in 1 2 and 3 dimesions through the Chombo library Chombo provides a distributed infrastructure for parallel calculations over block structured adaptively refined grids PLUTO Chombo is compatible with any of the available physics modules i e HD MHD RHD RMHD and grid refinement is supported in all coordinate systems Moreover grid zones are no longer constrained to be equilateral but sides can have different lengths Magnetic fields are evolved using cell centered schemes i e either Powell s EIGHT_WAVES or DIV_CLEANING Constrained transport is not yet available I O is provided by the Hierarchical Data Format HDF5 library designed to store and organize large amounts of numerical data A detailed presentation of the implementation method together with an extensive numerical test suite may be found in MZT 12 For compatibility reasons not all the algorithms available with the static grid version of PLUTO have been extended to the AMR version The AMR implementation of PLUTO is not compatible at present with e constrained transport MHD e finite difference schemes e the ShearingBox module 10 1 e the FARGO module e Super Time Stepping integration for dif
162. rho included by default are changed to represent a cut in log scale red colormap in the xz plane at the point coordinate y 1 1 in the y direction Note that the default for dbl of dbl h5 datasets should never be changed since restarting from a given file requires ALL variables being evolved in time CHAPTER 11 OUTPUT AND VISUALIZATION 98 Table 11 2 Output data formats and supported graphic visualization packages 11 3 Visualization PLUTO data files can be read with a variety of commercial and open source packages In what follows we describe how PLUTO data files can be read and visualized with IDI VisIf ParaView pyPLUTO 11 3 3 Mathematica and Gnuplof Table show some of the visualization softwares supporting different output formats We recall that reading of dbl or flt files must be complemented by grid information which is stored in a separate file grid out On the other hand VTK and HDF5 files xmf h5 vtk or hdf5 are stand alones in the sense that they embed grid information and can be opened alone 11 3 1 Visualization with Gnuplot Gnuplot can be used to visualize relatively small or moderately large 1 or 2D datasets written with the tabulated tab or binary data formats dbl or tity Gnuplot can be started at the command line by simply typing gt gnuplot In the following we give a short summary of the available options while a more detailed documen tation can be found in Do
163. rm X is an array of fractional abundances typically ion or molecule number fractions and S is a reaction source term The right hand side of Equations is implemented in the function Src Cooling lt COOLING gt radiat c of each corresponding cooling module except for the POWER_LAW cooling where integration is performed analytically The user can select one among several different cooling module by setting the COOLING flag during the python script e POWER_LAW power law cooling see e TABULATED only the equation for the internal energy with a tabulated cooling function A T is provided No chemical network see e SNEq cooling function for atomic hydrogen X X 1 including ionization recombination and collisionally excited emission lines e H2_COOL cooling function for atomic and molecular and atomic hydrogen X X o X 1 XaH11 including ionization recombination and collisionally excited emission lines e MINEg cooling function for atomic and molecular and atomic hydrogen treating the time dependent ionization state of the plasma X Xy XHe Xc Xn Xne Xo Xs see Cooling modules are implemented inside the Src Cooling directory and require three dimensional constants to be correctly initialized Dimensional constants are essential to scale data values to cgs physical units as explained in 95 1 1 Other variables are introduced to control crucial parameters such as the maximum allowed cooling rate in each time step or the cuto
164. rocessors used by mpirun or an error will occurr Use fname as initialization file instead of pluto ini Restart computations from the n th output file in HDF5 double preci sion format dbl h5 only for static grids The input data files are read from the directory specified by the output dir variables in pluto ini default is current working directory With Chombo AMR this switch is equivalent to restart Generate grid only do not start computations Stop computations after n steps Do not write data to disk Do not perform parallel domain decomposition along the x1 x2 or x3 direction respectively Restart computations from the n th output file in double in precision format dbl for static grid or Chombo checkpoint file chk nnnn hdf5 for Chombo AMR For the static grid input data files are read from the directory specified by the output dir variables in pluto ini de fault is current working directory Show domain decomposition when running in parallel mode Exclude from integration regions of zero pressure gradient that ex tends up to the end of the domain in the x1 x2 or x3 direction re spectively This option is specifically designed for jets propagating along one of the coordinate axis In parallel mode parallel decompo sition is not performed along the selected direction Set the grid resolution in the x1 direction to n1 zones by overriding pluto ini Cell aspect ratio is preserved by modifying the grid r
165. roner C D Munz T Schnitzer and M Wesenberg Hyperbolic Divergence Cleaning for the MHD Equations Journal of Computational Physics 175 2002 645 673 L Del Zanna O Zanotti N Bucciantini and P Londrillo ECHO a Eulerian conservative high order scheme for general relativistic magnetohydrodynamics and magnetodynamics A amp A473 2007 11 30 D Galli and F Palla The chemistry of the early Universe A amp A335 1998 403 420 T A Gardiner and J M Stone An unsplit Godunov method for ideal MHD via constrained transport Journal of Computational Physics 205 2005 509 539 A Harten B Engquist S Osher and S R Chakravarthy Uniformly high order accurate essentially non oscillatory schemes III Journal of Computational Physics 71 1987 231 303 D Hollenbach and C F McKee Molecule formation and infrared emission in fast interstellar shocks I Phys ical processes ApJS41 1979 555 592 G S Jiang and C W Shu Efficient Implementation of Weighted ENO Schemes Journal of Computational Physics 126 1996 202 228 124 Kle98 Ld04 Lio96 LL87 LOC94 Loh87 MB05 MBM 07 MFS 12 Mig07 Mig14 MKO05 MM96 MMO07 MPB05 MT10 MTB10 MUB09 MZT 12 OBR 08 Pow94 PRL 99 RBMF97 Roe81 Sal94 SG10 SH97 W Kley On the treatment of the Coriolis force in computational astrophysics A amp A338 1998 L37 L41 P Londrillo and L del Zanna On th
166. s v VX1 sqrt g_gammax v PRS v RHO in units of 1 Km s Mxcs in units of 15 cs 1 Em sec x Assign a magnetic field of 10 5 Gauss x v BX1 1 e 5 sqrt 4 0 CONST_PI UNIT_DENSITY UNIT_VELOCITY The CompEquil function is not strictly necessary but it has been introduced to compute ioniza tion equilibrium values for a given reference temperature and number density Its implementation may differ depending on the cooling module In alternative the fraction of neutrals could have been specified directly e g v X HI 0 6 Finally we notice that it is customary sometimes to assign magnetic field values in terms of the plasma 8 2p B Since 8 is already a dimensionless parameter one should not worry about proper dimensionalization and the line defining the magnetic field must be replaced by CHAPTER 5 INITIAL AND BOUNDARY CONDITIONS INIT C 47 4 0 this is my plasma beta 2p B 2 x sqrt 2 0 v PRS beta in units of v_01sqrt 4ipilrho_0 x 5 1 3 Assigning Initial Conditions from Input Files It is possible to assign initial conditions from user supplied binary data by providing i a grid data file and ii a single raw binary file containing the variables to be read The size dimensions and even the geometry of the input grid may be different from the actual grid employed by PLUTO as long as the coordinate transformation is supported This provides a flexible and efficient tool
167. s lossless data compression It requires libpng to be installed on your system Different images are associated with different variables and can have different sets of attributes de fined by the Image structure An image structure has the following customizable elements e slice plane a label X12_PLANE X13_PLANE X23_PLANE setting the slicing 2D plane e slice_coord a real number specifying the coordinate orthogonal to slice plane e max min the maximum and minimum values to which the image is scaled to If max min au toscaling is used CHAPTER 11 OUTPUT AND VISUALIZATION 95 e logscale an integer 0 or 1 specifying a linear or logarithmic scale e colormap the colormap Available options are red red map br blue red bw black and white blue blue green green In 2D the default is always slice plane X12_PLANE and slice_coord 0 Image attributes can be set independently for each variable in the function ChangeDumpVar in Src userdef_output c see q11 2 1 Default By default only density is written 11 1 6 The grid out output file The grid out file contains information about the computational grid used during the simulation It is an ASCII file starting with a comment header containing the creation date dimension and geometry of the erid PLUTO 4 2 Grid File Generated on lt date gt 4 4 DIMENSIONS lt DIMENSIONS gt GEOMETRY lt GEOMETRY gt X1 lt xl_
168. scale all quantities to non dimensional units in order to avoid occurrences of extremely small lt 1079 or large 10 2 numbers that may be misinterpreted by numerical algorithms For simple adiabatic simulations involving no source term the dimensionalization process can be avoided since the HD or MHD equations are scale invariant However dimensionalization is strictly necessary when specific length time or energy scales are introduced in the problem and they must compare to the dynamical advection scales For such problems PLUTO requires three fundamental units to be specified through the definitions of the following symbolic constants UNIT_DENSITY po sets the reference density in gr cm UNIT_LENGTH Lo sets the reference length in cm UNIT_VELOCITY vo sets the reference velocity in cm s All other units are derived from a combination of the previous three time is measured in units of to Lo vo pressure in units of po pov while magnetic field for the MHD module only see in units of Bo VoV AT po i e Pegs Ucgs Pegs Doge pP J U 9 P J B AA 5 1 po VO povs Ar povi If not specified the default values of UNIT_DENSITY UNIT_LENGTH UNIT_VELOCITY are respectively po 1 mp gr cm Lo 1 AU and vw 1 Km s The values of the three fundamental units can be changed by redefining them in your definitions h e g x Beg user defined constants do not change this line define UNIT
169. set color to black and white gnuplot gt splot data 0001 db1 bin array 200x200 format Sdouble If you have IDL installed you can easily display pressure from the last written output files with IDL gt pload 1 IDL gt display xl x1l x2 x2 prs Several other visualization options are described in more details in 911 3 You need to include PLUTO Tools IDL into your IDL search path 11 3 2 CHAPTER 0 QUICK START 8 0 4 Setting up your own test problem As an illustrative example we show how PLUTO can be configured to run a 2D Cartesian hydrody namic blast wave from scratch We assume that you have already followed the steps in 1 First in your home or work directory you need to create a folder which will contain the necessary files for the test For instance gt mkdir Blastwave gt cd Blastwave 2 You can now start the setup process by invoking the Python script to set dimensions geometry numerical scheme and so on Blastwave gt python SPLUTO_DIR setup py and select Setup problem from the main menu Using the arrows keys make the following changes set DIMENSIONS and COMPONENTS to 2 USER_DEF_PARAMETERS to 3 and leave the other fields as they are User defined parameters will be used later in the initial condition routine Press enter to confirm the changes and proceed to the following screen menu Since we don t have to change anything here you can press enter once more 3 We
170. sions Different builds are necessary for serial or parallel execution and since in both cases library names are the same by default it is advisable to store them in separate locations On a single processor machine serial libraries can be built for example using gt configure prefix usr local lib HDF5 serial gt make gt make check optional gt make install This will install the libraries under usr local lib HDF5 serial lib If you do not have root privileges choose a different location in your home directory e g PLUTO_DIR Lib HDF5 serial Note The I O of Chombo 3 2 has been updated to use the HDF5 1 8 API Howerer if HDF5 1 6 x is installed on your system the support for the 1 6 API is still provided by adding the DH5_USE_16_APT flag to the HDF INCFLAGS variable inside your Make defs local see 12 1 2 Nev ertheless it is not guaranteed that the HDF5 1 6 API will be supported in future Chombo releases On multiple processor architectures parallel libraries can be built by specifying the name of the mpicc compiler in the CC variable and invoking configure with the enable_parallel switch e g CC mpicc configure prefix usr local lib HDF5 parallel enable parallel bash user make make check optional gt gt gt gt make install This will install both shared dynamic so and static a libraries If you build shared libraries the environment variable LD_LIBRARY_PATH s
171. son and more information on their implementation with the MHD GLM formulation consult MTB10 11 Output and Visualization In this Chapter we describe the data formats supported by the static grid version of PLUTO and how they can be read and visualized with some popular visualization packages 11 1 Output Data Formats With the static version of PLUTO data can be dumped to disk in a variety of different formats The majority of them is supported on serial as well as parallel systems The available formats are classified based on their file extensions dbl double precision 8 byte binary data serial parallel flt single precision 4 byte binary data serial parallel dbl h5 double precision 8 byte HDF5 data serial parallel fith5 single precision 4 byte HDF5 data serial parallel _vik VTK legacy file format using structured or rectilinear grids serial parallel tab tabulated multi column ascii format serial only ppm portable pixmap color images of 2D data slices serial parallel png portable network graphicg color images of 2D data slices serial parallel Output files are named as base nnnn ext where base is either data when all variables are written to a single file or the name of the corresponding variable when each variable is written to a different file see Table 11 1 nnnn is a four digit zero padded integer counting the output number and ext is the corresponding fi
172. stant angular velocity 2 around the vertical polar axis z This feature should be enabled only when GEOMETRY is one of CYLINDRICAL POLAR or SPHERICAL The value of Q is specified using the global variable y OmegaZ inside your Init function The discretization of the angular momentum and energy equations is then done in a conservative fashion K1e98 MFS 12 For example in polar geometry we solve Op _ plug tw O p gt Ro vg wz V Rp ve w v 96 0 2 7 2 2 y E Sot wp V Fo Spo opos 0 where w RQ R is the cylindrical radius and F g is the standard energy flux and body force terms have been omitted only for the sake of exposition Note that the source term in the radial component of the momentum equation implicitly contains the Coriolis force and centrifugal terms w pu aL 2pvgQz P R 2 8 On the other hand the azimuthal component of the Coriolis force has been incorporated directly into the fluxes using the conservation form An example of such a configuration in polar or spherical geometry may be found in the directory Test_Problems HD Disk_Planet 2 2 7 THERMAL CONDUCTION Include thermal conduction effects see The available options are e NO thermal condution is not included e EXPLICIT thermal conduction is treated explicitly 98 4 1 e SUPER_TIME_STEPPING thermal conduction is treated using super time stepping 48 4 2 2 2 8 VISCOSITY Include viscous terms see Options a
173. t this type at this moment This particular boundary condition can be used only if the ShearingBox module described in is enabled userder User supplied boundary conditions it requires coding your own boundary conditions in the func tion UserDefBoundary in init c see 5 2 Like the Grid block you should include the x3 boundaries for 2D runs even if they will not be consid ered CHAPTER 4 RUNTIME INTITALIZATION FILE PLUTO INI 39 4 6 The Static Grid Output Block This block controls the output options in the static grid version of the code Output files are written at specific times in a specific directory local working directory by default using different file formats described in Chapter 11 The available fields are uservar integer stringl string2 Defines supplementary variables to be written to disk in any of the format described below The first integer represents the number of supplementary variables When greater than zero it must be followed by as many variable names separated by spaces see Chapter 11 output_dir string Sets the directory name for writing and reading simulation data When writing any of the data formats available below and the corresponding out log files will be written in the directory spec ified by output dir The pluto log log file only if PRINT_TO_FILE has been set to YES will also be written to the same directory When reading during restarts dbl or dbl h5 files and
174. tch only for the IDEAL or TAUB EOS the equation of entropy is added to the system of conservation laws 00 Ot where 0 po for HD or MHD or oe Do for relativistic flows o is the entropy e g o p p for the IDEAL EOS and S is a source term accounting for dissipative terms in the HD or MHD modules see Eq 3 27 of BS03 V ov So 2 3 Ov Ox j l Equation is solved in a conservative way and entropy is essentially treated as a passive scalar At each integration step gas pressure can be recovered selectively from total energy or conserved entropy in the following fashion gt p v RVD 937 pl 2 4 l p p E if F o 0 2 5 P p Fe if F 0 1 Here F 0 tells if a zone has been flagged with FLAG_ENTROPY at the beginning of the step and it depends on the value assigned to ENTROPY_SWITCH e NO the entropy equation is not included eo SELECTIVE by enabling this option a selective update is employed By default all zones are flagged to be updated using the entropy equation unless they lie in proximity of a shock wave Thus for each zone we evaluate 0 if V v lt 0 and IVpl p gt e w Vll p gt amp ae 1 otherwise where V is a three point undivided difference operator and e sets the shock strength threshold see EPS_PSHOCK ENTROPY in Appendix B 3 The first criterion acts as shock detector The implementation can be found in the source file Src flag_shock c Note that by ena
175. the corresponding out must be present in this directory or an error will occur If output di r is not specified the directory from which pluto is executed is assumed dbl double integer string string Assigns the output intervals for double precision 8 bytes binary data A negative value sup presses output with this format The first field double specifies the time interval in code units between consecutive out puts The second field can be an integer giving the number of steps between consecutive outputs ora string giving the wall clock time between consecutive outputs A value for instance of 2 40h tells PLUTO to write one dbl file every two hours and 40 minutes Negative values will be ignored for this control parameter The last field st ring can be either single_file one single output file per time step con taining all of the variables or mult iple_files different variables are written to different files When asynchronous I O is enabled 43 1 1 a third option single_file_async is permitted for flt or dbl binary files to specify that asynchrounous binary output has to be performed Double precision format files can be used to restart the code using the restart n command line argument flt double integer string string cgs Like db1 but for single precision 4 bytes data files The last value cgs is optional and can be given to save datafiles directly in cgs physical units rather th
176. tion e g 02 and copy the corresponding configuration files i e PLUTO Test_Problems MHD Orszag_Tang gt cp definitions_02 h definitions h PLUTO Test_Problems MHD Orszag_Tang gt cp pluto_02 ini pluto ini 3 Run the Python script PLUTO Test_Problems MHD Orszag_Tang gt python SPLUTO_DIR setup py select Setup problem and choose the default setting by pressing enter 4 Once you return to the main menu select Change makefile and choose a suitable makefile e g Linux gcc defs and press enter 5 Exit from the main menu Quit or press q Edit pluto ini and under the Grid block lower the resolution from 512 to 200 in both directions X1 grid and X2 grid Change single_file in the dbl output under the Uniform Grid Output block to multiple_files Finally edit definitions h and change PRINT_TO_FILE from YES to NO 6 Compile the code PLUTO Test_Problems MHD Orszag_Tang gt make 7 If compilation was successful you can now run the code by typing PLUTO Test_Problems MHD Orszag_Tang gt pluto At this point PLUTO reads the initialization file pluto ini and starts integrating The run should take a few minutes depending on the machine you re running on and the integration log should be dumped to screen You can display data e g density with Gnuplot v 4 2 or higher from the last written file using gnuplot gt set pm3d map set map style drawing gnuplot gt set palette gray
177. tion of the diffusion terms thermal conduction viscosity and magnetic resistivity and cooling 6 1 The HD Module The HD module implements and solves the Euler or Navier Stokes equations of classical fluid dynamics The relevant source files and definitions for this module can be found in the Src HD directory 6 1 1 Equations With the HD module PLUTO evolves in time following system of conservation laws T p pv 0 O a m V mv pl pVO pg 6 1 E p E p p v m g where pis the mass density m pv is the momentum density v is the velocity p is the thermal pressure and F is the total energy density 2 m E pe 20 6 2 An equation of state provides the closure pe pe p p The source term on the right includes contributions from body forces and is written in terms of the time independent gravitational potential and and the acceleration vector g 45 3 The right hand side of the system of Eqns is implemented in the RightHandSide function inside Src MHD rhs c employing a conservative discretization that closely follows the expression given in A 1 2Jand A 1 3 for Cartesian polar and spherical geometries without magnetic fields Primitive variables are defined by V p v p where v m p while p p p pe depends on the equation of state see Chapter 7 The maps U V and its inverse are provided by the functions PrimToCons and ConsToPrim Primitive variables are generally more con
178. tion of it This can be accomplished by specifying sub domain through the keywords xl range x2range and x3range For instance IDL gt PLOAD hdfb 5 lev 6 xlrange 0 25 0 75 x2range 0 715 1 25 will load data 0005 2d hdf5 ref level 6 T Dut only inside the region x in 0 45 0 75 y im 10 75 L 23 IDL gt DISPLAY x1 x1 x2 x2 rho nwin 1 imax 1 1 imin 0 65 IDL gt OPLOTBOX ctab 3 CHAPTER 12 ADAPTIVE MESH REFINEMENT AMR 112 12 4 2 Visualization with VisIt VisIt can read Chombo HDF5 datafiles Individual hdf5 files or databases can be opened and visualized from the GUI in exactly the same way as vtk of h5 files and level plots can be over imposed from Add Subset Levels If you are using curvilinear coordinates or cartesian coordinates with an origin offset i e the do main s lower corner 0 0 0 and or different grid spacings along different directions the correct coordinate transformation can be done by applying the Displacement operator Example e Select a valid data hdf5 database by clicking on Open e Add Pseudocolor rho e Operators gt Transform Displace e Click on the Displace operator to set the attributes Displacement variable gt Default gt Vectors Displacement As an example we show in Fig a close up of the final solution obtained with configuration 08 of the Test_Problems HD Disk_Planet problem 0 1149 Max 55 19 2 Min 0 1149 Subset Var levels
179. tly describe discontinuities the advantages offered by their employment are more evident in presence of smooth flows 10 3 1 WENO schemes The WENO schemes are based on the essentially non oscillatory ENO schemes originally developed by using a finite volume formulation and later improved by into a finite difference form Unlike TVD schemes that degenerate to first order at smooth extrema ENO schemes main tain their accuracy successfully suppressing spurious oscillations This is accomplished utilizing the smoothest stencil among a number of candidates to compute fluxes at the cell faces WENO schemes are the natural evolution of ENO schemes where a weighted average is taken from all the stencil candidates Weights are adjusted by local smoothness indicators Originally de veloped by for 1 D finite volume formulation WENO schemes were then implemented in multi dimensional FD by JS96 optimizing the original weighing for accuracy Currently the available WENO schemes in PLUTO are the 5 order WENOZ of which improves over the original one in that it is less dissipative and provide better resolution at critical points at a very modest additional computational cost A third order WENO scheme is also provided namely WENO 3 of YC09 More details can be found in the paper by Mignone Tzeferacos amp Bodo IMTB10 CHAPTER 10 ADDITIONAL MODULES 91 Figure 10 2 Long term numerical de cay of a circularly polarized Alfven wave on a 2D peri
180. to assign initial conditions by mapping data values originally defined on different computational domains For instance you can map a 2D spherical grid onto a 2D axisymmetric cylindrical domain generate a 3D Cartesian domain by rotating any 2D axisymmetric data and so forth The module is initialized by calling the function InputDataSet which reads and stores input grid information such as size number of variables geometry and dimensions This function should be called only once from your Init function InputDataSet grid_file input_var where the first argument grid_file isa string giving the name of the input grid file while input_var is an array of integers specifying which variables are contained in the input data file The input grid file should be written using the same format employed by PLUTO see After initialization any subsequent call to InputDataRead data_file string will read and store into memory the values of the variables contained in the input data file data_file The second argument specifies the endianity of the input data files i e string little big An empty string does not change anything Please note that unless the input data file is changed this function should also be called only once The data file should be written using binary format using either single or double precision with extensions flt for the former or dbl for the latter Variables should be stored sequentially
181. ts u is defined by the GetMu function CHAPTER 7 EQUATION OF STATE 70 Figure 7 1 Density plot for the Sod shock tube test att 0 2 obtained with the IDEAL EOS red and the PVTE_LAW EOS green with reference density 10 mp and reference velocity 10 cm s en eee e 006 Figure 7 2 Equivalent T p pe 1 for T K the PVTE_LAW EOS of a partially ionized hy drogen gas 7 3 2 Analytic vs tabulated approach As PLUTO performs conversions between primitive e g density and pressure and conservative vari ables e g total energy and momentum during a single update step the PVTE_LAW Eqns 7 8 must often be inverted to obtain the temperature from pressure or internal energy Since Eqns 7 8 specially the second one can be nonlinear functions of T the inversion must be taken numerically using a root finder and this can be an expensive task In the equilibrium case however a faster and often more convenient approach is to have PLUTO pre compute tabulated versions of the EOS so as to replace expensive function evaluations with tables In this case no root finder is used and computations involving EOS require direct or inverse simpler lookup table operations and cubic linear interpolation see 3 2 2 of VMBM15 This feature is always turned on by default but can be overridden through the user defined constants PV_TEMPERATURE_TABLE and or TV_ENERGY_TABLE in your definitions h see 92 3 For example the following definiti
182. uilibrium case there s no chemical network and fractions are not evolved independently but are computed when necessary using some sort of equilibrium assumptions such as Saha LTE valid in the high density limit or collisional ionization equilibrium CIE valid at low densities This cor responds to express fractions as X X T p so that quantities depending on X become functions of T p For example the thermal EOS becomes p p p T while internal energy becomes a func tion of two variables e e T p In this case the inverse functions T T p p and T T e p are computed by finding the roots of nonlinear equations The implementation of the PVTE_LAW EOS can be found in the Src EOS PVTE directory The source file pvte_law c or pvte_law_template c if you are starting from scratch provides the interface between the user implementation and the module through the following functions e InternalEnergyFunc compute and return the internal energy density pe where e e T X in non equilibrium chemistry or e e T p in the equilibrium case e GetMu compute the mean molecular weight u p X or u T p e Gammal compute the value of the first adiabatic index Fee log u XT p 1 o log T logT ie EA e where 7 9 cy pT _ Ologp _ _ logu a Dlogp gt O log p For a thermally ideal gas it can be shown that the specific internal energy e is a function of the temperature T and chemi
183. uire a reduction of the parabolic Courant number Cp see 48 4 and a more tight coupling be tween the hydrodynamical and thermal conduction scale The latter condition may be accom plished by lowering the rmax_par parameter 44 3 which controls the ratio between the cur rent time step and the diffusion time scale see also An example problem can be found in Test Problems MHD Thermal conduction Blast 8 3 1 Dimensions Equations 8 6 8 8 are solved in dimensionless form by expressing energy and time in units of p9ve and Lo vo respectively and by writing temperature as T p p K ys where p and p are in code units and pis CHAPTER 8 DISSIPATIVE EFFECTS 76 the mean molecular weight Here po vo Lo are the unit density velocity length while K is the KELVIN macro see 45 1 1 The thermal conduction coefficients must be properly defined by re absorbing the correct normalization constants in the TC_kappa function as follows My 8 9 povoLokB where for instance one may use Kegs y 5 6 1077T5 2 and keys 1 3 3 10716n VT Bss both in units of ergs7 K cm 1 while B2 4mpov B An example of such dimensionalization can be found in Test_Problems MHD Thermal_Conduction Blast CHAPTER 8 DISSIPATIVE EFFECTS 77 8 4 Numerical Integration of Diffusion Terms 8 4 1 Explicit Time Stepping With the explicit time integration parabolic contributions are added to the upwind hyperbolic fluxes at the same time in
184. ules algorithms and Riemann solvers in the code 7 3 The PVTE_LAW Equation of State The PVTE Pressure Volume Temperature Energy EOS allows the user to specify the internal energy as a general function of the temperature T and chemical fractions or concentrations X as described in the paper by VMBM15 The thermal EOS together with the caloric EOS link the five quantities p p T X and e and are used by the code to compute two of them given the remaining three a ee TS 7 8 e e T X where m is the atomic mass unit and u X the mean molecular weight depends on the gas compo sition The PVTE_LAW EOS allows the user to provide explicit definitions for u X and e T X ina thermodynamically consistent wa The implementation of this EOS depends on how chemical fractions are computed and a major dis tinction should be made between non equilibrium and equilibrium cases Non equilibrium case a chemical network is used to evolve X t through rate equations under non equilibrium conditions see This occurs for example when this EOS is used in conjunction with a cooling module that includes a time dependent reaction network In this case species are evolved independently and their value is at disposal when performing conversion between pressure temperature and internal energy In particular recovering temperature from internal energy T T e X requires inverting a nonlinear equation by means of an iterative root finder Eq
185. undary condition for a jet model The computational domain is a 2D box in cylindrical geometry so that x R x2 z A constant inflow is prescribed a the jet nozzle located at the z 0 boundary for R lt 1 while reflective boundary conditions are assigned for R gt 1 The inflow values are specified as P 1 p R x p R z UR 0 uR R z uR R z for R lt 1 for R gt 1 i M v R z v R z P 1 T p R 2 p R z where M is the Mach number void UserDefBoundary const Data d RBox int side Grid grid int le Je Ky ny double xl xX2 X3 xl grid IDIR x array pointer to xl coordinate x2 grid JDIR x array pointer to 22 coordinate x3 grid KDIR x array pointer to x3 coordinate lt if side X2_BEG select the boundary side x BOX_LOOP box k j i Loop over boundary zones x if xl i lt 1 set jet values for r lt 1 x d gt Vc RHO 1 0 d gt Vc iVR 0 0 d gt Vc ivZ g_inputParam MACH d gt Vc PRS 1 0 gmm Jelse x reflective boundary for d gt Vc RHO d gt Vc RHO k 2 IBEG j d gt Vc iVR d gt Vc iVR 2xJBEG jJ k d gt Vc ivVZ i d gt Vc iVZ k 2 JBEG k d gt Vc PRS j d gt Vc PRS 2xJBEG The previous piece of code is executed only if you have selected userdef at the X2 beg boundary inside your pluto ini The macro BOX_LOOP box k j i performs
186. uz Arguments e v a pointer to a vector of primitive variables e x1 x2 x3 local spatial coordinates e xnul a pointer to the 1st viscous coefficient shear e xnu2 a pointer to the 2nd viscous coefficient bulk Even though the behaviour of these coefficients is arbitrary according to the user s needs for mono atomic gases Molecular Theory gives v 0 The coefficient of shear viscosity v on the other hand is usually specified with a power law behaviour with respect to the temperature e g the Sutherland formula For more information on the analytical and numerical treatment of viscosity see LL87 and Tor97 It should be noted nonetheless that both transport coefficients must have dimensions of p x length time for the correct control of the timestep according to the stability condition discussed at the beginning of this chapter CHAPTER 8 DISSIPATIVE EFFECTS 74 8 2 Resistivity The resistive module is enabled by setting RESISTIVITY to either EXPLICIT for time explicit com putations or to SUPER_TIME_STEPPING to accelerate explicit computations from the Python menu Magnetic field dissipation is modeled by introducing the resistivity tensor 7 so that the electric field be comes v x B n J where J V x B is the current density The induction and energy equations gain extra terms on the right hand sides OB Vx vx B Vx n J Ot 95 8 3 a FVW E p B v B V n J x B Similarly t
187. vels elements or an error will result Refine_thresh double Sets the threshold value x above which cells are tagged for refinement during the grid generation process see 412 3 When x U gt x Refine thresh the cell is tagged to be refined Tag_buffer_size integer Sets the amount by which to grow tags as a safety factor before passing to MeshRefine Block_factor integer Sets the number of times that grids will be coarsenable by a factor of 2 A higher number produces blockier grids Max_grid_size integer Sets the largest allowable size of a grid in any direction Any boxes larger than that will be split up to satisfy this constraint Fill_ratio double A real number between 0 and 1 used to set the efficiency of the grid generation process Lower number means that more extra cells which are not tagged for refinement wind up being refined along with tagged cells The tradeoff is that higher fill ratios lead to more complicated grids and the extra coarse fine interface work may outweigh the savings due to the reduced number of fine level cells CHAPTER 4 RUNTIME INTITALIZATION FILE PLUTO INI 36 4 3 The Time Block This section specifies some adjustable time marching parameters CFL double Courant number it controls the time step length and in general it must be less than 1 The actual limit can be inferred from Table In the case of unsplit Runge Kutta time stepping for instance CFL S 1 Naim where N
188. venient and preferred when assigning initial boundary conditions and in the interpolation algorithms The vector of primitive quantities V obeys the quasi linear form of the equations Op ap UV poe pro 0 Ov Vp q eo O Lao Vp p2ev v TL Ot where cs yT p p is the adiabatic speed of sound for an ideal EOS The quasi linear form is implemented in the Src HD prim_eqn c source file and it is required during the predictor stages of the HANCOCK or CHARACTERISTIC_TRACING time stepping schemes V g 6 3 57 CHAPTER 6 BASIC PHYSICS MODULES 58 6 2 The MHD Module The MHD module is suitable for the solution of ideal or resistive non relativistic magnetohydrody namical equations Source and definition files are located inside the Src MHD directory 6 2 1 Equations With the MHD module PLUTO solves the following system of conservation laws p pu j 0 m mv BB lp pV pg aa HV 6 4 t E p E p pw B v B m g B vB Bv 0 where p is the mass density m pv is the momentum density v is the velocity p p B 2 is the total pressure thermal magnetic B is the magnetic field land E is the total energy density 2 2 bon 6 5 20 2 where an additional equation of state provides the closure pe pe p p see Chapter 7 The source term on the right includes contributions from body forces and is written in terms of the time independent gravitational potential and and the acce
189. way from the origin The default value is No CHOMBO_CONS_AM YES NO Sets the name of the conservative variable used by Chombo when tagging zones for refinement Allowed values are RHO for density ENG for total energy MX1 for normal com ponent of momentum etc The default value is total en ergy density or density when there s no energy equation The special value 1 can be given to supply a user defined variable instead e g pressure or kinetic energy using the CHOMBO_REF_VAR lt vname gt computeRefVar function See 12 3 for more detail No tice that since CHOMBO_REF_VAR is one of the conservative vari ables used to perform prolongation and restriction operations if CHOMBO_CONS_AM is set to YES see B 3 iMPHI stands for the conserved angular momentum Important owing to the different type of conserved variable names the CHOMBO_REF_VAR should never be used inside pre processor conditional statements Continued on next page 120 Table B 1 Continued from previous page Name CAIR REF STATE EPS PSHOCK ENTROPY EPS PSHOCK FLATTEN FARGO AVERAGE VELOCITY FARGO NSTEP AVERAGE FARGO_ORDER GLM_ALPHA GLM_EXTENDED H_MASS_FRAC He_MASS_FRAC Continued on next page Value 27273 real real YES NO int int real YES NO double double Description Set the reference state used during the characteristic tracing step see Src States char_traci
190. whereas to exclude staggered magnetic field components if any from the fIt h5 format CHAPTER 11 OUTPUT AND VISUALIZATION 94 11 1 3 VTK Output vtk data format VTK from the Visualization ToolKit format output follows essentially the same conventions used for the dbl or flt outputs Single or multiple VTK files can be written by specifying either single file or multiple_files in your pluto ini and data values are always written using single precision with byte order set to big endian The mesh topology uses a rectilinear grid format for CARTESIAN or CYLINDRICAL geometry while a structured grid format is employed for POLAR or SPHERICAL geometries Data is written with the CELL_DATA attribute and grid nodes or vertices are used to store the mesh The following symbolic constants 2 3 can be used to control some options of the output vtk files e VIK_TIME_INFO when set to YES time information will be added to the header section of the Vik file Beware that standard VTK files do not have a specific construct for adding time informa tion and by doing so this information will be available only to the VisIt visualisation software see which implements a convention where CYCLE and TIME values can be specified as FieldData in the file If you re using Paraview or other visualisation software different from Visit enabling this option will most likely result in a software crash e VIK_VECTOR_DUMP by default all flow quantities e g dens

User's Guide

Contents

Download Pdf Manuals

Related Search

Related Contents