[FLASH-USERS] hydro/MHD: clarification on requirements for initial and boundary conditions

[Jonathan Thurgood]

Dear All,

I'm a new FLASH user and I'm setting up my first problem. I'm looking for clarification that I'm correct on a few issues before pressing ahead with more complicated setups because I worry about accidently introducing a problem that will sabotage me down the line.  Any comments from experienced users would be much appreciated.

Question 1:  Which variables must be specified at a minimum for a self-consistent MHD initial condition and boundary condition?

I think it is the case that solution variables in Table 14.2 (dens, velx, vely, velz, pres, ener, temp) and Table 14.7 (magx, magy, magz,  magp, divb) must be specified in the initial condition in a self-consistent manner, and any special boundary conditions must correctly account for these variables.

Am I correct in thinking that the above *must* be specified, and that the only consequence of not setting additional variables like current at t=0 is that they will just give you the "wrong" output at t=0, but then will be correctly derived from the solution variables for t>0 output?

Question 2: What should ENER_VAR be given as?

If using one of the MHD solvers, should the user set in Simulation_initBlock.F90 (and in boundary conditions if using a specialised BC):

(1) ENER_VAR = ekin + eint

Or,

(2) ENER_VAR = ekin + eint + emag

The MHD examples provided of course all seem to follow (1) during initialisation, but this message http://flash.uchicago.edu/pipermail/flash-users/2012-April/001079.html threw me. Is the correct way to interpret  this that any manipulation of ENER_VAR outside of hy_* files such as the initial conditions and specialist boundary conditions should follow (1)?

Question 3: Correct logic for specifying spatially dependent, face-centred variables in initial conditions?

I assume that if you fill local variables with a loop over block indices that you run the risk of not filling the 'leftmost' or rightmost cell faced-variables and have to account for it. Can someone please see if this looks ok?

!loop is over cell indices, i.e. for face-centred values if not careful might miss leftmost or rightmost face

  do k = blkLimitsGC(LOW,KAXIS),blkLimitsGC(HIGH,KAXIS)
  do j = blkLimitsGC(LOW,JAXIS),blkLimitsGC(HIGH,JAXIS)
  do i = blkLimitsGC(LOW,IAXIS),blkLimitsGC(HIGH,IAXIS)

    xx = xCoord(i) !cell-center position
    yy = yCoord(j)

#if NFACE_VARS > 0
    if (sim_killdivb) then

    xxL = xCoordL(i)  !left face
    xxR = xCoordR(i)  !right face
    yyL = yCoordL(j)
    yyR = yCoordR(j)

    endif
#endif

    !calc. local values (<var>Zone) for the cell centered
    !variables

    !the following  are user-set

    eintZone = 1e-6  !set eint as small for zero beta?
    densZone = 1.0

    velxZone = 0.0
    velyZone = 0.0
    velzZone = 0.0

    magxZone = xx      !don't forget - any changes to the analytical form (e.g, if rotate B=[x,-y] -> B=[y,x])
    magyZone = -yy     !must be appropriately replicated on later lines
    magzZone = 0.0      !for cell-face values


    !the following are derived from the above

    presZone = enerZone * (sim_gamma-1.) * densZone

    ekinZone = 0.5 * &
      (velxZone**2. + velyZone**2. + velzZone**2.)

    enerZone = eintZone + ekinZone

    !asign to cell-centered data block

    solnData(VELX_VAR,i,j,k) = velxZone
    solnData(VELY_VAR,i,j,k) = velyZone
    solnData(VELZ_VAR,i,j,k) = velzZone

    solnData(DENS_VAR,i,j,k) = densZone
    solnData(PRES_VAR,i,j,k) = presZone
!   solnData(TEMP_VAR,i,j,k) =  !required?
    solnData(ENER_VAR,i,j,k) = enerZone
    solnData(EINT_VAR,i,j,k) = eintZone

    solnData(MAGX_VAR,i,j,k) = magxZone
    solnData(MAGY_VAR,i,j,k) = magyZone
    solnData(MAGZ_VAR,i,j,k) = magzZone

    solnData(MAGP_VAR,i,j,k) = 0.5 * dot_product( &
                                  solnData(MAGX_VAR:MAGZ_VAR,i,j,k),&
                                  solnData(MAGX_VAR:MAGZ_VAR,i,j,k))
    solnData(DIVB_VAR,i,j,k) = 0.0

    solnData(GAMC_VAR,i,j,k) = sim_gamma
    solnData(GAME_VAR,i,j,k) = sim_gamma


    !now do cell face values if needed for USM
#if NFACE_VARS > 0

    if (sim_killdivb) then

    magxZone = xxL
    magyZone = -yyL  !fill from left, account for rightmost later
    magzZone = 0.0

    facexData(MAG_FACE_VAR,i,j,k) =  magxZone
    faceyData(MAG_FACE_VAR,i,j,k) =  magyZone

    ! extra specification for "rightmost" cell-face if at RHS of a block

    if (i == blkLimitsGC(HIGH,IAXIS) ) then
    magxZone = xxR
    facexData(MAG_FACE_VAR,i+1,j,k) =  magxZone
    endif
    if (j == blkLimitsGC(HIGH,JAXIS)  ) then
    magyZone = -yyR
    facexData(MAG_FACE_VAR,i,j+1,k) =  magyZone
    endif

    !re-specify cell-centred values by averaging face-values across faces.
    !this should happen after filling "extra" face cell (here rightmost)
    !so that i+1 isn't missed out if i == imax

    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))
    solnData(MAGP_VAR,i,j,k) = 0.5 * dot_product( &
                                  solnData(MAGX_VAR:MAGZ_VAR,i,j,k),&
                                  solnData(MAGX_VAR:MAGZ_VAR,i,j,k))

    endif
#endif


  enddo
  enddo
  enddo


Question 4:

In my initial condition I'm setting up a spatially dependent field from an analytical prescription (a 2D null point, B=[x,-y]). I would like to use a special boundary condition whereby the boundary values are held fixed at their t=0 position.

It seems that my way forward is that I must modify Grid_bcApplyToRegion.F90 and, specifically for the B-field (a function of space), that I must use gr-extendGetCellCoords to get the x and y position of a cell to correctly specify spatially dependent values.

Is there a more simple way to realise this, or must I proceed as above? E.g., is there a way I can simply set the boundary info to be that from t=0?

---

Thank you in advance for taking the time to have a look at this.

Regards,

Jonathan

-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://flash.rochester.edu/pipermail/flash-users/attachments/20160418/829a4dc1/attachment.htm>