JHUAPL/kaiju
1
0
Fork 0
mirror of https://github.com/JHUAPL/kaiju.git synced 2026-01-09 18:57:53 -05:00
Code Issues Packages Projects Releases 1 Wiki Activity

Merge remote-tracking branch 'origin/raijudev' into raijudev

Browse Source
This commit is contained in:
Anthony M. Sciola
2024-08-22 13:17:06 -04:00
parent bcd00abbc6 61e439c510
commit 373d8289bc
16 changed files with 474 additions and 75 deletions
Download Patch File Download Diff File Expand all files Collapse all files

10
src/base/defs/raijudefs.F90
View File

@@ -95,7 +95,13 @@ module raijudefs
real(rp), parameter :: def_ChorusEMin = 1.1 ! [keV]
! Coupling defaults
real(rp), parameter :: def_vaFracThresh = 0.10 ! min allowable fracton of Alfven to total tube velocity
real(rp), parameter :: def_bminThresh = 1.0 ! [nT] default allowable bmin strencgh
real(rp), parameter :: def_vaFracThresh = 0.10
!! min allowable fracton of Alfven to total tube velocity
real(rp), parameter :: def_PstdThresh = 0.10
!! Max allowable standard deviation of pressure along a field line
real(rp), parameter :: def_normAngle = 25
!! [deg] Max allowable angle between any two normals surrounding a cell corner
real(rp), parameter :: def_bminThresh = 1.0
!! [nT] default allowable bmin strencgh
end module raijudefs

1
src/base/types/ebtypes.F90
View File

@@ -76,6 +76,7 @@ module ebtypes
type ebTrc_T
real(rp) :: OCb !Topology
real(rp) :: dvB,bD,bP,bS !Flux-tube volume, averaged density/pressure, integrated entropy
real(rp) :: stdD, stdP ! Standard deviation of density and pressure
real(rp), dimension(NDIM) :: MagEQ, xEPm,xEPp !xyz of equator/ -/+ field endpoints
real(rp) :: bMin !Minimum B (@ equator)
end type ebTrc_T

8
src/base/types/raijuTypes.F90
View File

