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