[FLASH-USERS] USM Custom BC / Grid_bcApplyToRegionSpecialized.F90 / Spurious Current at boundaries for magnetic null point
Jonathan Thurgood
jonathan.thurgood at northumbria.ac.uk
Wed May 4 08:15:39 EDT 2016
Dear Flash Dev's and Users,
I've been setting up a USM simulation where the equilibrium is a linear 2D magnetic X-point of the form B=[x,-y,0]) (e.g., http://ukads.nottingham.ac.uk/abs/2011A%26A...533A..18M , http://ukads.nottingham.ac.uk/abs/2004A%26A...420.1129M).
The default boundary options don't play well with the null due to the non-constant / curved field through the boundary. The various specifications end up implying some spurious current on the boundaries, and so launch waves into the domain. To overcome this I need to implement my own boundary conditions in Grid_bcApplyToRegionSpecialised.
One particular approach I have tried is to force the B-field variables to be their equilibrium values on the boundary (from the analytical prescription), which I believe should properly maintain the equilibrium. Specifically, this boundary condition is intended to set velocity to zero, other UNK variables similar to "reflecting" (without the sign change) and the B field to B=[x,-y,0].
This doesn't seem to be working as I expected - after only coding up the custom BC for the lower face (but not the high faces), I found that that the result is that the current along each boundary is the same at a given output (as opposed to the lower left corner being different to the upper right, say). This suggests that my "implementation" has some silly mistake and key lines are not being executed when I think they should be.
I further tested this by going to the lines I had expected to modify the B-field variables and setting the variable to be a large value (1e6) to see if it changes (breaks) the simulation. Many of the lines actually appear to do nothing / not be executed and I am not sure why, likely I have not properly understood the logic used in their execution. I have tagged these "!*1e6 does nothing" or "!* 1e6 breaks simulation" below.
If anyone could give me some pointers with this I would be extremely grateful.
Best regards,
Jonathan
(Post-doc working in magnetic reconnection, Northumbria University Solar Physics Group, UK)
subroutine Grid_bcApplyToRegionSpecialized(bcType,gridDataStruct,&
guard,axis,face,regionData,regionSize,mask,applied,&
blockHandle,secondDir,thirdDir,endPoints,blkLimitsGC, idest)
use Grid_data, ONLY : gr_meshMe
#include "constants.h"
implicit none
integer, intent(IN) :: bcType,axis,face,guard,gridDataStruct
integer,dimension(REGION_DIM),intent(IN) :: regionSize
real,dimension(regionSize(BC_DIR),&
regionSize(SECOND_DIR),&
regionSize(THIRD_DIR),&
regionSize(STRUCTSIZE)),intent(INOUT)::regionData
logical,intent(IN),dimension(regionSize(STRUCTSIZE)):: mask
logical, intent(OUT) :: applied
integer,intent(IN) :: blockHandle
integer,intent(IN) :: secondDir,thirdDir
integer,intent(IN),dimension(LOW:HIGH,MDIM) :: endPoints, blkLimitsGC
integer,intent(IN),OPTIONAL:: idest
!!!jonathan's added declarations
integer :: i,j, k,ivar,je,ke,n,varCount,bcTypeActual
logical :: isFace
integer :: sizeX, sizeY, sizeZ
real, allocatable,dimension(:) :: xCoord, yCoord, XCoordF, yCoordF
! stub provided
applied = .false.
!added lines to the stub below this comment
select case (bcType)
case(USER_DEFINED)
applied = .true.
case default
applied = .false.
return
end select
je=regionSize(SECOND_DIR)
ke=regionSize(THIRD_DIR)
varCount=regionSize(STRUCTSIZE)
isFace = (gridDataStruct==FACEX).and.(axis==IAXIS)
isFace = isFace.or.((gridDataStruct==FACEY).and.(axis==JAXIS))
isFace = isFace.or.((gridDataStruct==FACEZ).and.(axis==KAXIS))
!some allocation and calls to populate coordinate arrays to allow analytical prescription of B on boundary
sizeX = blkLimitsGC(HIGH,IAXIS)-blkLimitsGC(LOW,IAXIS)+1 !double check sizes
sizeY = blkLimitsGC(HIGH,JAXIS)-blkLimitsGC(LOW,JAXIS)+1
sizeZ = blkLimitsGC(HIGH,KAXIS)-blkLimitsGC(LOW,KAXIS)+1
allocate(xCoord(sizeX) )
allocate(yCoord(sizeY) )
allocate(xCoordF(sizeX+1) )
allocate(yCoordF(sizeY+1) )
xCoord = 0.0
yCoord = 0.0
xCoordF = 0.0
yCoordF = 0.0
call gr_extendedGetCellCoords(IAXIS, blockHandle, gr_meshMe, CENTER , .true. , xCoord, sizeX)
call gr_extendedGetCellCoords(JAXIS, blockHandle, gr_meshMe, CENTER , .true. , yCoord, sizeY)
call gr_extendedGetCellCoords(IAXIS, blockHandle, gr_meshMe, FACES, .true. , xCoordF, sizeX+1)
call gr_extendedGetCellCoords(JAXIS, blockHandle, gr_meshMe, FACES, .true. , yCoordF, sizeY+1)
!main loop over variables
do ivar = 1,varCount
if(mask(ivar)) then
call gr_bcMapBcType(bcTypeActual,bcType,ivar,gridDataStruct,axis,face,idest)
if(face==LOW) then
select case (bcTypeActual)
case(USER_DEFINED)
k = 2*guard+1
if(isFace)k=k+1
do i = 1,guard
regionData(i,1:je,1:ke,ivar) = regionData(k-i,1:je,1:ke,ivar) ! do znd-type for all variables first, then overwrite
! e.g. b field with more specific values
!*1e6 breaks simulation
#ifdef VELX_VAR
if ( ivar==VELX_VAR ) regionData(i,1:je,1:ke,ivar) = 0.0 !overwrite previously assigned value with zero for velocities !*1e6 does nothing
#endif
#ifdef VELY_VAR
if ( ivar==VELY_VAR ) regionData(i,1:je,1:ke,ivar) = 0.0 !*1e6 does nothing
#endif
#ifdef VELZ_VAR
if ( ivar==VELZ_VAR ) regionData(i,1:je,1:ke,ivar) = 0.0 !*1e6 does nothing
#endif
!set B=[x,-y,0] on lower face to make current free boundary
#ifdef MAGX_VAR
if( (ivar==MAGX_VAR).AND.(axis=iaxis) ) regionData(i,1:je,1:ke,ivar) = xCoord(i) !bx = x !*1e6 does nothing
if( (ivar==MAGX_VAR).AND.(axis=jaxis) ) then
do j=1,je
regionData(i,j,1:ke,ivar) = xCoord(j) !*1e6 does nothing
enddo
endif
#endif
#ifdef MAGY_VAR
if( (ivar==MAGY_VAR).AND.(axis=jaxis) ) regionData(i,1:je,1:ke,ivar) = -yCoord(i) !by = -y !*1e6 does nothing
if( (ivar==MAGY_VAR).AND.(axis=iaxis) ) then
do j=1,je
regionData(i,j,1:ke,ivar) = -yCoord(j) !*1e6 does nothing
enddo
endif
#endif
#ifdef MAG_FACE_VAR
if( (ivar==MAG_FACE_VAR).AND.(axis=iaxis) ) regionData(i,1:je,1:ke,ivar) = xCoordF(i) !*1e6 does nothing
if( (ivar==MAG_FACE_VAR).AND.(axis=jaxis) ) regionData(i,1:je,1:ke,ivar) = -yCoordF(i) !*1e6 does nothing
#endif
!should also do intermediate step mag field on face if doing a custom BC
!http://flash.uchicago.edu/pipermail/flash-users/2012-February/001021.html
#ifdef MAGI_FACE_VAR
if( (ivar==MAGI_FACE_VAR).AND.(axis=iaxis) ) regionData(i,1:je,1:ke,ivar) = xCoordF(i) !*1e6 does nothing
if( (ivar==MAGI_FACE_VAR).AND.(axis=jaxis) ) regionData(i,1:je,1:ke,ivar) = -yCoordF(i) !*1e6 does nothing
#endif
end do
end select
else !face==HIGH
select case (bcTypeActual)
case(USER_DEFINED)
k = 2*guard+1
if(isFace)k=k+1
do i = 1,guard
regionData(k-i,1:je,1:ke,ivar)= regionData(i,1:je,1:ke,ivar) !znd for all first
#ifdef VELX_VAR
if ( ivar==VELX_VAR ) regionData(k-i,1:je,1:ke,ivar) = 0.0 !over-write with zero for velocities
#endif
#ifdef VELY_VAR
if ( ivar==VELY_VAR ) regionData(k-i,1:je,1:ke,ivar) = 0.0 !over-write with zero for velocities
#endif
#ifdef VELZ_VAR
if ( ivar==VELZ_VAR ) regionData(k-i,1:je,1:ke,ivar) = 0.0 !over-write with zero for velocities
#endif
end do
end select
endif
endif
enddo !end ivar do loop
deallocate(xCoord)
deallocate(yCoord)
deallocate(xCoordF)
deallocate(yCoordF)
return
end subroutine Grid_bcApplyToRegionSpecialized
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://flash.rochester.edu/pipermail/flash-users/attachments/20160504/f8f469d3/attachment.htm>
More information about the flash-users
mailing list