@@ -244,7 +244,7 @@ module raijutypes
integer :: nFluidIn = 0
type(mhd2raiSpcMap_T), dimension(:), allocatable :: fluidInMaps
! Coupling-related knobs
real(rp) :: vaFracThresh, bminThresh
real(rp) :: vaFracThresh, bminThresh, normAngThresh, PstdThresh
character(len=strLen) :: icStr
procedure(raijuStateIC_T ), pointer, nopass :: initState => NULL()
@@ -360,6 +360,10 @@ module raijutypes
!! (Ngi, Ngj, Ns) [nPa] Average pressure along flux tube
real(rp), dimension(:,:,:), allocatable :: Davg
!! (Ngi, Ngj, Ns) [#/cc] Average density along flux tube
real(rp), dimension(:,:,:), allocatable :: Pstd
!! (Ngi, Ngj, Ns) Normalized standard deviation of the species pressure along the field line
real(rp), dimension(:,:,:), allocatable :: Dstd
!! (Ngi, Ngj, Ns) Normalized standard deviation of the species density along the field line
real(rp), dimension(:,:,:), allocatable :: Bmin
!! (Ngi+1, Ngj+1, NDIM) [nT] Bmin vector
real(rp), dimension(:,:,:), allocatable :: xyzMin
@@ -384,6 +388,8 @@ module raijutypes
real(rp), dimension(:,:), allocatable :: espot
!! (Ngi+1, Ngj+1) [kV] electro-static potential
integer, dimension(:), allocatable :: bndLoc
!! (Ngi) i value of boundary between inactive and buffer domain
! (Ngi, Ngj) cell-centered values
integer , dimension(:,:), allocatable :: active
!! (Ngi, Ngj) (-1=inactive, 0=buffer, 1=active)

2
src/base/types/voltCplTypes.F90
View File

@@ -16,6 +16,8 @@ module voltCplTypes
!! Average plasma beta
real(rp), dimension(0:MAXTUBEFLUIDS) :: Pave, Nave
!! average pressure, average density
real(rp), dimension(0:MAXTUBEFLUIDS) :: Pstd, Nstd
!! standard deviation of pressure and densiy along the field line
real(rp) :: X_bmin(NDIM)
integer(ip) :: topo
real(rp) :: latc,lonc !Conjugate lat/lon

2
src/chimp/sliceio.F90
View File

@@ -380,6 +380,8 @@ module sliceio
call AddOutVar(IOVars,"bP" ,ebTrcIJ(:,:)%bP )
call AddOutVar(IOVars,"bS" ,ebTrcIJ(:,:)%bS )
call AddOutVar(IOVars,"bMin",ebTrcIJ(:,:)%bMin)
call AddOutVar(IOVars,"stdD",ebTrcIJ(:,:)%stdD)
call AddOutVar(IOVars,"stdP",ebTrcIJ(:,:)%stdP)
!Equator and end-points
call AddOutVar(IOVars,"xBEQ",ebTrcIJ(:,:)%MagEQ(XDIR))

97
src/chimp/streamline.F90
View File

@@ -155,6 +155,8 @@ module streamline
ebTrc%bP = 0.0
ebTrc%bS = 0.0
ebTrc%bMin = 0.0
ebTrc%stdP = 0.0
ebTrc%stdD = 0.0
ebTrc%MagEQ(:) = 0.0
ebTrc%xEPm (:) = 0.0
@@ -169,6 +171,8 @@ module streamline
if (ebTrc%OCb > 0) then
!Get flux-tube integrals
call FLThermo(Model,ebState%ebGr,bTrc,ebTrc%bD,ebTrc%bP,ebTrc%dvB)
call FLStdev(Model, ebState%ebGr,bTrc,ebTrc%bD,ebTrc%bP,ebTrc%stdD,ebTrc%stdP)
!call FLStdev(Model, ebState%ebGr,bTrc,ebTrc%stdD,ebTrc%stdP)
ebTrc%bS = FLEntropy(Model,ebState%ebGr,bTrc)
!Get magnetic equator info
@@ -272,7 +276,7 @@ module streamline
integer , intent(in ), optional :: sOpt
integer :: k,s0
real(rp) :: bMag,dl,eP,eD,ePb !Edge-centered values
real(rp) :: dvB_active,bMag,dl,eP,eD,ePb !Edge-centered values
real(rp) :: bPb,bBeta
!Zero out accumulators
@@ -281,6 +285,7 @@ module streamline
dvB = 0.0
bPb = 0.0
bBeta = 0.0
dvB_active = 0.0
if (.not. bTrc%isGood) return
@@ -309,20 +314,26 @@ module streamline
ePb = 1.0e+14*(bMag*oBScl*1.0e-9/0.501)**2.0 !Edge mag pressure in nPa
!Now accumulate into flux-tube integrals
dvB = dvB + dl/bMag
bD = bD + eD*dl/bMag
bP = bP + eP*dl/bMag
bPb = bPb + ePb*dl/bMag
bBeta = bBeta + (eP/ePb)*dl/bMag
! Always accumulate total flux tube volume
dvB = dvB + dl/bMag
! Only accumulate other quantities if in active domain, se we have more sensible flux tube averages
if (bTrc%ijk(k,XDIR) > ebGr%is+4 .and. bTrc%ijk(k+1,XDIR) > ebGr%is+4) then
bD = bD + eD*dl/bMag
bP = bP + eP*dl/bMag
bPb = bPb + ePb*dl/bMag
bBeta = bBeta + (eP/ePb)*dl/bMag
dvB_active = dvB_active + dl/bMag
endif
enddo
dvB_active = max(dvB_active, TINY) ! Max just do prevent nans for lines that were completely under the gap+padding
!Now turn flux-tube integrals of quantities into flux-tube averages
bD = bD/dvB
bP = bP/dvB
bPb = bPb/dvB
bD = bD /dvB_active
bP = bP /dvB_active
bPb = bPb/dvB_active
!bBeta = bP/bPb
bBeta = bBeta/dvB
bBeta = bBeta/dvB_active
if (present(bBetaO)) then
@@ -331,7 +342,73 @@ module streamline
end associate
end subroutine FLThermo
subroutine FLStdev(Model, ebGr, bTrc, bD, bP, stdD, stdP, sOpt)
!subroutine FLStdev(Model, ebGr, bTrc, stdD, stdP, sOpt)
!! Computes volume-weighted standard deviation of thermo quantities
type(chmpModel_T), intent(in) :: Model
type(ebGrid_T) , intent(in) :: ebGr
type(magLine_T ), intent(in) :: bTrc
real(rp), intent(in) :: bD, bP
! !! Field line averages calculated from FLThermo
real(rp), intent(out) :: stdD, stdP
!! Standard deviation of pressure and density we calculate and return
integer , intent(in ), optional :: sOpt
!! Which fluid to use
integer :: k,s0
real(rp) :: bMag,dl,eP,eD !Edge-centered values
real(rp) :: dvB
real(rp), dimension(NDIM) :: xyzTmp
!Zero out accumulators
stdP = 0.0
stdD = 0.0
dvB = 0.0
xyzTmp = 0.0
if (.not. bTrc%isGood) return
if (present(sOpt)) then
s0 = min(sOpt,Model%nSpc)
else
s0 = BLK
endif
associate(Np=>bTrc%Np,Nm=>bTrc%Nm)
! Now to stdev
!Loop over edges
do k=-Nm,Np-1
if ( (bTrc%ijk(k,XDIR) .le. ebGr%is+4) .or. (bTrc%ijk(k+1,XDIR) .le. ebGr%is+4) ) then
cycle
endif
!Get edge-centered quantities
dl = norm2(bTrc%xyz(k+1,:) - bTrc%xyz(k,:)) !Edge length
eD = 0.5*( bTrc%Gas(k+1,DEN ,s0) + bTrc%Gas(k,DEN ,s0) )
eP = 0.5*( bTrc%Gas(k+1,PRESSURE,s0) + bTrc%Gas(k,PRESSURE,s0) )
bMag = 0.5* (bTrc%magB(k+1) + bTrc%magB(k))
dvB = dvB + dl/bMag
! Perform volume-weighted summations
stdP = stdP + (eP - bP)**2 * dl/bMag
stdD = stdD + (eD - bD)**2 * dl/bMag
! Keep track of flux tube volume to normalize by later
enddo
dvB = max(dvB, TINY)
! Turn summations into standard deviations
stdP = sqrt(stdP / dvB)
stdD = sqrt(stdD / dvB)
end associate
end subroutine FLStdev
!Flux tube entropy (from BLK)
function FLEntropy(Model,ebGr,bTrc,GamO) result(S)
type(chmpModel_T), intent(in) :: Model

24
src/drivers/raijuOWDx.F90
View File

@@ -107,28 +107,38 @@ program raijuOWDx
call Tic("Output")
! Output if ready
if (raiApp%State%IO%doRestart(raiApp%State%t)) then
call Tic("RAIJU Restart")
call raijuResOutput(raiApp%Model,raiApp%Grid,raiApp%State)
call Toc("RAIJU Restart")
endif
if (raiApp%State%IO%doOutput(raiApp%State%t)) then
call Tic("RAIJU Output")
call raijuOutput(raiApp%Model,raiApp%Grid,raiApp%State)
call Toc("RAIJU Output")
if (ebModel%doEBOut) then
! Write eb at the same output cadence as imag
write(gStr,'(A,I0)') "Step#", ebModel%nOut
call Tic("CHIMP EB3D")
call writeEB3D(ebModel,ebState,gStr)
ebModel%nOut = ebModel%nOut + 1
ebModel%tOut = inTscl*raiApp%State%IO%tOut ! Idk if we need to set this since chimp isn't controlling its own output
call Toc("CHIMP EB3D")
endif
if (doFLOut) then
call Tic("RCM FLs")
call WriteRCMFLs(raijuCplBase%fromV%magLines,raiApp%State%IO%nOut, &
raiApp%State%mjd,raiApp%State%t, &
raiApp%Grid%shGrid%Nt,raiApp%Grid%shGrid%Np)
call Toc("RCM FLs")
endif
write(gStr,'(A,I0)') "Step#", raiApp%State%IO%nOut-1 ! nOut got advanced by raijuOutput above
call Tic("Remix rmReader")
call outputRMSG(rmReader,"rmReader.h5",.false., gStr)
call Toc("Remix rmReader")
endif
call Toc("Output")
@@ -158,13 +168,13 @@ program raijuOWDx
call Toc("RAIJU Advance")
isFirstCpl = .false.
write(*,*)raiApp%State%t
write(*,*)ebModel%t/inTscl
write(*,*)rmReader%time
write(*,*)"MJDs: "
write(*,*)raiApp%State%mjd
write(*,*)T2MJD((ebModel%t - ebState%ebTab%times(1))/inTscl, ebState%ebTab%MJDs(1))
write(*,*)T2MJD((rmReader%time - rmReader%rmTab%times(1)), rmReader%rmTab%MJDs(1))
!write(*,*)raiApp%State%t
!write(*,*)ebModel%t/inTscl
!write(*,*)rmReader%time
!write(*,*)"MJDs: "
!write(*,*)raiApp%State%mjd
!write(*,*)T2MJD((ebModel%t - ebState%ebTab%times(1))/inTscl, ebState%ebTab%MJDs(1))
!write(*,*)T2MJD((rmReader%time - rmReader%rmTab%times(1)), rmReader%rmTab%MJDs(1))
! Advance model times
raiApp%State%t = raiApp%State%t + raiApp%Model%dt

3
src/raiju/raijuAdvancer.F90
View File

@@ -116,7 +116,8 @@ module raijuAdvancer
! Calc next time step
dt = activeDt(Model, Grid, State, k)
!dt = activeDt(Model, Grid, State, k)
dt = activeDt_LR(Model, Grid, State, k)
! If needed, reduce dt to fit within remaining time
! This way, all channels will end exactly at tEnd

226
src/raiju/raijuDomain.F90
View File

@@ -2,6 +2,8 @@ module raijuDomain
use raijudefs
use raijuTypes
use raijugrids
use math
implicit none
@@ -34,51 +36,160 @@ module raijuDomain
end subroutine setActiveDomain
subroutine getInactiveCells(Model, Grid, State, isInactive)
subroutine getInactiveCells(Model, Grid, State, isInactive, isCoreInactiveO)
!! Applies series of criteria to determine which cells should be marked as inactive
type(raijuModel_T), intent(in) :: Model
type(raijuGrid_T ), intent(in) :: Grid
type(raijuState_T), intent(in) :: State
logical, dimension(Grid%shGrid%isg:Grid%shGrid%ieg,Grid%shGrid%jsg:Grid%shGrid%jeg), intent(inout) :: isInactive
!! The ultimate domain we want to call inactive, don't trust at all
logical, dimension(Grid%shGrid%isg:Grid%shGrid%ieg,Grid%shGrid%jsg:Grid%shGrid%jeg), &
optional, intent(inout) :: isCoreInactiveO
!! Debug option: Contians the inactive region according to only core physical constraints, no extra trimming
integer :: i,j
real(rp) :: xyMin
real(rp), dimension(2,2) :: bminSquare
logical :: iShellHasCheck
real(rp) :: bndRateLim
integer :: n_bndLim, n_contig
integer, dimension(Grid%shGrid%jsg:Grid%shGrid%jeg) :: bndLoc, bndR, bndL
real(rp), dimension(Grid%shGrid%isg:Grid%shGrid%ieg+1,Grid%shGrid%jsg:Grid%shGrid%jeg+1) :: cornerNormAngle
integer, dimension(Grid%shGrid%isg:Grid%shGrid%ieg,Grid%shGrid%jsg:Grid%shGrid%jeg) :: ocbDist
logical, dimension(Grid%shGrid%isg:Grid%shGrid%ieg,Grid%shGrid%jsg:Grid%shGrid%jeg) :: isCoreInactive, checkMask, isInactive_tmp
isInactive = .false.
do j=Grid%shGrid%jsg, Grid%shGrid%jeg
!do i=Grid%shGrid%ieg, Grid%shGrid%isg,-1
!! NOTE: idir in descending order
!! This means we are going from inner R boundary to outer R
do i=Grid%shGrid%isg, Grid%shGrid%ieg
checkMask = .false.
bndLoc = Grid%shGrid%isg
call calcCoreConstraints(Model, Grid, State, isCoreInactive)
isInactive = isCoreInactive
! Any cell with an open corner is bad
if (any(State%topo(i:i+1,j:j+1) .eq. RAIJUOPEN)) then
isInactive(i,j) = .true.
! Prep trim criteria
call calcCornerNormAngles(Model, Grid, State, cornerNormAngle)
! Now do trimming
associate(sh=>Grid%shGrid)
! Get dist from OCb
call CalcOCBDist(sh, isCoreInactive, 4, ocbDist)
! Highlight areas to start checking for extra constraints
where(ocbDist > 0 .and. ocbDist < 5)
checkMask = .true.
end where
do i=sh%isg, sh%ie ! NOTE: Not touching upper theta ghosts, we shouldn't reach there anyways. If we do then we have bigger problems
iShellHasCheck = .false.
if ( all(State%active(i,:) == RAIJUINACTIVE) ) then
cycle
endif
do j=sh%js, sh%je ! NOTE: Not touching j ghosts
if (checkMask(i,j)) then
iShellHasCheck = .true.
! Criteria check
if ( any(State%Pstd(i:i+1,j:j+1,0) > Model%PstdThresh) .or. any(cornerNormAngle(i:i+1,j:j+1) < Model%normAngThresh)) then
isInactive(i,j) = .true.
if (i > bndLoc(j)) then
bndLoc(j) = i
endif
! Flag points a little depper into the domain as places to check for criteria pass
checkMask(i+1,j-2:j+2) = .true.
checkMask(i+2,j-1:j+1) = .true.
checkMask(i+3,j ) = .true.
endif
endif
enddo
if (.not. iShellHasCheck) then
exit
endif
enddo
! Wrap in J
!call wrapJcc(sh, isInactive)
isInactive(:, sh%jsg:sh%js-1) = isInactive(:, sh%je-sh%Ngw+1:sh%je)
isInactive(:, sh%je+1:sh%jeg) = isInactive(:, sh%js:sh%js+sh%Nge-1)
bndLoc(sh%jsg:sh%js-1) = bndLoc(sh%je-sh%Ngw+1:sh%je)
bndLoc(sh%je+1:sh%jeg) = bndLoc(sh%js:sh%js+sh%Nge-1)
if (any(State%vaFrac(i:i+1,j:j+1) .le. Model%vaFracThresh)) then
isInactive(i,j) = .true.
bndL = bndLoc
bndR = bndLoc
bndRateLim = 0.45 ! del(theta) / del(phi)
n_bndLim = 3
!!NOTE: FIXME!!
!! Just hard-setting number of cells for testing purposed
write(*,*)"Bad hard-coded cell num for bndLoc smooth"
! Right sweep
do j=sh%js,sh%je
bndR(j) = max(bndR(j), bndR(j-1)-n_bndLim)
enddo
bndR(sh%je+1:sh%jeg) = bndR(sh%je) ! Extend final right sweep result into j ghosts
! Left sweep
do j=sh%je,sh%js, -1
bndL(j) = max(bndL(j), bndL(j+1)-n_bndLim)
enddo
bndL(sh%isg:sh%is-1) = bndL(sh%is) ! Extend final left sweep into ghosts
! Now combine back into bndLoc
do j=sh%jsg,sh%jeg
bndLoc(j) = max(bndL(j), bndR(j))
isInactive(sh%isg:bndLoc(j), j) = .true.
enddo
! Finally, kick out any stumpy regions
n_contig=4
isInactive_tmp = isInactive
do i=sh%isg, sh%ieg
do j=sh%js, sh%je
if (isInactive(i,j) .and. isInactive(i+n_contig,j)) then
isInactive_tmp(i:i+n_contig, j) = .true.
endif
bminSquare(1,1) = norm2(State%Bmin(i ,j ,:))
bminSquare(2,1) = norm2(State%Bmin(i+1,j ,:))
bminSquare(1,2) = norm2(State%Bmin(i ,j+1,:))
bminSquare(2,2) = norm2(State%Bmin(i+1,j+1,:))
if (any( bminSquare .le. Model%bminThresh) ) then
!if (any( norm2(State%Bmin(i:i+1,j:j+1,:),dim=3) .le. Model%bminThresh) ) then
isInactive(i,j) = .true.
if (isInactive(i,j) .and. isInactive(i,j+n_contig)) then
isInactive_tmp(i,j:j+n_contig) = .true.
endif
xyMin = sqrt(State%xyzMin(i,j,XDIR)**2 + State%xyzMin(i,j,YDIR)**2)
! Simple circle limit
if ( (xyMin > Model%maxTail_buffer) .or. xyMin > Model%maxSun_buffer) then
isInactive(i,j) = .true.
endif
enddo
enddo
isInactive = isInactive_tmp
end associate
contains
subroutine calcCoreConstraints(Model, Grid, State, isCoreInactive)
type(raijuModel_T), intent(in) :: Model
type(raijuGrid_T ), intent(in) :: Grid
type(raijuState_T), intent(in) :: State
logical, dimension(Grid%shGrid%isg:Grid%shGrid%ieg,Grid%shGrid%jsg:Grid%shGrid%jeg), intent(inout) :: isCoreInactive
integer :: i,j
real(rp) :: xyMin
real(rp), dimension(2,2) :: bminSquare
isCoreInactive = .false.
do j=Grid%shGrid%jsg, Grid%shGrid%jeg
do i=Grid%shGrid%isg, Grid%shGrid%ieg
! Any cell with an open corner is bad
if (any(State%topo(i:i+1,j:j+1) .eq. RAIJUOPEN)) then
isCoreInactive(i,j) = .true.
endif
if (any(State%vaFrac(i:i+1,j:j+1) .le. Model%vaFracThresh)) then
isCoreInactive(i,j) = .true.
endif
bminSquare(1,1) = norm2(State%Bmin(i ,j ,:))
bminSquare(2,1) = norm2(State%Bmin(i+1,j ,:))
bminSquare(1,2) = norm2(State%Bmin(i ,j+1,:))
bminSquare(2,2) = norm2(State%Bmin(i+1,j+1,:))
if (any( bminSquare .le. Model%bminThresh) ) then
!if (any( norm2(State%Bmin(i:i+1,j:j+1,:),dim=3) .le. Model%bminThresh) ) then
isCoreInactive(i,j) = .true.
endif
xyMin = sqrt(State%xyzMin(i,j,XDIR)**2 + State%xyzMin(i,j,YDIR)**2)
! Simple circle limit
if ( (xyMin > Model%maxTail_buffer) .or. xyMin > Model%maxSun_buffer) then
isCoreInactive(i,j) = .true.
endif
enddo
enddo
end subroutine calcCoreConstraints
end subroutine getInactiveCells
@@ -220,4 +331,61 @@ module raijuDomain
enddo
end subroutine calcMapJacNorm
subroutine calcCornerNormAngles(Model, Grid, State, normAngle)
!! For each cell corner, calculate the maximum angle between the normals of the 4 triangles its a part of
type(raijuModel_T), intent(in) :: Model
type(raijuGrid_T ), intent(in) :: Grid
type(raijuState_T), intent(in) :: State
real(rp), dimension(Grid%shGrid%isg:Grid%shGrid%ieg+1,Grid%shGrid%jsg:Grid%shGrid%jeg+1), intent(inout) :: normAngle
integer :: i,j,d,u,v
real(rp), dimension(3) :: v1, v2, v3, v4
real(rp), dimension(4,3) :: crosses
real(rp), dimension(Grid%shGrid%isg:Grid%shGrid%ieg+1,Grid%shGrid%jsg:Grid%shGrid%jeg+1,3) :: xyz_sm
normAngle = 1.0
! First, smooth xyzmins
xyz_sm = State%xyzMin(:,:,:)
do i=Grid%shGrid%isg+1,Grid%shGrid%ieg ! Ignore last row and column, no smoothing for you
do j=Grid%shGrid%jsg+1,Grid%shGrid%jeg
do d=XDIR,ZDIR
xyz_sm(i,j,d) = SmoothOperator33(State%xyzMin(i-1:i+1,j-1:j+1,d))
enddo
enddo
enddo
do i=Grid%shGrid%isg+1,Grid%shGrid%ieg
do j=Grid%shGrid%jsg+1,Grid%shGrid%jeg
! Build vectors from corner i,j
do d=XDIR,ZDIR
v1(d) = xyz_sm(i,j,d) - xyz_sm(i ,j-1,d)
v2(d) = xyz_sm(i,j,d) - xyz_sm(i-1,j ,d)
v3(d) = xyz_sm(i,j,d) - xyz_sm(i ,j+1,d)
v4(d) = xyz_sm(i,j,d) - xyz_sm(i+1,j ,d)
enddo
! Do cross products
crosses(1,:) = cross(v2, v1)
crosses(2,:) = cross(v3, v2)
crosses(3,:) = cross(v4, v3)
crosses(4,:) = cross(v1, v4)
crosses(1,:) = crosses(1,:)/max(norm2(crosses(1,:)), TINY)
crosses(2,:) = crosses(2,:)/max(norm2(crosses(2,:)), TINY)
crosses(3,:) = crosses(3,:)/max(norm2(crosses(3,:)), TINY)
crosses(4,:) = crosses(4,:)/max(norm2(crosses(4,:)), TINY)
! Calculate min
do u=1,3
do v=u+1,4
! We are doing dot product, so 1 = zero angle and -1 means 180 deg angle
normAngle(i,j) = min(normAngle(i,j), dot_product(crosses(u,:), crosses(v,:)))
enddo
enddo
enddo
enddo
end subroutine calcCornerNormAngles
end module raijuDomain

3
src/raiju/raijuEtautils.F90
View File

@@ -35,6 +35,9 @@ module raijuetautils
State%vAvg = 0.0
associate (shG => Grid%shGrid, spc => Grid%spc)
!$OMP PARALLEL DO default(shared) &
!$OMP schedule(dynamic) &
!$OMP private(i,j,s)
do j=shG%jsg,shG%jeg
do i=shG%isg,shG%ieg
do s=1,Grid%nSpc

19
src/raiju/raijuIO.F90
View File

@@ -99,6 +99,8 @@ module raijuIO
call AddOutVar(IOVars,"doExcesstoPsph",Model%doExcesstoPsph) ! Attr
call AddOutVar(IOVars,"doPlasmasphere",Model%doPlasmasphere) ! Attr
call AddOutVar(IOVars,"doActiveShell",Model%doActiveShell) ! Attr
call AddOutVar(IOVars,"nSpcIn",Model%nSpcMHD) ! Attr
call AddOutVar(IOVars,"nSpcRAIJU",Model%nSpc) ! Attr
call WriteVars(IOVars,.true.,Model%raijuH5)
! Root data
@@ -259,6 +261,8 @@ module raijuIO
call AddOutVar(IOVars,"vaFrac" ,State%vaFrac (is:ie+1,js:je+1) ,uStr="fraction",dStr="Fraction of Alfven speed over magnetofast + ExB speed")
call AddOutVar(IOVars,"Pavg_in",State%Pavg (is:ie,js:je, :) ,uStr="nPa" ,dStr="Pressures from imagtubes")
call AddOutVar(IOVars,"Davg_in",State%Davg (is:ie,js:je, :) ,uStr="#/cc",dStr="Densities from imagtubes")
call AddOutVar(IOVars,"Pstd_in",State%Pstd (is:ie,js:je, :) ,uStr="normalized" ,dStr="Std. dev. of species pressure from imagtubes")
call AddOutVar(IOVars,"Dstd_in",State%Dstd (is:ie,js:je, :) ,uStr="normalized" ,dStr="Std. dev. of species density from imagtubes")
! Idk about you but I did not expect true to equal -1
!allocate(outActiveShell(is:ie, Grid%Nk))
@@ -297,12 +301,15 @@ module raijuIO
!allocate(outEnt(is:ie,js:je))
allocate(outTmp2D(is:ie,js:je))
outTmp2D = -1.0
do i=is,ie
do j=js,je
if (all(State%bVol(i:i+1,j:j+1) > 0)) then
bvol_cc = 0.25*(State%bVol(i,j) + State%bVol(i+1,j) + State%bVol(i,j+1) + State%bVol(i+1,j+1))
outTmp2D(i,j) = State%Press(i,j,1)*bVol_cc**(5./3.)
endif
!$OMP PARALLEL DO default(shared) &
!$OMP schedule(dynamic) &
!$OMP private(i,j)
do j=js,je
do i=is,ie
!if (all(State%bVol(i:i+1,j:j+1) > 0)) then
!bvol_cc = 0.25*(State%bVol(i,j) + State%bVol(i+1,j) + State%bVol(i,j+1) + State%bVol(i+1,j+1))
outTmp2D(i,j) = State%Press(i,j,1)*State%bvol_cc(i,j)**(5./3.)
!endif
enddo
enddo
call AddOutVar(IOVars,"FTEntropy",outTmp2D,uStr="nPa*(Rp/nT)^(5/3)")

96
src/raiju/raijuPreAdvancer.F90
View File

@@ -75,6 +75,7 @@ module raijuPreAdvancer
! Calc sub-time step. Each channel will do this on its own, but this way we can output the step sizes everyone is using
!State%dtk(k) = activeDt(Model, Grid, State, k)
call reconVelocityLRs(Model, Grid, State, k, State%iVelL(:,:,k,:), State%iVelR(:,:,k,:))
State%dtk(k) = activeDt_LR(Model, Grid, State, k)
enddo
call Toc("Calc cell-centered velocities")
@@ -799,6 +800,7 @@ module raijuPreAdvancer
Grid%shGrid%jsg:Grid%shGrid%jeg) :: isGCC
real(rp), dimension(Grid%shGrid%isg:Grid%shGrid%ieg+1,&
Grid%shGrid%jsg:Grid%shGrid%jeg+1,2) :: tmpVelL, tmpVelR
integer :: i,j
iVelL = 0.0
iVelR = 0.0
@@ -819,6 +821,25 @@ module raijuPreAdvancer
iVelR(:,:,RAI_PH) = tmpVelR(:,:,RAI_PH)
! Hax
!!$OMP PARALLEL DO default(shared) &
!!$OMP schedule(dynamic) &
!!$OMP private(i,j)
!do j = Grid%shGrid%js,Grid%shGrid%je+1
! do i = Grid%shGrid%is,Grid%shGrid%ie+1
! if (State%active(i,j) .ne. RAIJUACTIVE) then
! cycle
! endif
! if (State%active(i-1,j) .eq. RAIJUBUFFER) then
! iVelL(i,j,RAI_TH) = 0.0
! iVelR(i,j,RAI_TH) = 0.0
! endif
! if (State%active(i,j-1) .eq. RAIJUBUFFER) then
! iVelL(i,j,RAI_PH) = 0.0
! iVelR(i,j,RAI_PH) = 0.0
! endif
! enddo
!enddo
end subroutine reconVelocityLRs
!------
@@ -842,9 +863,9 @@ module raijuPreAdvancer
dtArr = HUGE
!$OMP PARALLEL DO default(shared) &
!$OMP schedule(dynamic) &
!$OMP private(i,j, velij_th, velij_ph, dl)
!!$OMP PARALLEL DO default(shared) &
!!$OMP schedule(dynamic) &
!!$OMP private(i,j, velij_th, velij_ph, dl)
do j=sh%js,sh%je+1
do i=sh%is,sh%ie+1
velij_th = TINY
@@ -895,6 +916,75 @@ module raijuPreAdvancer
end function activeDt
function activeDt_LR(Model, Grid, State, k) result(dt)
!! Calculates min dt needed to safele evolve active domain for given energy invariant channel k
!! TODO: Consider dynamic CFL factor
type(raijuModel_T), intent(in) :: Model
type(raijuGrid_T ), intent(in) :: Grid
type(raijuState_T), intent(in) :: State
integer, intent(in) :: k
integer :: i,j
real(rp) :: dl, velij_th, velij_ph
real(rp), dimension(Grid%shGrid%is:Grid%shGrid%ie+1,Grid%shGrid%js:Grid%shGrid%je+1) :: dtArr
real(rp) :: dt
associate (sh => Grid%shGrid)
dtArr = HUGE
!!$OMP PARALLEL DO default(shared) &
!!$OMP schedule(dynamic) &
!!$OMP private(i,j, velij_th, velij_ph, dl)
do j=sh%js,sh%je+1
do i=sh%is,sh%ie+1
velij_th = TINY
velij_ph = TINY
! NOTE: We are only checking faces bordering non-ghost cells because those are the only ones we use for evolution
! TODO: Strictly speaking, there are some faces that are included here that shouldn't be because they're never used to evolve anything
! so we should make the loop js:je, is:ie, and handle the last row and column afterwards
! We are responsible for face (i,j,TH) and (i,j,PH)
! Theta faces first
! This is a weird pattern. Basically, we can't cycle if the overall desired condition isn't met, because we are doing theta and phi directions in the same loop
! At the same time, we don't want to write one massive if condition with all options covered, or a bunch of nested if statements
! So instead, we break them up, such that if a single if condition is true it means we have failed the physical condition for us to calculate a valid timestep
! If all if conditions fail then the physical conditions are met and we can do our calculations in the else block
if (Model%doActiveShell .and. ( .not. State%activeShells(i-1,k) .or. .not. State%activeShells(i,k)) ) then
! In order for a theta-dir face to be usable, we need both sides to be active
continue
else if ( State%active(i-1,j) .ne. RAIJUACTIVE .or. State%active(i,j) .ne. RAIJUACTIVE ) then
continue
else
! We are good lets calculate a dt for this face
velij_th = max(abs(State%iVelL(i,j,k,RAI_TH)), abs(State%iVelR(i,j,k,RAI_TH)), TINY) ! [m/s]
!dtArr(i,j,RAI_TH) = ( Grid%delTh(i) * Model%planet%ri_m ) / velij ! [s]
endif
! Phi faces
if (Model%doActiveShell .and. .not. State%activeShells(i,k)) then
! In order for a phi-dir face to be usable, we just need this i shell to be active
continue
else if (State%active(i,j-1) .ne. RAIJUACTIVE .or. State%active(i,j) .ne. RAIJUACTIVE ) then
continue
else
velij_ph = max(abs(State%iVelL(i,j,k,RAI_PH)), abs(State%iVelR(i,j,k,RAI_PH)), TINY)
!dtArr(i,j,RAI_PH) = ( Grid%delPh(j) * Model%planet%ri_m ) / velij ! [s]
endif
dl = abs(min(Grid%delTh(i), Grid%delPh(j)*sin(sh%thc(i))))
dtArr(i,j) = (dl*Model%planet%ri_m) / sqrt(velij_th**2 + velij_ph**2)
enddo
enddo
dt = Model%CFL*minval(dtArr)
end associate
end function activeDt_LR
!------
! Extras
!------

