[FLASH-USERS] User Boundary conditions for USM - WARNING after gc filling: min. unk(DENS_VAR)=-1
Jonathan Thurgood
jonathan.thurgood at northumbria.ac.uk
Fri May 6 08:00:55 EDT 2016
Dear All,
In case anyone ever has the same sort issue in the future I thought I would submit the solution to the mailing list for archiving.
This problem came about from an ambiguity in my given attempt resulting from a failure to explicitly distinguish between gridDataStruct == CENTER or FACEX,FACEY,FACEZ
If you look at Flash.h in the object directory, you see that the cell-centred variables are numbered independently of the face variables (e.g., CURZ_VAR=1, DENS_VAR=2,.. etc, and MAG_FACE_VAR=1,MAGI_FACE_VAR=2). Failing to distinguish between the two allowed the "if ivar = MAGI_FACE_VAR" statement to be evaluated as true for the density also.
If I had looked more carefully at the examples for the outflow/diode/reflecting boundary conditions, I would have noticed that this distinction is actually made when selecting which sign = +-1 is applied. Schoolboy error!
Cheers,
Jonathan
From: flash-users-bounces at flash.uchicago.edu [mailto:flash-users-bounces at flash.uchicago.edu] On Behalf Of Jonathan Thurgood
Sent: 05 May 2016 13:41
To: flash-users at flash.uchicago.edu
Subject: [FLASH-USERS] User Boundary conditions for USM - WARNING after gc filling: min. unk(DENS_VAR)=-1
Dear All,
I spotted some of the sillier errors I sent you all in my message about problems at a 2D null point yesterday (apologies) and got further with the problem.
I now suffer from an error message indicating that the boundary conditions (BC) I attempted somehow forces a -ve density:
(e.g. WARNING after gc filling: min. unk(DENS_VAR)=-1.000000000000000000000 PE=0 block=1 type=1 )
I've been able to reproduce this error in a simplified problem of a 2D homogeneous plasma equilibrium (dens = 1 eint=1e-12 vel=0) with a uniform field of Bx = -1, By = 0.
I attempt to set the B values in the boundaries to their correct value of B=[-1,0] "manually" and get the above error, indicating that I have somehow assigned the value intended for the field to the density.
The relevant bits of my Grid_bcApplyToRegionSpecialized.F90 file are attached below, and commented where appropriate.
I have
(1) set only the magx_var and magy_var (cell-centered) manually whilst changing nothing else - this works
(2) set only the MAG_FACE_VAR values manually - this produces a running code, but introduces an artefact which propagates in from the boundary of the order 1e-7
(3) set only the MAGI_FACE_VAR values manually - this line seems to be responsible for forcing a -ve density In particular. If this value is changed, say from the "real" value of -1 to -2, the error will report that gc filling set DENS_VAR=given negative value there.
Could anyone let me know why (2) introduces the artefact (and if possible how to avoid), and why (3) is apparently setting values for variables other than MAGI_FACE_VAR?
More generally, any extra advice about setting user BC (particularly for USM) would be appreciated also. As far as I know, for a custom BC in USM I need only to set velocities (velx_var, vely_var, velz_var), cell-centred fields (magx_var, magy_var, magz_var), dens, eint (or other depending on EOS) and face-fields (MAG_FACE_VAR and MAGI_FACE_VAR).
Thank you,
Jonathan
!!!!!!!!! Bits of Code Below here !!!!!!!!
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) !all variables are uniform and in equilibrium - this should set
!correct BC values for all variables
!as test, then going to overwrite this to set some variables (namely, fields) manually
!to try and pinpoint the error
!the idea is eventually to set the field in the boundary as f(x,y)
!1. the following two specifications for the cell-centre fields work as expected - no artefacts introduced
!
#ifdef MAGX_VAR
if( ivar==MAGX_VAR ) regionData(i,1:je,1:ke,ivar) = -1.0
#endif
#ifdef MAGY_VAR
if( ivar==MAGY_VAR ) regionData(i,1:je,1:ke,ivar) = 0.0
#endif
!2. The following specifies bx on the i-drn face as -1 and by on the j-drn face as 0.
! the code will run with these lines un-commented but with
! a small propagating artefact is launched from the lower corner
! error in magx in domain is of order 1e-7
#ifdef MAG_FACE_VAR
if( (ivar==MAG_FACE_VAR).AND.(axis==IAXIS) ) regionData(i,1:je,1:ke,ivar) = -1.0
if( (ivar==MAG_FACE_VAR).AND.(axis==JAXIS) ) regionData(i,1:je,1:ke,ivar) = 0.0
#endif
!3. Same thing with MAGI
! This produces a crash as follows
!
! Zero or imaginary sound speed has obtained! Please try other (more diffusive) slope limiter, flux, order, cfl, etc.
! WARNING after gc filling: min. unk(DENS_VAR)=-1.000000000000000000000 PE=0 block=1 type=1
! WARNING after gc filling: min. unk(DENS_VAR)=0.000000000000000000000 PE=4 block=1 type=1
! (etc for other blocks)
!
! Looks as if this has actually set dens_var = -1 or =0 rather than magi_face_var
!
#ifdef MAGI_FACE_VAR
if( (ivar==MAGI_FACE_VAR).AND.(axis==IAXIS) ) regionData(i,1:je,1:ke,ivar) = -1.0
if( (ivar==MAGI_FACE_VAR).AND.(axis==JAXIS) ) regionData(i,1:je,1:ke,ivar) = 0.0
#endif
end do
end select
else !face==HIGH
! (I only tested the above in the lower-facing boundaries)
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)
end do
end select
endif
endif
enddo !end ivar loop
return
end subroutine Grid_bcApplyToRegionSpecialized
From: flash-users-bounces at flash.uchicago.edu<mailto:flash-users-bounces at flash.uchicago.edu> [mailto:flash-users-bounces at flash.uchicago.edu] On Behalf Of Jonathan Thurgood
Sent: 04 May 2016 13:16
To: flash-users at flash.uchicago.edu<mailto:flash-users at flash.uchicago.edu>
Subject: [FLASH-USERS] USM Custom BC / Grid_bcApplyToRegionSpecialized.F90 / Spurious Current at boundaries for magnetic null point
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/20160506/290105ce/attachment.htm>
More information about the flash-users
mailing list