[FLASH-USERS] 3DMHD +usm, jumps of divb at refinement jumps
Martin Huarte-Espinosa
mh475 at cam.ac.uk
Fri Feb 20 04:56:30 EST 2009
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
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://flash.rochester.edu/pipermail/flash-users/attachments/20090220/29f32838/attachment.htm>
More information about the flash-users
mailing list