14
src/raiju/raijuRecon.F90
View File

@@ -319,10 +319,10 @@ module raijuRecon
if (present(QreconLO) .and. present(QreconRO)) then
QreconLO = 0.0
QreconRO = 0.0
!$OMP PARALLEL DO default(shared) &
!$OMP schedule(dynamic) &
!$OMP private(i,j) &
!$OMP IF(doOMP)
!!$OMP PARALLEL DO default(shared) &
!!$OMP schedule(dynamic) &
!!$OMP private(i,j) &
!!$OMP IF(doOMP)
do j=Grid%shGrid%js,Grid%shGrid%je+1
do i=Grid%shGrid%is,Grid%shGrid%ie+1
! Theta dir
@@ -344,6 +344,10 @@ module raijuRecon
enddo
else
!!$OMP PARALLEL DO default(shared) &
!!$OMP schedule(dynamic) &
!!$OMP private(i,j) &
!!$OMP IF(doOMP)
do j=Grid%shGrid%js,Grid%shGrid%je+1
do i=Grid%shGrid%is,Grid%shGrid%ie+1
! Theta dir
@@ -460,7 +464,7 @@ module raijuRecon
if (Model%doUseVelLRs) then
do j=Grid%shGrid%js,Grid%shGrid%je+1
do i=Grid%shGrid%is,Grid%shGrid%ie+1
! Only allow elocities passing through the interface
! Only allow velocities passing through the interface
QfluxL = QfaceL(i,j,RAI_TH)*max(0.0, State%iVelL(i,j,k,RAI_TH))
QfluxR = QfaceR(i,j,RAI_TH)*min(0.0, State%iVelR(i,j,k,RAI_TH))
Qflux(i,j,RAI_TH) = QfluxL + QfluxR

