!****************************************************************************** ! ! Subroutine: init_block() ! ! Description: Initializes fluid data for a specified block. ! This version sets up the magnetic rotator MHD problem. ! ! References: Mouschovias, T. Ch., & Paleologou, E. V. 1980, ApJ, 237, 877 ! Stone, J. M., & Norman, M. L. 1992, ApJS, 80, 791 ! Stone, J. M., Hawley, J. F., Evans, C. R., & Norman, M. L. ! 1992 , ApJ, 388, 415 ! ! Parameters: block_no The number of the block to initialize ! ! subroutine init_block (block_no) ! !============================================================================== use logfile, ONLY: write_logfile use physical_constants use runtime_parameters, ONLY : get_parm_from_context use dBase, ONLY: iu_bnd, ju_bnd, ku_bnd, & & ndim, & & nguard, ionmax, & & dBaseKeyNumber, dBaseSpecies, & & dBaseGetCoords, dBasePutData, & & dBasePropertyInteger implicit none integer, INTENT(in) :: block_no integer, save :: MyPE, MasterPE integer, save :: idens, ipres, iener, igamc, igame integer, save :: ivelx, ively, ivelz integer, save :: imagx, imagy, imagz integer, save :: iPoint integer, save :: iXcoord, iYcoord, iZcoord integer, save :: izn, iznl, iznr real, save :: const_PI real, save :: gamma real, save :: rho_disk, rho_amb, B_z, z_disk, omega_disk, const_press real, save :: r_disk real :: rho, pres, eng, gamc, game real :: vx, vy, vz, bx, by, bz integer, parameter :: q = max(iu_bnd,ju_bnd,ku_bnd) real, DIMENSION(q) :: x, y, z, zL, zR real :: xx, yy, zz, zzL, zzR, cyl_r, theta integer :: i, j, k, n logical, save:: first_call = .true. !============================================================================== !------------------------------------------------------------------------------ ! Compute direction cosines for the shock normal. if (first_call) then MyPE = dBasePropertyInteger('MyProcessor') MasterPE = dBasePropertyInteger('MasterProcessor') call get_constant_from_db("pi", const_PI) idens = dBaseKeyNumber('dens') ivelx = dBaseKeyNumber('velx') ively = dBaseKeyNumber('vely') ivelz = dBaseKeyNumber('velz') imagx = dBaseKeyNumber('magx') imagy = dBaseKeyNumber('magy') imagz = dBaseKeyNumber('magz') iener = dBaseKeyNumber('ener') ipres = dBaseKeyNumber('pres') igamc = dBaseKeyNumber('gamc') igame = dBaseKeyNumber('game') iPoint = dBaseKeyNumber('Point') iXcoord = dBaseKeyNumber('xCoord') iYcoord = dBaseKeyNumber('yCoord') iZcoord = dBaseKeyNumber('zCoord') izn = dBaseKeyNumber('zn' ) ! iznl = dBaseKeyNumber('znl' ) ! iznr = dBaseKeyNumber('znr' ) call get_parm_from_context('rho_disk', rho_disk) call get_parm_from_context('rho_amb', rho_amb) call get_parm_from_context('B_z', B_z) call get_parm_from_context('z_disk', z_disk) call get_parm_from_context('r_disk', r_disk) call get_parm_from_context('omega_disk', omega_disk) call get_parm_from_context('const_press', const_press) call get_parm_from_context('gamma', gamma) !------------------------------------------------------------------------------ ! Write a message to stdout describing the problem setup. if (MyPE .eq. MasterPE) then call write_logfile("flash: initializing the magnetic rotator problem.") write (*,*) write (*,*) 'Flash: initializing the magnetic rotator problem.' write (*,*) write (*,1) 'rho_disk = ', rho_disk, 'rho_amb = ', rho_amb write (*,1) 'B_z = ', B_z, 'z_disk = ', z_disk, & & 'omega_disk = ', omega_disk write (*,2) 'gamma = ', gamma, 'ndim = ', ndim write (*,*) 1 format (1X, 1P, 5(A7, E13.6, :, 1X)) 2 format (1X, 1P, 2(A7, E13.6, 1X), A7, I13) endif first_call = .false. endif !------------------------------------------------------------------------------ ! Loop through the block and initialize all values. x(:) = 0.0 y(:) = 0.0 z(:) = 0.0 ! zL(:) = 0.0 ! zR(:) = 0.0 ! If simulation is 3D, initialze, with rotation axis along z. if (ndim == 3) then !***** Get coordinates. I think this gets the center cell values. call dBaseGetCoords(izn, iXcoord, block_no, x) call dBaseGetCoords(izn, iYcoord, block_no, y) call dBaseGetCoords(izn, iZcoord, block_no, z) ! call dBaseGetCoords(iznl, iZcoord, block_no, zL) ! call dBaseGetCoords(iznr, iZcoord, block_no, zR) !***** Loop through the block. do k = 1, ku_bnd zz = z(k) ! zzL= zL(k) ! zzR= zR(k) do j = 1, ju_bnd yy = y(j) do i = 1, iu_bnd xx = x (i) !***** Calculate the cylindrical r and azimuthal angle. cyl_r = sqrt(xx*xx + yy*yy) if (cyl_r > 0.0) then theta = ACOS(xx/cyl_r) if (yy < 0.0) theta = const_PI*2.0 - theta else theta = 0.0 endif !***** Set initial values. game = gamma gamc = gamma rho = rho_amb pres = const_press vx = 0.0 vy = 0.0 vz = 0.0 bx = 0.0 by = 0.0 bz = B_z !***** Create the rotating disk. if ((ABS(zz) < z_disk) .and. (cyl_r < r_disk)) then rho = rho_disk vx =-omega_disk * cyl_r * SIN(theta) vy = omega_disk * cyl_r * COS(theta) ! elseif ((((ABS(zzL) < z_disk) .and. (ABS(zzR) > z_disk)) .or. & ! & ((ABS(zzL) > z_disk) .and. (ABS(zzR) < z_disk))) & ! & .and. (cyl_r < r_disk)) then ! ! rho = 0.25 * rho_disk endif eng = 0.5 * (vx*vx + vy*vy + vz*vz) + & & pres / (gamma - 1.0) / rho call dBasePutData(idens, iPoint, i, j, k, block_no, rho) call dBasePutData(ipres, iPoint, i, j, k, block_no, pres) call dBasePutData(iener, iPoint, i, j, k, block_no, eng) call dBasePutData(ivelx, iPoint, i, j, k, block_no, vx ) call dBasePutData(ively, iPoint, i, j, k, block_no, vy ) call dBasePutData(ivelz, iPoint, i, j, k, block_no, vz ) call dBasePutData(imagx, iPoint, i, j, k, block_no, bx ) call dBasePutData(imagy, iPoint, i, j, k, block_no, by ) call dBasePutData(imagz, iPoint, i, j, k, block_no, bz ) call dBasePutData(igamc, iPoint, i, j, k, block_no, gamc) call dBasePutData(igame, iPoint, i, j, k, block_no, game) do n=1,ionmax call dBasePutData(dBaseSpecies(n),iPoint,i,j,k,block_no, 1.0) enddo enddo enddo enddo ! If simulation is 1D, print error message. elseif (ndim < 2) then if (MyPE .eq. MasterPE) then write(*,*) 'Simulation must have at least 2 dimensions!' endif ! If simulation is 2D, initialize axisymmetric setup, using x,y as r,z. !***** AXISYMMETRY DOES NOT WORK !!!! ***** elseif (ndim == 2) then if (MyPE .eq. MasterPE) then write(*,*) 'Simulation does not work in 2 dimensions!' endif !***** Get coordinates. I think (hope) this gets the center cell !***** values. call dBaseGetCoords(izn, iXcoord, block_no, x) call dBaseGetCoords(izn, iYcoord, block_no, y) !***** Loop through the block. do k = 1, ku_bnd do j = 1, ju_bnd yy = y(j) do i = 1, iu_bnd xx = x (i) !***** Set initial values. game = gamma gamc = gamma rho = rho_amb pres = const_press vx = 0.0 vy = 0.0 vz = 0.0 bx = 0.0 by = B_z bz = 0.0 !***** Create the disk. if (yy .lt. z_disk) then rho = rho_disk vz = omega_disk * xx endif eng = 0.5 * (vx*vx + vy*vy + vz*vz) + & & pres / (gamma - 1.0) / rho call dBasePutData(idens, iPoint, i, j, k, block_no, rho) call dBasePutData(ipres, iPoint, i, j, k, block_no, pres) call dBasePutData(iener, iPoint, i, j, k, block_no, eng) call dBasePutData(ivelx, iPoint, i, j, k, block_no, vx ) call dBasePutData(ively, iPoint, i, j, k, block_no, vy ) call dBasePutData(ivelz, iPoint, i, j, k, block_no, vz ) call dBasePutData(imagx, iPoint, i, j, k, block_no, bx ) call dBasePutData(imagy, iPoint, i, j, k, block_no, by ) call dBasePutData(imagz, iPoint, i, j, k, block_no, bz ) call dBasePutData(igamc, iPoint, i, j, k, block_no, gamc) call dBasePutData(igame, iPoint, i, j, k, block_no, game) do n=1,ionmax call dBasePutData(dBaseSpecies(n),iPoint,i,j,k,block_no, 1.0) enddo enddo enddo enddo else if (MyPE .eq. MasterPE) then write(*,*) 'Simulations must have 1, 2, or 3 dimensions!' endif endif !============================================================================== end subroutine init_block