Good day Flash community:<br><br>Two questions. I'll be really thankful if you had any hints about them.<br><br>1) I'm trying to complie flash3.1 for a uniform grid with no fixed block size:<br><br>./setup mhd_ug_nofbs -3d +usm -auto -opt +nofbs -site=apgrid<br>
<br>and I get:<br><br><i>-L/mraosw/data1/mh475/krause/lib/hdf/5-1.6.5-amd-icc/lib -lhdf5 -lz -lhdf5_fortran -L/mraosw/data1/mh475/krause/lib/mpich-1.2.6/lib -lmpich <br>/home/krause/dataw/lib/mpich-1.2.4/lib/libmpich.a(p4_secure.o)(.text+0x91): In function `start_slave':<br>
: warning: Using 'getpwuid' in statically linked applications requires at runtime the shared libraries from the glibc version used for linking<br>/mraosw/data1/mh475/krause/lib/hdf/5-1.6.5-amd-icc/lib/libhdf5.a(H5FDstream.o)(.text+0x723): In function `H5FD_stream_open_socket':<br>
: warning: Using 'gethostbyname' in statically linked applications requires at runtime the shared libraries from the glibc version used for linking<br>io_h5file_interface.o(.text+0x6d): In function `io_h5init_file_':<br>
: undefined reference to `H5Pset_fapl_mpio'<br>io_h5file_interface.o(.text+0x12f): In function `io_h5open_file_for_read_':<br>: undefined reference to `H5Pset_fapl_mpio'<br>make: *** [flash3] Error 1</i><br><br>
I've tried, unsuccessfully, to edit some parts of the cluster's Makefile.h.<br><br><br>2) I'm implementing a gas with a random magnetic distribution. I generate the magnetic fields outside flash, making sure their divb~1.e-14. I then read these fields into Flash3.1 and they work well for a uniform grid or with the AMR with lrefine_max=lrefine_min. However, if I give lrefine_max different from lrefine_min, divb appears to be fine at time=0, but after the first timestep it gets jumps up to ~ 0.1 at cells adjacent to refinement jumps only. I'm attaching important bits of info below.<br>
<br>Thanks a lot.<br><br>./setup MySims/mhd_01 -3d +usm -auto -opt -maxblocks=400 -site=darwin<br><br>------------------------<br><div style="margin-left: 40px;"><i>subroutine Simulation_initBlock(blockID, myPE, F) ! following the setup from the Fields loop test problem</i><br>
<br><i> use Simulation_data, ONLY : sim_gCell, sim_gamma, &</i><br><i> sim_smallX, sim_smallP, &</i><br><i> sim_beta, sim_B_cluster, &</i><br>
<i> sim_killdivb, sim_rc, &</i><br><i>sim_xMin, sim_xMax, sim_yMin, &</i><br><i>sim_yMax, sim_zMin, sim_zMax, &</i><br><i>sim_n_cube</i><br><br><i>! use Grid_data, ONLY : gr_nBlockX, gr_nBlockY, gr_nBlockZ</i><br>
<br><i>! use tree, ONLY : lrefine_max, lrefine_min</i><br><br><i> use Grid_interface, ONLY : Grid_getBlkIndexLimits, &</i><br><i> Grid_getCellCoords, &</i><br><i> Grid_getBlkPtr, &</i><br>
<i> Grid_releaseBlkPtr</i><br><br><i> use Driver_interface, ONLY : Driver_abortFlash</i><br><i> </i><br><i> implicit none</i><br><br><i>#include "constants.h"</i><br><i>#include "Flash.h"</i><br>
<br><i> !!$ Arguments -----------------------</i><br><i> integer, intent(in) :: blockID, myPE</i><br><i> !!$ ---------------------------------</i><br><br><i> integer :: i, j, k, n, istat, sizeX, sizeY, sizeZ</i><br><i> integer, dimension(2,MDIM) :: blkLimits, blkLimitsGC</i><br>
<i> real :: enerZone, ekinZone, eintZone</i><br><i> real :: rot, radius, dx, dy, dz, r0, taper, vel_factor</i><br><i> real, allocatable,dimension(:) :: xCoord,xCoordL,xCoordR,&</i><br><i> yCoord,yCoordL,yCoordR,&</i><br>
<i> zCoord,zCoordL,zCoordR</i><br><i> real, dimension(MDIM) :: del</i><br><i> real, pointer, dimension(:,:,:,:) :: solnData, facexData, faceyData, facezData</i><br><i> real :: xx,yy,zz</i><br>
<i>#ifdef FIXEDBLOCKSIZE</i><br><i> real, dimension(GRID_IHI_GC+1,GRID_JHI_GC+1,GRID_KHI_GC+1) :: Az,Ax,Ay</i><br><i> </i><i>#else </i><br><i> real, allocatable, dimension(:,:,:) :: Az,Ax,Ay</i><br><i> </i><i>#endif </i><br>
<br><i>logical, save :: once=.true., once2=.true., VecPotA=.true.</i><br><i> integer :: ii, jj, kk</i><br><i> real :: B0, x_slope, y_slope, z_slope, &</i><br><i> idl_x_ini, idl_x_end, &</i><br>
<i> idl_y_ini, idl_y_end, &</i><br><i> idl_z_ini, idl_z_end, &</i><br><i>largestCell, minimumX, maximumX, &</i><br><i> minimumY, maximumY ,&</i><br><i> minimumZ, maximumZ, &</i><br>
<i> stdev_b, max_b, mean_pres, mean_b, norm_b<br><br></i><i style="color: rgb(255, 0, 0);"> real, dimension(3,sim_n_cube,sim_n_cube,sim_n_cube), intent(in) :: F <br> !contains the previously-generated magnetic fields<br>
</i><i>!!</i><br><br><br><i> ! dump some output to stdout listing the paramters</i><br><i> if (myPE == MASTER_PE) then</i><br><i>1 format (1X, 1P, 4(A7, E13.7, :, 1X))</i><br><i>2 format (1X, 1P, 2(A7, E13.7, 1X), A7, I13)</i><br>
<i> endif</i><br><br><i> call Grid_getBlkIndexLimits(blockId,blkLimits,blkLimitsGC)</i><br><br><i> sizeX = blkLimitsGC(HIGH,IAXIS)-blkLimitsGC(LOW,IAXIS)+1</i><br><i> sizeY = blkLimitsGC(HIGH,JAXIS)-blkLimitsGC(LOW,JAXIS)+1</i><br>
<i> sizeZ = blkLimitsGC(HIGH,KAXIS)-blkLimitsGC(LOW,KAXIS)+1</i><br><br><i> allocate(xCoord(sizeX), stat=istat)</i><br><i> allocate(xCoordL(sizeX),stat=istat)</i><br><i> allocate(xCoordR(sizeX),stat=istat)</i><br><br>
<i> allocate(yCoord(sizeY), stat=istat)</i><br><i> allocate(yCoordL(sizeY),stat=istat)</i><br><i> allocate(yCoordR(sizeY),stat=istat)</i><br><br><i> allocate(zCoord(sizeZ), stat=istat)</i><br><i> allocate(zCoordL(sizeZ),stat=istat)</i><br>
<i> allocate(zCoordR(sizeZ),stat=istat)</i><br><br><i> xCoord = 0.0</i><br><i> xCoordL = 0.0</i><br><i> xCoordR = 0.0</i><br><br><i> yCoord = 0.0</i><br><i> yCoordL = 0.0</i><br><i> yCoordR = 0.0</i><br><br><i> zCoord = 0.0</i><br>
<i> zCoordL = 0.0</i><br><i> zCoordR = 0.0</i><br><br><i>#ifndef FIXEDBLOCKSIZE</i><br><i> if (NDIM == 2) then</i><br><i> allocate(Ax(sizeX+1,sizeY+1,1),stat=istat)</i><br><i> allocate(Ay(sizeX+1,sizeY+1,1),stat=istat)</i><br>
<i> allocate(Az(sizeX+1,sizeY+1,1),stat=istat)</i><br><i> </i><i> elseif (NDIM == 3) then</i><br><i> allocate(Ax(sizeX+1,sizeY+1,sizeZ+1),stat=istat)</i><br><i> allocate(Ay(sizeX+1,sizeY+1,sizeZ+1),stat=istat)</i><br>
<i> allocate(Az(sizeX+1,sizeY+1,sizeZ+1),stat=istat)</i><br><i> </i><i> endif</i><br><i>#endif</i><br><br><br><i> if (NDIM == 3) then</i><br><i> call Grid_getCellCoords(KAXIS,blockId,CENTER, sim_gCell,zCoord, sizeZ)</i><br>
<i> call Grid_getCellCoords(KAXIS,blockId,LEFT_EDGE, sim_gCell,zCoordL,sizeZ)</i><br><i> call Grid_getCellCoords(KAXIS,blockId,RIGHT_EDGE,sim_gCell,zCoordR,sizeZ)</i><br><i> endif</i><br><i> if (NDIM >= 2) then</i><br>
<i> call Grid_getCellCoords(JAXIS,blockId,CENTER, sim_gCell,yCoord, sizeY)</i><br><i> call Grid_getCellCoords(JAXIS,blockId,LEFT_EDGE, sim_gCell,yCoordL,sizeY)</i><br><i> call Grid_getCellCoords(JAXIS,blockId,RIGHT_EDGE,sim_gCell,yCoordR,sizeY)</i><br>
<i> endif</i><br><br><i> call Grid_getCellCoords(IAXIS,blockId,CENTER, sim_gCell,xCoord, sizeX)</i><br><i> call Grid_getCellCoords(IAXIS,blockId,LEFT_EDGE, sim_gCell,xCoordL,sizeX)</i><br><i> call Grid_getCellCoords(IAXIS,blockId,RIGHT_EDGE,sim_gCell,xCoordR,sizeX)</i><br>
<br><i> call Grid_getDeltas(blockID,del)</i><br><i> dx = del(1)</i><br><i> dy = del(2)</i><br><i> dz = del(3)</i><br><br><i> !------------------------------------------------------------------------------</i><br><i> ! Construct Az at each cell corner</i><br>
<i> ! Bx = dAz/dy - dAy/dz</i><br><i> ! By = dAx/dz - dAz/dx</i><br><i> ! Bz = dAy/dx - dAx/dy</i><br><i> Az = 0.</i><br><i> Ax = 0.</i><br><i> Ay = 0.<br><br> x_ini = 1. <br> y_ini = 1. <br> z_ini = 1. <br> x_end = real(sim_n_cube) <br>
y_end = real(sim_n_cube) <br> z_end = real(sim_n_cube) <br> <br> minimumX = sim_xMin - 4.d0*dx <br> minimumY = sim_yMin - 4.d0*dy <br> minimumZ = sim_zMin - 4.d0*dz <br> maximumX = sim_xMax + 4.d0*dx <br>
maximumY = sim_yMax + 4.d0*dy <br> maximumZ = sim_zMax + 4.d0*dz </i><br><i> </i><br><i> x_slope = (x_end-x_ini)/(maximumX-minimumX)</i><br><i> y_slope = (y_end-y_ini)/(maximumY-minimumY) </i><br>
<i> z_slope = (z_end-z_ini)/(maximumZ-minimumZ) </i><br><i> </i><br><i></i><i> do k = blkLimitsGC(LOW,KAXIS),blkLimitsGC(HIGH,KAXIS)+1</i><br><i> do j = blkLimitsGC(LOW,JAXIS),blkLimitsGC(HIGH,JAXIS)+1</i><br>
<i> do i = blkLimitsGC(LOW,IAXIS),blkLimitsGC(HIGH,IAXIS)+1</i><br><br><br><i> ! x Coord at cell corner</i><br><i> if (i <=blkLimitsGC(HIGH,IAXIS)) then</i><br><i> xx = xCoordL(i)</i><br>
<i> else</i><br><i> xx = xCoordR(i-1)</i><br><i> endif</i><br><br><i> ! y Coord at cell corner</i><br><i> if (j <=blkLimitsGC(HIGH,JAXIS)) then</i><br>
<i> yy = yCoordL(j)</i><br><i> else</i><br><i> yy = yCoordR(j-1)</i><br><i> endif</i><br><br><i> ! z Coord at cell corner</i><br><i> if (k <=blkLimitsGC(HIGH,KAXIS)) then</i><br>
<i> zz = zCoordL(k)</i><br><i> else</i><br><i> zz = zCoordR(k-1)</i><br><i> endif</i><br><br><i></i><br><i></i><i> ii = nint( x_slope*(xx-minimumX)+idl_x_ini ) </i><br>
<i> jj = nint( y_slope*(yy-minimumY)+idl_y_ini )</i><br><i> kk = nint( z_slope*(zz-minimumZ)+idl_z_ini )</i><br><br><i></i><i>!--------------------------------------------------------------</i><br>
<i>if ( (ii.lt.idl_x_ini).or.(ii.gt.idl_x_end).or.&</i><br><i> (jj.lt.idl_y_ini).or.(jj.gt.idl_y_end).or.& </i><br><i> (kk.lt.idl_z_ini).or.(kk.gt.idl_z_end) ) &</i><br><i></i><i>call Driver_abortFlash("*** init_block, IDL-Flash coords. error ***")</i><br>
<br><i>!!--------------------------------------------------------------</i><br><br><div style="margin-left: 80px;"><font color="#ff0000">!Read the outise-fields into flash</font><br></div><i> Ax(i,j,k) = F(1,ii,jj,kk)</i><br>
<i> Ay(i,j,k) = F(2,ii,jj,kk)</i><br><i> Az(i,j,k) = F(3,ii,jj,kk)</i><br><i> </i><br><i> enddo</i><br><i> enddo</i><br><i> enddo</i><br><br><br><i></i><i>!Initial conditions:</i><br>
<br><i> call Grid_getBlkPtr(blockID,solnData,CENTER)</i><br><br><i>#if NFACE_VARS > 0</i><br><i> if (sim_killdivb) then</i><br><i> call Grid_getBlkPtr(blockID,facexData,FACEX)</i><br><i> call Grid_getBlkPtr(blockID,faceyData,FACEY)</i><br>
<i> if (NDIM == 3) call Grid_getBlkPtr(blockID,facezData,FACEZ)</i><br><i> endif</i><br><i>#endif</i><br><br><i></i><br><i> ! Loop over cells within the block.</i><br><i> do k = blkLimitsGC(LOW,KAXIS),blkLimitsGC(HIGH,KAXIS)</i><br>
<i> do j = blkLimitsGC(LOW,JAXIS),blkLimitsGC(HIGH,JAXIS)</i><br><i> do i = blkLimitsGC(LOW,IAXIS),blkLimitsGC(HIGH,IAXIS)</i><br><br><i> solnData(SPECIES_BEGIN,i,j,k)=1.0e0-(NSPECIES-1)*sim_smallX</i><br>
<i> do n=SPECIES_BEGIN,SPECIES_END</i><br><i> solnData(n,i,j,k)=sim_smallX</i><br><i> enddo</i><br><i> </i><br><br><i></i><i> solnData(DENS_VAR,i,j,k)= 1.</i><br><i> solnData(PRES_VAR,i,j,k)= solnData(DENS_VAR,i,j,k)/sim_gamma </i><br>
<br><br><i> solnData(VELX_VAR,i,j,k)= 0.d0</i><br><i> solnData(VELY_VAR,i,j,k)= 0.d0</i><br><i> solnData(VELZ_VAR,i,j,k)= 0.d0</i><br><br><i> </i><br><i> </i><br><br><i> ! Compute the gas energy and set the gamma-values needed for the EOS</i><br>
<i> ekinZone = 0.5 * dot_product(solnData(VELX_VAR:VELZ_VAR,i,j,k),&</i><br><i> solnData(VELX_VAR:VELZ_VAR,i,j,k))</i><br><br><i> ! specific internal energy</i><br>
<i> eintZone = solnData(PRES_VAR,i,j,k)/(sim_gamma-1.)/solnData(DENS_VAR,i,j,k)</i><br><br><i> ! total specific gas energy</i><br><i> enerZone = eintZone + ekinZone</i><br><br><i> ! Take a limit value</i><br>
<i> enerZone = max(enerZone, sim_smallP)</i><br><br><i> solnData(ENER_VAR,i,j,k)=enerZone</i><br><i> solnData(EINT_VAR,i,j,k)=eintZone</i><br><i> solnData(GAMC_VAR,i,j,k)=sim_gamma</i><br>
<i> solnData(GAME_VAR,i,j,k)=sim_gamma</i><br><br><i></i><br><br><i> enddo</i><br><i> enddo</i><br><i> enddo</i><br><br><br><i>!! CURL:</i><br><i> do k = blkLimitsGC(LOW,KAXIS),blkLimitsGC(HIGH,KAXIS)+1</i><br>
<i> do j = blkLimitsGC(LOW,JAXIS),blkLimitsGC(HIGH,JAXIS)+1</i><br><i> do i = blkLimitsGC(LOW,IAXIS),blkLimitsGC(HIGH,IAXIS)+1</i><br><br><i></i><br><i> if ( (j <=blkLimitsGC(HIGH,JAXIS)).and.&</i><br>
<i> (k <=blkLimitsGC(HIGH,KAXIS)) ) then</i><br><i> facexData(MAG_FACE_VAR,i,j,k)= -(Ay(i ,j ,k+1)-Ay(i,j,k))/dz &</i><br><i>+(Az(i ,j+1,k )-Az(i,j,k))/dy</i><br><i> end if</i><br>
<br><i> if ( (i <=blkLimitsGC(HIGH,IAXIS)).and.&</i><br><i> (k <=blkLimitsGC(HIGH,KAXIS)) ) then</i><br><i> faceyData(MAG_FACE_VAR,i,j,k)= (Ax(i ,j ,k+1)-Ax(i,j,k))/dz &</i><br>
<i>-(Az(i+1,j ,k )-Az(i,j,k))/dx</i><br><i> end if</i><br><br><i> if ( (j <=blkLimitsGC(HIGH,JAXIS)).and.&</i><br><i> (i <=blkLimitsGC(HIGH,IAXIS)) ) then</i><br><i> facezData(MAG_FACE_VAR,i,j,k)= -(Ax(i ,j+1,k )-Ax(i,j,k))/dy &</i><br>
<i>+(Ay(i+1,j ,k )-Ay(i,j,k))/dx</i><br><i> end if</i><br><i>!N end</i><br><br><i> enddo</i><br><i> enddo</i><br><i> enddo</i><br><br><i> </i>!!B, DIVB AND MAGP:<br><i> do k=blkLimitsGC(LOW,KAXIS),blkLimitsGC(HIGH,KAXIS)</i><br>
<i> do j = blkLimitsGC(LOW,JAXIS),blkLimitsGC(HIGH,JAXIS)</i><br><i> do i = blkLimitsGC(LOW,IAXIS),blkLimitsGC(HIGH,IAXIS)<br><br></i><br><i> solnData(MAGX_VAR,i,j,k) = 0.5*(facexData(MAG_FACE_VAR,i,j,k)+facexData(MAG_FACE_VAR,i+1,j,k))</i><br>
<i> solnData(MAGY_VAR,i,j,k) = 0.5*(faceyData(MAG_FACE_VAR,i,j,k)+faceyData(MAG_FACE_VAR,i,j+1,k))</i><br><i> if (NDIM == 3) then</i><br><i> solnData(MAGZ_VAR,i,j,k) = 0.5*(facezData(MAG_FACE_VAR,i,j,k)+facezData(MAG_FACE_VAR,i,j,k+1))</i><br>
<i> endif</i><br><br><br><i> </i><br><i></i><div style="margin-left: 40px;"><i> solnData(DIVB_VAR,i,j,k)= &</i><br></div><i> (facexData(MAG_FACE_VAR,i+1,j, k ) - facexData(MAG_FACE_VAR,i,j,k))/dx &</i><br>
<i> + (faceyData(MAG_FACE_VAR,i, j+1,k ) - faceyData(MAG_FACE_VAR,i,j,k))/dy</i><br><br><i> solnData(DIVB_VAR,i,j,k)= solnData(DIVB_VAR,i,j,k) &</i><br><i> + (facezData(MAG_FACE_VAR,i, j, k+1) - facezData(MAG_FACE_VAR,i,j,k))/dz</i><br>
<i></i><br><i></i><br><i> solnData(MAGP_VAR,i,j,k) = .5*dot_product(solnData(MAGX_VAR:MAGZ_VAR,i,j,k),&</i><br><i> solnData(MAGX_VAR:MAGZ_VAR,i,j,k))</i><br>
<br><br><i> enddo</i><br><i> enddo</i><br><i> enddo</i><br><br><br><br><i> ! Release pointer</i><br><i> call Grid_releaseBlkPtr(blockID,solnData,CENTER)</i><br><br><i>#if NFACE_VARS > 0</i><br><i> if (sim_killdivb) then</i><br>
<i> call Grid_releaseBlkPtr(blockID,facexData,FACEX)</i><br><i> call Grid_releaseBlkPtr(blockID,faceyData,FACEY)</i><br><i> if (NDIM == 3) call Grid_releaseBlkPtr(blockID,facezData,FACEZ)</i><br><i> endif</i><br>
<i>#endif</i><br><br><i>deallocate(xCoord)</i><br><i>deallocate(xCoordL)</i><br><i>deallocate(xCoordR)</i><br><i>!</i><br><i>deallocate(yCoord)</i><br><i>deallocate(yCoordL)</i><br><i>deallocate(yCoordR)</i><br><br><i>deallocate(zCoord)</i><br>
<i>deallocate(zCoordL)</i><br><i>deallocate(zCoordR)</i><br><br><i>#ifndef FIXEDBLOCKSIZE</i><br><i> deallocate(Az)</i><br><i> deallocate(Ax)</i><br><i> deallocate(Ay)</i><br><i></i><i>#endif</i><br><br><i>end subroutine Simulation_initBlock</i><br>
</div><br>---------------------------------<br>flash.par:<br><br><br><div style="margin-left: 40px;">sim_n_cube = 73 #= 64+8guard+1curl<br>nblockx = 1 <br>nblocky = 1 <br>nblockz = 1 <br>lrefine_max = 5 <br>lrefine_min = 3 <br>
refine_var_1 = "velx"<br>nrefs = 2<br><br>run_comment = "14Feb09-1"<br>log_file = "14Feb09-1.log"<br>basenm = "14Feb09-1-"<br>restart = .false.<br>
nend = 10<br>checkPointFileNumber = 1<br>tmax = 4.0 <br>checkpointFileIntervalTime = 1.0 <br>plotFileIntervalTime = 8.e-4 #2.5e-2 <br>dtmax = 0.02 #0.01<br>plotFileNumber = 0<br>
<br>sim_beta = 100.0 #Magnetic beta<br><br>gamma = 1.666666666667<br><br>geometry = "cartesian"<br><br>xmin = -0.5<br>xmax = 0.5<br>ymin = -0.5 <br>
ymax = 0.5<br>zmin = -0.5<br>zmax = 0.5<br>xl_boundary_type = "outflow" <br>xr_boundary_type = "outflow"<br>yl_boundary_type = "outflow"<br>yr_boundary_type = "outflow"<br>
zl_boundary_type = "outflow"<br>zr_boundary_type = "outflow"<br><br># Simulation (grid, time, I/O) parameters<br><br>cfl = 0.3<br>dtini = 1.e-12<br>plot_var_1 = "dens"<br>
plot_var_2 = "velx"<br>plot_var_3 = "vely"<br>plot_var_4 = "velz"<br>plot_var_5 = "pres"<br>plot_var_6 = "magx"<br>plot_var_7 = "magy"<br>
plot_var_8 = "magz"<br>plot_var_9 = "divb"<br>plot_var_10 = "magp"<br><br>convertToConsvdInMeshInterp = .true.<br><br>eintSwitch = 1.e-6<br><br>#MHD:<br>ForceHydroLimit = .false. #(set .true. when using hydro limit B=0)<br>
sim_B_cluster = 1.0<br>UnitSystem = "none"<br>resistive_mhd = .false.<br>killdivb = .true. <br><br>flux_correct = .true.<br><br>order = 2 #(first / second order scheme)<br>slopeLimiter = "vanLeer" #(minmod, mc, vanLeer, hybrid, limited)<br>
LimitedSlopeBeta= 1. # only needed for the limited slope by Toro<br>charLimiting = .true. #(.false. will give primitive limitings)<br>E_modification = .true. #(.false. will use simple arithmetic avg)<br>energyFix = .true. #(.false. will not fix total energy due to div-free magnetic pressure)<br>
facevar2ndOrder = .true. #(.false. will give less accurate solution but fast performance)<br><br>RiemannSolver = "hll" <br> <br># CTU integrator<br>CTU = .false. #(6-solve CTU scheme, not recommended)<br>
<br># Prolongation method of Facevars<br>prolMethod = "injection_prol" #(injecton_prol, balsara_prol)<br><br>#ICM initial conditions with a central initial jet /10:<br>smlrho = 1.0e-4 # dens_j/10.<br>
smallp = 1.e-3 # ~ min(pres)_ICM/10<br>smalle = 0.00900001 #(1/gamma)/((gamma-1)*1)/10: initial IMC Therm Energy/10<br> #15.0<br>#smallt = 1.000<br><br>#Total energy minimum limit<br>
smallu = 0.000900001 #[(1/gamma)/((gamma-1)*1) + 0.]/10 initial IMC Total Energy/10 <br> #15.5 #smallp/((gamma-1)smlrho)+.5(1) <br> #1815.0 #smallp/((gamma-1)smlrho)+.5(vel_j)^2<br>
<br>tiny = 1.e-12<br></div><br><br>Best,<br><br>Martin HE.<br><a href="mailto:mh475@cam.ac.uk">mh475@cam.ac.uk</a><br>