12
src/raiju/raijuStarter.F90
View File

@@ -238,8 +238,12 @@ module raijustarter
! Domain determination
call iXML%Set_Val(Model%vaFracThresh, "cpl/vaFracThresh", def_vaFracThresh)
call iXML%Set_Val(Model%bminThresh, "cpl/bminThresh", def_bminThresh)
call iXML%Set_Val(Model%vaFracThresh , "cpl/vaFracThresh" , def_vaFracThresh)
call iXML%Set_Val(Model%bminThresh , "cpl/bminThresh" , def_bminThresh)
call iXML%Set_Val(Model%PstdThresh , "cpl/Pstd" , def_PstdThresh)
call iXML%Set_Val(Model%normAngThresh, "cpl/normAngThresh", def_normAngle)
! Convert degrees to dot product value
Model%normAngThresh = cos(Model%normAngThresh*PI/180.0)
! Fluid mapping
call iXML%Set_Val(nFluids, "cpl/nFluidsIn",0)
@@ -366,6 +370,8 @@ module raijustarter
! Coupling input moments
allocate( State%Pavg(sh%isg:sh%ieg , sh%jsg:sh%jeg, 0:Grid%nFluidIn) )
allocate( State%Davg(sh%isg:sh%ieg , sh%jsg:sh%jeg, 0:Grid%nFluidIn) )
allocate( State%Pstd(sh%isg:sh%ieg , sh%jsg:sh%jeg, 0:Grid%nFluidIn) )
allocate( State%Dstd(sh%isg:sh%ieg , sh%jsg:sh%jeg, 0:Grid%nFluidIn) )
! Bmin surface
allocate( State%Bmin (sh%isg:sh%ieg+1, sh%jsg:sh%jeg+1, 3 ) )
allocate( State%xyzMin (sh%isg:sh%ieg+1, sh%jsg:sh%jeg+1, 3 ) )
@@ -378,6 +384,8 @@ module raijustarter
allocate( State%bvol (sh%isg:sh%ieg+1, sh%jsg:sh%jeg+1) )
allocate( State%bvol_cc(sh%isg:sh%ieg , sh%jsg:sh%jeg ) )
allocate( State%vaFrac(sh%isg:sh%ieg+1, sh%jsg:sh%jeg+1) )
! 1D cell-centered quantities
allocate( State%bndLoc(sh%jsg:sh%jeg) )
! 2D cell-centered quantities
allocate( State%active (sh%isg:sh%ieg, sh%jsg:sh%jeg) )
allocate( State%active_last (sh%isg:sh%ieg, sh%jsg:sh%jeg) )

9
src/voltron/imagtubes.F90
View File

@@ -136,7 +136,7 @@ module imagtubes
real(rp), dimension(NDIM) :: x0, bEq, xyzIon
real(rp), dimension(NDIM) :: xyzC,xyzIonC
integer :: OCb
real(rp) :: bD,bP,dvB,bBeta,rCurv, bDRC, bPRC, bBetaRC
real(rp) :: bD,bP,dvB,bBeta,rCurv, bDRC, bPRC, bBetaRC, stdP, stdD
real(rp) :: VaMKS,CsMKS,VebMKS !Speeds in km/s
real(rp) :: TiEV !Temperature in ev
@@ -187,6 +187,8 @@ module imagtubes
ijTube%Vol = -1.0
ijTube%Pave = 0.0
ijTube%Nave = 0.0
ijTube%Pstd = 0.0
ijTube%Nstd = 0.0
ijTube%beta_average = 0.0
ijTube%latc = 0.0
ijTube%lonc = 0.0
@@ -210,11 +212,14 @@ module imagtubes
!dvB = Flux-tube volume (Re/EB)
!write(*,*)"FLThermo, s=",s
call FLThermo(ebModel,ebGr,bTrc,bD,bP,dvB,bBeta,s)
call FLStdev (ebModel,ebGr,bTrc,bD,bP,stdD,stdP,s)
!call FLStdev (ebModel,ebGr,bTrc,stdD,stdP,s)
bP = bP*1.0e-9 !nPa=>Pa
bD = bD*1.0e+6 !#/cc => #/m3
ijTube%Pave(s) = bP
ijTube%Nave(s) = bD
ijTube%Pstd(s) = stdP * 1.0e-9 ! nPa=>Pa
ijTube%Nstd(s) = stdD * 1.0e+6 ! #/cc=>#/m3
if (s .eq. RCFLUID) then
bPRC = bP
bDRC = bD

23
src/voltron/raijuInterfaces/raijuCplHelpers.F90
View File

@@ -24,7 +24,7 @@ module raijuCplHelpers
integer :: i,j,s
real(rp) :: P, D
real(rp) :: P, D, Pstd, Dstd
real(rp), dimension(Grid%shGrid%isg:Grid%shGrid%ieg+1) :: bVol_dip_corner
real(rp), dimension(Grid%shGrid%isg:Grid%shGrid%ieg) :: bVol_dip_cc
real(rp), dimension(2,2) :: dBVol
@@ -88,19 +88,28 @@ module raijuCplHelpers
+ ijTubes(i+1,j)%Pave(s) + ijTubes(i+1,j+1)%Pave(s)) * 1.0D+9 ! [Pa -> nPa]
D = 0.25*(ijTubes(i ,j)%Nave(s) + ijTubes(i ,j+1)%Nave(s) &
+ ijTubes(i+1,j)%Nave(s) + ijTubes(i+1,j+1)%Nave(s)) * 1.0D-6 ! [#/m^3 --> #/cc]
Pstd = 0.25*(ijTubes(i ,j)%Pstd(s) + ijTubes(i ,j+1)%Pstd(s) &
+ ijTubes(i+1,j)%Pstd(s) + ijTubes(i+1,j+1)%Pstd(s)) * 1.0D+9 ! [Pa -> nPa]
Dstd = 0.25*(ijTubes(i ,j)%Nstd(s) + ijTubes(i ,j+1)%Nstd(s) &
+ ijTubes(i+1,j)%Nstd(s) + ijTubes(i+1,j+1)%Nstd(s)) * 1.0D+9 ! [Pa -> nPa]
State%Pavg(i,j,s) = P
State%Davg(i,j,s) = D
State%Pstd(i,j,s) = Pstd / max(P, TINY)
State%Dstd(i,j,s) = Dstd / max(D, TINY)
!State%bvol_cc(i,j) = toCenter2D(State%bvol(i:i+1,j:j+1))
! Before we average, take dipole out so we don't get poor averaging of background
dBVol = State%bvol(i:i+1,j:j+1)
dBVol(:,1) = dBVol(:,1) - bVol_dip_corner(i:i+1)
dBVol(:,2) = dBVol(:,2) - bVol_dip_corner(i:i+1)
State%bvol_cc(i,j) = toCenter2D(dBVol) + bVol_dip_cc(i)
enddo
dBVol = State%bvol(i:i+1,j:j+1)
dBVol(:,1) = dBVol(:,1) - bVol_dip_corner(i:i+1)
dBVol(:,2) = dBVol(:,2) - bVol_dip_corner(i:i+1)
State%bvol_cc(i,j) = toCenter2D(dBVol) + bVol_dip_cc(i)
! Do our own "wIMAG" calculation here so we ensure we use the fluid we want to (Bulk)
! Calculate the fraction of Alfven speed to total velocity
!VaMKS = flux tube arc length [km] / Alfven crossing time [s]