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

Merged development into dev_atiken_update

Browse Source
This commit is contained in:
Nikhil Rao
2025-11-17 14:38:27 +00:00
parent a386f4d9b9 76a31f7526
commit 646ed1513c
25 changed files with 1113 additions and 921 deletions
Download Patch File Download Diff File Expand all files Collapse all files

3
.gitmodules vendored
View File

@@ -1,3 +0,0 @@
[submodule "wiki"]
path = wiki
url = https://bitbucket.org/aplkaiju/kaiju/wiki

5
.mailmap
View File

@@ -1,5 +0,0 @@
Kareem <kareem.sorathia@gmail.com>
Kareem <kareem.sorathia@gmail.com> <soratka1@jhuapl.edu>
Jeffrey Garretson <Jeffrey.Garretson@jhuapl.edu> <jeffreyg@cheyenne1.ib0.cheyenne.ucar.edu>
Jeffrey Garretson <Jeffrey.Garretson@jhuapl.edu> <jeffreyg@cheyenne3.ib0.cheyenne.ucar.edu>
Jeffrey Garretson <Jeffrey.Garretson@jhuapl.edu> <jeffreyg@cheyenne4.ib0.cheyenne.ucar.edu>

3
CMakeLists.txt
View File

@@ -1,5 +1,6 @@
cmake_minimum_required(VERSION 3.20.2)
project(Kaiju Fortran)
project(Kaiju Fortran C)
#K: Adding C to project to deal w/ issues w/ most recent HDF (10/13/25)
# add and search for pfunit (fingers crossed)
list(APPEND CMAKE_PREFIX_PATH "./external")

24
cmake/compilers.cmake
View File

@@ -83,7 +83,12 @@ if(CMAKE_Fortran_COMPILER_ID MATCHES Intel)
#Base
string(APPEND CMAKE_Fortran_FLAGS " -fPIC -fpconstant")
#Production
set(PROD "-align array64byte -align rec32byte -no-prec-div -fast-transcendentals")
set(PROD "-align array64byte -align rec32byte -no-prec-div")
if (CMAKE_Fortran_COMPILER_VERSION VERSION_LESS 23.1)
#Fast transcendentals removed in ifx
string(APPEND PROD " -fast-transcendentals")
endif()
#Production with Debug Info
set(PRODWITHDEBUGINFO "-traceback -debug all -align array64byte -align rec32byte -no-prec-div -fast-transcendentals")
#Debug
@@ -123,11 +128,24 @@ if(CMAKE_Fortran_COMPILER_ID MATCHES Intel)
endif()
string(APPEND PROD " -march=core-avx2")
string(APPEND PRODWITHDEBUGINFO " -march=core-avx2")
elseif (HOST MATCHES stampede3)
message("You're on Stampede3!")
if (ENABLE_MKL)
string(APPEND CMAKE_Fortran_FLAGS " -qmkl")
endif()
string(APPEND PROD " -xCORE-AVX512")
string(APPEND PRODWITHDEBUGINFO " -xCORE-AVX512")
endif()
#Check Intel Fortran version
if(NOT ALLOW_INVALID_COMPILERS AND CMAKE_Fortran_COMPILER_VERSION VERSION_GREATER "2021.9")
message(FATAL_ERROR "Intel Fortran compilers newer than 2023 (version 2021.8) are not supported. Set the ALLOW_INVALID_COMPILERS variable to ON to force compilation at your own risk.")
if(NOT ALLOW_INVALID_COMPILERS AND CMAKE_Fortran_COMPILER_VERSION VERSION_GREATER "2021.9" AND CMAKE_Fortran_COMPILER_VERSION VERSION_LESS "2025.1")
message(FATAL_ERROR "Intel OneAPI compilers between 2022-2024 have compiler bugs which cause weird numerical errors in our code. You can set the ALLOW_INVALID_COMPILERS variable to ON to force compilation at your own risk. You'll probably get what you deserve.")
endif()
#Check for MKL + intel25
if (ENABLE_MKL AND (CMAKE_Fortran_COMPILER_VERSION VERSION_GREATER "2024"))
message(WARNING "Intel OneAPI MKL has been found to fail in weird ways and should probably be avoided. But hey, do what you want. I'm a warning message, not a cop.")
endif()
elseif(CMAKE_Fortran_COMPILER_ID MATCHES GNU)

24
src/base/earthhelper.F90
View File

@@ -61,13 +61,24 @@ module earthhelper
real(rp), intent(in) :: kp
!Can only be 1,5
if ( (kp>=1.0) .and. (kp<=5.0) ) then
kpDefault = kp
else if (kp>5.0) then
kpDefault = 5.0
endif
kpDefault = GallagherKp(kp)
end subroutine SetKp0
! Routine to enforce kp to Gallagher's valid range [1,5]
pure function GallagherKp(kp) result(validKp)
real(rp), intent(in) :: kp
real(rp) :: validKp
if (kp < 1.0_rp) then
validKp = 1.0_rp
else if (kp > 5.0_rp) then
validKp = 5.0_rp
else
validKp = kp
endif
end function GallagherKp
!Return 2D gallagher density afa r,phi (rad)
function GallagherRP(r,phi,kpO) result(D)
real(rp), intent(in) :: r,phi
@@ -97,7 +108,8 @@ module earthhelper
D = 0.0
if (present(kpO)) then
kp = nint(kpO) ! Gallagher model takes integer Kp
!Can only be 1,5
kp = nint(GallagherKp(kpO))
else
kp = nint(kpDefault)
endif

71
src/chimp/chopio.F90
View File

@@ -16,6 +16,8 @@ module chopio
integer, parameter :: MAXEBVS = 50
integer :: Nx1 = 64, Nx2 = 64, Nx3 = 64
real(rp) :: originOff = 0.0 !Default origin offset
real(rp) :: lookDir = 0.0 !Default CONE grid tilt
real(rp) :: xyzMax = 10.0 !Default bound
real(rp), dimension(:,:,:), allocatable :: xxi,yyi,zzi,xxc,yyc,zzc
@@ -37,6 +39,10 @@ module chopio
integer :: iS,iE
real(rp) :: xcc(NDIM)
logical :: doLogR
real(rp) :: xOff,yOff,zOff !Cartesian offsets for CONE grid origin
real(rp) :: theta_pitch,phi_yaw !Theta and phi for y and z rotations in 3D space
real(rp) :: xxi_temp,yyi_temp,zzi_temp !Temp variables for 3D rotations
real(rp) :: xxi_tilt,yyi_tilt,zzi_tilt,xxi_pan,yyi_pan,zzi_pan !More temp variables for 3D rotations
!Check for time parallelism
call InitParInTime(Model,inpXML,"eb3",eb3DOutF)
@@ -61,6 +67,15 @@ module chopio
call inpXML%Set_Val(Nx2,'chop/Nx2',Nx2)
call inpXML%Set_Val(Nx3,'chop/Nx3',Nx3)
!Get origin offset for CONE chopping. Default to originOff (zero)
call inpXML%Set_Val(xOff,'chop/xOff',originOff)
call inpXML%Set_Val(yOff,'chop/yOff',originOff)
call inpXML%Set_Val(zOff,'chop/zOff',originOff)
!Get pitch and yaw angles for CONE chopping. Default to lookDir (zero)
call inpXML%Set_Val(theta_pitch,'chop/lookPol',lookDir)
call inpXML%Set_Val(phi_yaw,'chop/lookAzi',lookDir)
!for LFM grid, take all i-shells within x1Max, in the Sun direction
select case(trim(toUpper(idStr)))
!---------
@@ -120,11 +135,12 @@ module chopio
endif
dx2 = (x2Max-x2Min)*deg2rad/Nx2
dx3 = (x3Max-x3Min)*deg2rad/Nx3
!$OMP PARALLEL DO
!$OMP PARALLEL DO default(shared) &
!$OMP private(i,j,k,x1,x2,x3)
do k=1, Nx3+1
do j=1,Nx2+1
do i=1,Nx1+1
!x1 = Rin + (i-1)*dx1
!x1 = Rin(i,j,k) + (i-1)*dx1
if (doLogR) then
x1 = 10**( log10(x1Min) + (i-1)*dx1 )
else
@@ -141,6 +157,57 @@ module chopio
enddo
enddo
!---------
! Conical/Solid angle (based on RTP)
case("CONE")
!Check for log spacing in r
call inpXML%Set_Val(doLogR,'chop/doLogR',.false.)
if (doLogR) then
dx1 = ( log10(x1Max)-log10(x1Min) )/Nx1
else
dx1 = (x1Max - x1Min)/Nx1
endif
dx2 = (x2Max-x2Min)*deg2rad/Nx2
dx3 = (x3Max-x3Min)*deg2rad/Nx3
!$OMP PARALLEL DO default(shared) &
!$OMP private(i,j,k,x1,x2,x3,xxi_temp,yyi_temp,zzi_temp,xxi_tilt,yyi_tilt,zzi_tilt,xxi_pan,yyi_pan,zzi_pan)
do k=1,Nx3+1
do j=1,Nx2+1
do i=1,Nx1+1
!x1 = Rin(i,j,k) + (i-1)*dx1
if (doLogR) then
x1 = 10**( log10(x1Min) + (i-1)*dx1 )
else
x1 = x1Min + (i-1)*dx1
endif
x2 = x2Min*deg2rad + (j-1)*dx2 !Theta
x3 = x3Min*deg2rad + (k-1)*dx3 !Phi
! pre-rotation variables
xxi_temp = x1*cos(x3)*sin(x2)
yyi_temp = x1*sin(x3)*sin(x2)
zzi_temp = x1*cos(x2)
! first tilt (polar) then pan (azimuthal) (pitch then yaw)
! perform rotation about y axis (theta, polar angle)
xxi_tilt = (cos(theta_pitch*deg2rad)*xxi_temp) + (0) + (sin(theta_pitch*deg2rad)*zzi_temp)
yyi_tilt = (0) + (yyi_temp) + (0)
zzi_tilt = (-sin(theta_pitch*deg2rad)*xxi_temp) + (0) + (cos(theta_pitch*deg2rad)*zzi_temp)
! perform rotation about z axis (phi, azimuthal angle)
! apply to the tilted xxi/yyi/zzi
xxi_pan = (cos(phi_yaw*deg2rad)*xxi_tilt) + (-sin(phi_yaw*deg2rad)*yyi_tilt) + (0)
yyi_pan = (sin(phi_yaw*deg2rad)*xxi_tilt) + (cos(phi_yaw*deg2rad)*yyi_tilt) + (0)
zzi_pan = (0) + (0) + (zzi_tilt)
! post-rotation, add cartesian offsets
xxi(i,j,k) = xxi_pan + xOff
yyi(i,j,k) = yyi_pan + yOff
zzi(i,j,k) = zzi_pan + zOff
enddo
enddo
enddo
!---------
! only x1Max is used to set how far downtail you chop
case("LFM")
!$OMP PARALLEL DO

12
src/gamera/ringutils.F90
View File

@@ -491,8 +491,9 @@ module ringutils
!Map i to itself
ip = i
!Next do k, map via periodicity
!NOTE: This is assuming you have all
!Next do k, map via periodicity. Have to do k first otherwise have to deal w/ right hand becoming left
!NOTE: This is assuming you have all k cells (ie, no mpi decomp in KDIR)
if (k < Grid%ks) then
kp = k + Np
elseif (k > Grid%ke) then
@@ -501,18 +502,19 @@ module ringutils
kp = k
endif
!Finally do j
!Now handle ip,j,kp => ip,jp,kp
jp = j ! default value
if ( Model%Ring%doS .and. (j<Grid%js) ) then
!js-1 => js
jp = Grid%js + (Grid%js-j) - 1
kp = WrapK(k,Np)
kp = WrapK(kp,Np)
endif
if ( Model%Ring%doE .and. (j>Grid%je) ) then
!je+1 => je
jp = Grid%je - (j-Grid%je) + 1
kp = WrapK(k,Np)
kp = WrapK(kp,Np)
endif
end subroutine lfmIJKcc

54
src/voltron/modelInterfaces/imag2mhd_interface.F90
View File

@@ -79,13 +79,11 @@ module imag2mhd_interface
!Below inner boundary, do dipole projection
isGoodCC(i,j,k) = .true.
xyz = Gr%xyzcc(i,j,k,:) !Gamera grid center
call Proj2Rad(xyz,Rion,x1,x2)
Gr%Gas0(i,j,k,PROJLAT) = x1
call NHProj(xyz,x1,x2)
Gr%Gas0(i,j,k,PROJLAT) = x1 !Must project to NH
Gr%Gas0(i,j,k,PROJLON) = x2
else
!Get value from xyzsquish
if ( all(vApp%chmp2mhd%isGood(i:i+1,j:j+1,k:k+1)) ) then
!All values are good, so just do this thing
call SquishCorners(vApp%chmp2mhd%xyzSquish(i:i+1,j:j+1,k:k+1,1),Qs)
@@ -123,7 +121,7 @@ module imag2mhd_interface
if (i < Gr%is) then
!Use dipole projection
xyz = Gr%xyz(i,j,k,:) !Gamera grid corner
call Proj2Rad(xyz,Rion,x1,x2)
call NHProj(xyz,x1,x2)
isG = .true.
else
x1 = vApp%chmp2mhd%xyzSquish(i,j,k,1)
@@ -284,41 +282,43 @@ module imag2mhd_interface
real(rp), intent(inout) :: Q(Gr%isg:Gr%ieg,Gr%jsg:Gr%jeg,Gr%ksg:Gr%keg)
integer :: i,j,k,ip,jp,kp
logical :: isActive
!!$OMP PARALLEL DO default(shared) collapse(2) &
!!$OMP private(i,j,k,isActive,ip,jp,kp)
! this causes a race condition copying values between ghost cells
! probably a false positive since some of the cells are just copying
! values onto themselves, but easier to remove for now
!$OMP PARALLEL DO default(shared) collapse(2) &
!$OMP private(i,j,k,ip,jp,kp)
do k=Gr%ksg,Gr%keg
do j=Gr%jsg,Gr%jeg
do i=Gr%isg,Gr%ieg
isActive = (j >= Gr%js) .and. (j <= Gr%je) .and. &
(k >= Gr%ks) .and. (k <= Gr%ks)
if(.not. isActive) then
!If still here map this ghost to active and set value based on active
if (.not. isPhysical(Model,Gr,i,j,k)) then
!This is a geometric ghost so we can map it to a physical cell
call lfmIJKcc(Model,Gr,i,j,k,ip,jp,kp)
Q(i,j,k) = Q(ip,jp,kp)
endif
endif !isPhys
enddo
enddo !j
enddo !k
end subroutine FillGhostsCC
!Project xyz along dipole to R0 and return lat (x1) and lon (x2)
subroutine Proj2Rad(xyz,R0,x1,x2)
real(rp), intent(in ) :: xyz(NDIM), R0
!Checks if cell is "physical" as opposed to "geometric".
!Ie, i ghosts are physical but j/k ghosts are geometric (they point to other physical cells)
function isPhysical(Model,Gr,i,j,k)
type(Model_T), intent(in) :: Model
type(Grid_T) , intent(in) :: Gr
integer, intent(in) :: i,j,k
logical :: isPhysical
isPhysical = (j >= Gr%js) .and. (j <= Gr%je) .and. &
(k >= Gr%ks) .and. (k <= Gr%ke)
end function isPhysical
!Get azimuth and invariant latitude
subroutine NHProj(xyz,x1,x2)
real(rp), intent(in ) :: xyz(NDIM)
real(rp), intent(out) :: x1,x2
real(rp), dimension(NDIM) :: xyz0
xyz0 = DipoleShift(xyz,R0)
x1 = asin(xyz0(ZDIR)/R0) !Lat
x2 = katan2(xyz0(YDIR),xyz0(XDIR)) !katan => [0,2pi] instead of [-pi,pi]
end subroutine Proj2Rad
x1 = InvLatitude(xyz)
x2 = katan2(xyz(YDIR),xyz(XDIR)) !katan => [0,2pi] instead of [-pi,pi]
end subroutine NHProj
end subroutine

187
src/voltron/modelInterfaces/imag2mix_interface.F90
View File

@@ -10,59 +10,8 @@ module imag2mix_interface
implicit none
type(mixGrid_T), private :: rai2mixG
integer, private :: Np_mix, Nt_mix, Np_rai, Nt_rai, Npc_rai, Ntc_rai
logical, private :: isInit = .true.
contains
subroutine init_raiju_mix(imagApp,remixApp)
! called by subroutine initializeFromGamera in module voltapp
! type(voltApp_T), intent(inout) :: vApp
class(imagCoupler_T), intent(in) :: imagApp
type(mixApp_T), intent(inout) :: remixApp
real(rp), dimension(:,:), allocatable :: raijup, raijut ! Np x Nt, remix-style 2-D arrays to hold the RAIJU grid
integer :: i, j
! associate(remixApp=>vApp%remixApp, imagApp=>vApp%imagApp )
Nt_mix = remixApp%ion(NORTH)%G%Nt
Np_mix = remixApp%ion(NORTH)%G%Np
SELECT TYPE (imagA=>imagApp)
TYPE IS (raijuCoupler_T)
! in here you can treat imagType as if it is type raijuCoupler_T, and it points to vApp%imagApp
Np_rai = imagA%raiApp%Grid%shGrid%Np
Nt_rai = imagA%raiApp%Grid%shGrid%Nt
! Np x Nt, transposed relative to mix grid.
if (.not.allocated(raijup)) allocate(raijup(Np_rai, Nt_rai))
if (.not.allocated(raijut)) allocate(raijut(Np_rai, Nt_rai))
! construct the 2-D grid
!! thc/phc: (Nt or Np) [radians] grid centers
!! th (theta) is colatitude and runs from north pole toward south
do j=1,Np_rai
raijut(j,:) = imagA%raiApp%Grid%shGrid%thc
enddo
!! Phi is longitude, with zero/2pi at 12 MLT
do i=1,Nt_rai
raijup(:,i) = imagA%raiApp%Grid%shGrid%phc
enddo
! call remix grid constructor
call init_grid_fromTP(rai2mixG,raijut(1:Np_rai-1,:),raijup(1:Np_rai-1,:),isSolverGrid=.false.)
Npc_rai = rai2mixG%Np ! Np_rai-1
Ntc_rai = rai2mixG%Nt
CLASS DEFAULT
WRITE (*,*) "Imag Coupler is an unsupported type"
stop
END SELECT
! end associate
end subroutine init_raiju_mix
subroutine CoupleIMagToMix(vApp)
class(voltApp_T), intent(inout) :: vApp
integer :: i,j
@@ -265,141 +214,5 @@ module imag2mix_interface
end subroutine CoupleIMagToMix
subroutine mapRaijuToRemix(vApp)
type(voltApp_T), intent(inout) :: vApp
! type(raijuCoupler_T), intent(inout) :: imagApp
! type(mixApp_T), intent(inout) :: remixApp
real(rp), dimension(:,:,:), allocatable :: rai_fluxes, mix_fluxes
real(rp), dimension(:,:), allocatable :: mix_flux
real(rp), dimension(:,:), allocatable :: mixt, mixp
real(rp), dimension(:,:), allocatable :: raijup, raijut ! Np x Nt, remix-style 2-D arrays to hold the RAIJU grid
real(rp), dimension(:), allocatable :: phc, thc
integer :: Nf = nVars_imag2mix
integer :: i,j
type(Map_T) :: raiMap
! collect raiju fluxes.
! in getMomentsPrecip: allocate(rai_fluxes (is:ie,js:je,nVars_imag2mix)), (Nt_rai, Np_rai, Nf)
! call vApp%imagApp%getMomentsPrecip(rai_fluxes)
associate(remixApp=>vApp%remixApp ) !, imagApp=>vApp%imagApp
allocate(mix_fluxes(Np_mix,Nt_mix,Nf))
allocate(mix_flux(Np_mix,Nt_mix))
mix_fluxes = 0.0_rp
mix_flux = 0.0_rp
mixt = remixApp%ion(NORTH)%G%t ! G%t(G%Np,G%Nt)
mixp = remixApp%ion(NORTH)%G%p
! NH mapping to remix
call mix_set_map(rai2mixG,remixApp%ion(NORTH)%G,raiMap)
do i=1,Nf
call mix_map_grids(raiMap,transpose(rai_fluxes(:,1:Npc_rai,i)), mix_flux)
mix_fluxes(:,:,i) = mix_flux
enddo
remixApp%ion(NORTH)%St%Vars(:,:,IM_EAVG ) = mix_fluxes(:,:,RAI_EAVG ) ! [keV]
remixApp%ion(NORTH)%St%Vars(:,:,IM_ENFLX) = mix_fluxes(:,:,RAI_ENFLX) ! [#/cm^2/s]
remixApp%ion(NORTH)%St%Vars(:,:,IM_EFLUX) = mix_fluxes(:,:,RAI_EFLUX) ! [ergs/cm^2/s]
remixApp%ion(NORTH)%St%Vars(:,:,IM_GTYPE) = mix_fluxes(:,:,RAI_GTYPE) ! [0~1]
remixApp%ion(NORTH)%St%Vars(:,:,IM_EDEN ) = mix_fluxes(:,:,RAI_EDEN ) ! [#/m^3]
remixApp%ion(NORTH)%St%Vars(:,:,IM_EPRE ) = mix_fluxes(:,:,RAI_EPRE ) ! [Pa]
remixApp%ion(NORTH)%St%Vars(:,:,IM_NPSP ) = mix_fluxes(:,:,RAI_NPSP ) ! [#/m^3]
! SH mapping
mix_fluxes = 0.0_rp
call mapRaijuSToRemix(rai_fluxes,mixt,mixp,mix_fluxes)
remixApp%ion(SOUTH)%St%Vars(:,:,IM_EAVG ) = mix_fluxes(Np_mix:1:-1,:,RAI_EAVG ) ! [keV]
remixApp%ion(SOUTH)%St%Vars(:,:,IM_ENFLX) = mix_fluxes(Np_mix:1:-1,:,RAI_ENFLX) ! [#/cm^2/s]
remixApp%ion(SOUTH)%St%Vars(:,:,IM_EFLUX) = mix_fluxes(Np_mix:1:-1,:,RAI_EFLUX) ! [ergs/cm^2/s]
remixApp%ion(SOUTH)%St%Vars(:,:,IM_GTYPE) = mix_fluxes(Np_mix:1:-1,:,RAI_GTYPE) ! [0~1]
remixApp%ion(SOUTH)%St%Vars(:,:,IM_EDEN ) = mix_fluxes(Np_mix:1:-1,:,RAI_EDEN ) ! [#/m^3]
remixApp%ion(SOUTH)%St%Vars(:,:,IM_EPRE ) = mix_fluxes(Np_mix:1:-1,:,RAI_EPRE ) ! [Pa]
remixApp%ion(SOUTH)%St%Vars(:,:,IM_NPSP ) = mix_fluxes(Np_mix:1:-1,:,RAI_NPSP ) ! [#/m^3]
end associate
end subroutine mapRaijuToRemix
subroutine mapRaijuSToRemix(rai_fluxes,mixt,mixp,mix_fluxes)
! Directly map from irregular raiju SH grid to ReMIX.
real(rp), dimension(Nt_rai, Np_rai,nVars_imag2mix), intent(in) :: rai_fluxes
real(rp), dimension(Np_mix, Nt_mix), intent(in) :: mixt, mixp
real(rp), dimension(Np_mix, Nt_mix,nVars_imag2mix), intent(out) :: mix_fluxes
real(rp), dimension(Nt_rai, Np_rai) :: colatc, glongc
real(rp), dimension(:,:), allocatable :: mixtE, mixpE
real(rp), dimension(Np_mix, Nt_mix) :: Ainvdwgt2
real(rp) :: dlat, delt, delp, invdwgt
integer :: i, j, i0, j0, il, iu, jp, dj
! Source grid: latc is negative. colatc is positive from ~15 to 75 deg. Note latc=0 for open field lines.
colatc = PI-rai_fluxes(:,:,RAI_THCON) ! RAI_THCON is conjugate co-lat in radians pi/2-pi, PI - RAI_THCON is -> 0-pi/2
glongc = rai_fluxes(:,:,RAI_PHCON) ! conjugate long in radians, 0-2pi, need to double check if they are consistent at lon=0
! Destination grid: remix Grid.
dlat = mixt(1,2)-mixt(1,1)
! dj is the ratio of rcm dlon to remix dlon, i.e. min number of rcm/raiju cells to collect.
! now with 360 raiju/rcm longitudinal cells, dj is 1 with quad res remix grid.
dj = nint(dble(Np_mix)/dble(Np_rai))
! make an extended mixt/mixp grid for easier periodic boundary processing.
allocate(mixtE(1-dj:Np_mix+dj,1:Nt_mix))
allocate(mixpE(1-dj:Np_mix+dj,1:Nt_mix))
mixtE(1:Np_mix,:) = mixt
mixtE(1-dj:0,:) = mixt(Np_mix+1-dj:Np_mix,:)
mixtE(Np_mix+1:Np_mix+dj,:) = mixt(1:dj,:)
mixpE(1:Np_mix,:) = mixp
mixpE(1-dj:0,:) = mixp(Np_mix+1-dj:Np_mix,:)
mixpE(Np_mix+1:Np_mix+dj,:) = mixp(1:dj,:)
! Mapping: remix dlat is ~10x of rcm, dlon is ~1/3.6 of rcm. Remix lat is from 0-45 deg. RCM is from 15-75 deg.
! For each rcm SH point, find the nearest remix lat. If it's not too far away (within dlat) then
! find the nearest remix lon. Assign rcm contribution to the nearest lat shell within 2 rcm dlon.
! The difference is due to remix dlat is larger while dlon is smaller. Need to make sure all remix grids have some contribution from rcm.
! Lastly, normalize the contribution by total IDW.
Ainvdwgt2 = 0.0_rp
mix_fluxes = 0.0_rp
!$OMP PARALLEL DO default(shared) collapse(2) &
!$OMP private(i,j,i0,il,iu,j0,jp,delt,delp,invdwgt) &
!$OMP reduction(+:Ainvdwgt2,mix_fluxes)
do j=1,Np_rai
do i=1,Nt_rai
i0 = minloc(abs(mixt(1,:)-colatc(i,j)),1) ! Find the nearest remix colat index for rcm colatc(i,j)
if(mixt(1,i0)<=colatc(i,j)) then ! If the nearest remix colat is < rcm colatc, only collect rcm to this colat and its next grid.
il=i0
iu=min(i0+1,Nt_mix)
else ! Otherwise, collect from this point and its neighbor lat.
il=max(i0-1,1)
iu=i0
endif
do i0=il,iu
! For any remix grid, interpolate if rcm lat is within dlat away
if(abs(mixt(1,i0)-colatc(i,j))<dlat) then
jp = minloc(abs(mixp(:,1)-glongc(i,j)),1)
! 1 <= jp <= Np_mix
! jp-dj>= 1-dj; jp+dj<= Np_mix+dj
! mixtE/mixpE is from 1-dj:Np_mix+dj
do j0=jp-dj,jp+dj
delt = abs(mixtE(j0,i0)-colatc(i,j))
delp = abs((mixpE(j0,i0)-glongc(i,j)))*sin(mixtE(j0,i0))
invdwgt = 1.0_rp/sqrt(delt**2+delp**2)
mix_fluxes(j0,i0,:) = mix_fluxes(j0,i0,:) + rai_fluxes(i,j,:)*invdwgt
Ainvdwgt2(j0,i0) = Ainvdwgt2(j0,i0) + invdwgt
enddo
endif
enddo
! endif
enddo
enddo
!$OMP PARALLEL DO default(shared) collapse(2) &
!$OMP private(i0,j0)
do j0=1,Np_mix
do i0=1,Nt_mix
if(Ainvdwgt2(j0,i0)>0.0_rp) then
mix_fluxes(j0,i0,:) = mix_fluxes(j0,i0,:)/Ainvdwgt2(j0,i0)
endif
enddo
enddo
end subroutine mapRaijuSToRemix
end module

2
src/voltron/mpi/modelInterfaces/gamCouple_mpi_G2V.F90
View File

@@ -86,7 +86,7 @@ module gamCouple_mpi_G2V
! initialize F08 MPI objects
App%couplingComm = MPI_COMM_NULL
App%zeroArraytypes = (/ MPI_DATATYPE_NULL /)
App%zeroArraytypes(:) = (/ MPI_INTEGER /) ! = (/ MPI_DATATYPE_NULL /)
! split voltron helpers off of the communicator
! split couplingPoolComm into a communicator with only the non-helper voltron rank

8
src/voltron/mpi/modelInterfaces/gamCouple_mpi_V2G.F90
View File

@@ -559,13 +559,13 @@ module gamCouple_mpi_V2G
subroutine sendGameraCplDataMpi(gCplApp, CouplingTargetT)
class(gamCouplerMpi_volt_T), intent(inout) :: gCplApp
real(rp), intent(in) :: CouplingTargetT
real(rp), intent(inout) :: CouplingTargetT
call sendShallowCplDataMpi(gCplApp)
if(gCplApp%doDeep) call sendDeepCplDataMpi(gCplApp)
call sendCplTimeMpi(gCplApp, CouplingTargetT)
end subroutine
end subroutine sendGameraCplDataMpi
subroutine sendShallowCplDataMpi(gCplApp)
class(gamCouplerMpi_volt_T), intent(inout) :: gCplApp
@@ -628,14 +628,14 @@ module gamCouple_mpi_V2G
subroutine sendCplTimeMpi(gCplApp, CouplingTargetT)
class(gamCouplerMpi_volt_T), intent(inout) :: gCplApp
real(rp), intent(in) :: CouplingTargetT
real(rp), intent(inout) :: CouplingTargetT
integer :: ierr
! Send Target Time for next coupling
call mpi_bcast(CouplingTargetT,1,MPI_MYFLOAT, gCplApp%myRank, gCplApp%couplingComm, ierr)
end subroutine
end subroutine sendCplTimeMpi
end module

29
src/voltron/mpi/voltapp_mpi.F90
View File

@@ -391,6 +391,25 @@ module voltapp_mpi
call Toc("GCM2MIX")
end if
! tubes are only done after spinup
if(vApp%doDeep .and. vApp%time >= 0) then
if(vApp%useHelpers) call vhReqStep(vApp)
!Update i-shell to trace within in case rTrc has changed
vApp%iDeep = vApp%gApp%Grid%ie-1
!Pull in updated fields to CHIMP
call Tic("G2C")
call convertGameraToChimp(vApp%mhd2chmp,vApp%gApp,vApp%ebTrcApp)
call Toc("G2C")
if(vApp%useHelpers .and. vApp%doTubeHelp) then
call Tic("VoltTubes",.true.)
call VhReqTubeStart(vApp)
call Toc("VoltTubes",.true.)
endif
endif
! run remix
call Tic("ReMIX", .true.)
call runRemix(vApp)
@@ -405,19 +424,11 @@ module voltapp_mpi
! only do imag after spinup
if(vApp%doDeep .and. vApp%time >= 0) then
if(vApp%useHelpers) call vhReqStep(vApp)
! instead of PreDeep, use Tube Helpers and replicate other calls
!Update i-shell to trace within in case rTrc has changed
vApp%iDeep = vApp%gApp%Grid%ie-1
!Pull in updated fields to CHIMP
call Tic("G2C")
call convertGameraToChimp(vApp%mhd2chmp,vApp%gApp,vApp%ebTrcApp)
call Toc("G2C")
call Tic("VoltTubes",.true.)
if(vApp%useHelpers .and. vApp%doTubeHelp) then
call VhReqTubeStart(vApp)
! Tubes were started earlier
call vhReqTubeEnd(vApp)
! Now pack into tubeShell
call Tic("Tube2Shell")

15
src/voltron/mpi/volthelpers_mpi.F90
View File

@@ -107,7 +107,7 @@ module volthelpers_mpi
! chimp data update functions
subroutine sendChimpStateData(ebState, vHelpComm)
type(ebState_T), intent(in) :: ebState
type(ebState_T), intent(inout) :: ebState
type(MPI_Comm), intent(in) :: vHelpComm
integer :: ierr, length
@@ -195,7 +195,7 @@ module volthelpers_mpi
call mpi_Abort(MPI_COMM_WORLD, 1, ierr)
end if
end subroutine
end subroutine sendChimpStateData
subroutine recvChimpStateData(ebState, vHelpComm)
type(ebState_T), intent(inout) :: ebState
@@ -272,7 +272,7 @@ module volthelpers_mpi
end subroutine
subroutine sendChimpUpdate(vApp)
type(voltAppMpi_T), intent(in) :: vApp
type(voltAppMpi_T), intent(inout) :: vApp
integer :: ierr, length
character( len = MPI_MAX_ERROR_STRING) :: message
@@ -317,7 +317,7 @@ module volthelpers_mpi
call mpi_Abort(MPI_COMM_WORLD, 1, ierr)
end if
end subroutine
end subroutine sendChimpUpdate
subroutine recvChimpUpdate(vApp)
type(voltAppMpi_T), intent(inout) :: vApp
@@ -374,14 +374,15 @@ module volthelpers_mpi
type(voltAppMpi_T), intent(in) :: vApp
integer, intent(in) :: rType
integer :: ierr
integer :: ierr,wtf
type(MPI_Request) :: helpReq
wtf = rType
! async to match waiting helper nodes
call mpi_Ibcast(rType, 1, MPI_INTEGER, 0, vApp%vHelpComm, helpReq, ierr)
call mpi_Ibcast(wtf, 1, MPI_INTEGER, 0, vApp%vHelpComm, helpReq, ierr)
call mpi_wait(helpReq, MPI_STATUS_IGNORE, ierr)
end subroutine
end subroutine vhRequestType
subroutine vhReqStep(vApp)
type(voltAppMpi_T), intent(inout) :: vApp

15
src/voltron/voltapp.F90
View File

@@ -459,19 +459,8 @@ module voltapp
call init_volt2Chmp(vApp,gApp)
endif
!Ensure chimp and voltron restart numbers match
! Actually chimp doesn't write restart files right now
!if (isRestart .and. vApp%IO%nRes /= ebTrcApp%ebModel%IO%nRes) then
! write(*,*) "Voltron and Chimp disagree on restart number, you should sort that out."
! write(*,*) "Error code: A house divided cannot stand"
! write(*,*) " Voltron nRes = ", vApp%IO%nRes
! write(*,*) " Chimp nRes = ", ebTrcApp%ebModel%IO%nRes
! stop
!endif
call init_mhd2Chmp(vApp%mhd2chmp, gApp, vApp%ebTrcApp)
call init_chmp2Mhd(vApp%chmp2mhd, vApp%ebTrcApp, gApp)
call init_raiju_mix(vApp%imagApp,vApp%remixApp)
vApp%iDeep = gApp%Grid%ie-1
@@ -486,9 +475,7 @@ module voltapp
! convert gamera inputs to remix
call MJDRecalc(vApp%MJD)
if (vApp%doDeep) then
! call mapIMagToRemix(vApp%imag2mix,vApp%remixApp) ! rcm style
! call mapRaijuToRemix(vApp)
if ( vApp%doDeep .and. (vApp%time>0) ) then
call CoupleIMagToMix(vApp)
endif
call mapGameraToRemix(vApp%mhd2mix, vApp%remixApp)

2
testingScripts/build_mage-template.pbs
View File

@@ -9,8 +9,6 @@
#PBS -j oe
#PBS -m abe
# This script just builds the MAGE software.
echo "Job $PBS_JOBID started at `date` on `hostname` in directory `pwd`."
echo 'Loading modules.'

152
testingScripts/send_slack_message.py Normal file
View File

@@ -0,0 +1,152 @@
#!/usr/bin/env python
"""Send a message to Slack.
Send a message to Slack.
Authors
-------
Eric Winter
"""
# Import standard modules.
import datetime
# import glob
import os
import sys
# # Import 3rd-party modules.
# Import project modules.
import common
# Program constants
# Program description.
DESCRIPTION = "Send a message to Slack."
# # Root of directory tree for this set of tests.
# MAGE_TEST_SET_ROOT = os.environ['MAGE_TEST_SET_ROOT']
# # Directory for unit tests
# UNIT_TEST_DIRECTORY = os.path.join(MAGE_TEST_SET_ROOT, 'unitTest')
# # glob pattern for naming unit test directories
# UNIT_TEST_DIRECTORY_GLOB_PATTERN = 'unitTest_*'
# # Name of build subdirectory containing binaries
# BUILD_BIN_DIR = 'bin'
# # Name of file containing job IDs for each unit test directory.
# JOB_ID_LIST_FILE = 'jobs.txt'
def create_command_line_parser():
"""Create the command-line argument parser.
Create the parser for command-line arguments.
Parameters
----------
None
Returns
-------
parser : argparse.ArgumentParser
Command-line argument parser for this script.
Raises
------
None
"""
parser = common.create_command_line_parser(DESCRIPTION)
parser.add_argument(
"message",
default="",
help="Message to send to Slack (default: %(default)s)"
)
return parser
def send_slack_message(args: dict = None):
"""Send a message to Slack.
Send a message to Slack.
Parameters
----------
args : dict
Dictionary of program options.
Returns
-------
None
Raises
------
None
"""
# Local convenience variables.
debug = args["debug"]
be_loud = args["loud"]
slack_on_fail = args["slack_on_fail"]
is_test = args["test"]
verbose = args["verbose"]
message = args["message"]
# ------------------------------------------------------------------------
if debug:
print(f"Starting {sys.argv[0]} at {datetime.datetime.now()}")
print(f"Current directory is {os.getcwd()}")
# ------------------------------------------------------------------------
# Create the Slack client.
slack_client = common.slack_create_client()
slack_response_summary = common.slack_send_message(
slack_client, message, is_test=is_test
)
# ------------------------------------------------------------------------
if debug:
print(f"Ending {sys.argv[0]} at {datetime.datetime.now()}")
def main():
"""Begin main program.
This is the main program code.
Parameters
----------
None
Returns
-------
None
Raises
------
None
"""
# Set up the command-line parser.
parser = create_command_line_parser()
# Parse the command-line arguments.
args = parser.parse_args()
if args.debug:
print(f"args = {args}")
# ------------------------------------------------------------------------
# Call the main program logic. Note that the Namespace object (args)
# returned from the option parser is converted to a dict using vars().
send_slack_message(vars(args))
if __name__ == "__main__":
main()

44
testingScripts/unitTest-build-template.pbs Normal file
View File

@@ -0,0 +1,44 @@
#!/bin/bash
#PBS -N {{ job_name }}
#PBS -A {{ account }}
#PBS -q {{ queue }}
#PBS -l job_priority={{ job_priority }}
#PBS -l select=1:ncpus=128
#PBS -l walltime={{ walltime }}
#PBS -j oe
#PBS -m abe
# Abort script on any error.
set -e
echo "Job $PBS_JOBID started at `date` on `hostname` in directory `pwd`."
echo 'Loading modules.'
module --force purge
{%- for module in modules %}
module load {{ module }}
{%- endfor %}
echo 'The currently loaded modules are:'
module list
echo 'The active environment variables are:'
printenv
echo 'Copying pFUnit binaries.'
pfunit_dir="{{ mage_test_root }}/pfunit/pFUnit-4.2.0/ifort-23-mpich-derecho"
kaiju_external_dir="{{ kaijuhome }}/external"
cp -rp "${pfunit_dir}/FARGPARSE-1.1" "${kaiju_external_dir}/"
cp -rp "${pfunit_dir}/GFTL-1.3" "${kaiju_external_dir}/"
cp -rp "${pfunit_dir}/GFTL_SHARED-1.2" "${kaiju_external_dir}/"
cp -rp "${pfunit_dir}/PFUNIT-4.2" "${kaiju_external_dir}/"
# Build the code.
cmd="{{ cmake_cmd }} >& cmake.out"
echo $cmd
eval $cmd
cmd="{{ make_cmd }} >& make.out"
echo $cmd
eval $cmd
echo "Job $PBS_JOBID ended at `date` on `hostname` in directory `pwd`."

1276
testingScripts/unitTest.py
View File

File diff suppressed because it is too large Load Diff

39
testingScripts/unitTestReport.py
View File

@@ -111,8 +111,7 @@ def main():
print(f"Checking unit test results in {unit_test_directory}.")
# Move to the directory containing the unit test results.
path = os.path.join(UNIT_TEST_DIRECTORY, unit_test_directory,
BUILD_BIN_DIR)
path = os.path.join(UNIT_TEST_DIRECTORY, unit_test_directory)
if debug:
print(f"path = {path}")
os.chdir(path)
@@ -136,19 +135,27 @@ def main():
# NOTE: This needs to be reorganized.
# Compute the names of the job log files.
job_file_0 = f"genTestData.o{job_ids[0]}" # 0 OKs
job_file_1 = f"runCaseTests.o{job_ids[1]}" # 2 OKs
job_file_2 = f"runNonCaseTests1.o{job_ids[2]}" # 6 OKs
job_file_3 = f"runNonCaseTests2.o{job_ids[3]}" # 1 OK
if debug:
print(f"job_file_0 = {job_file_0}")
print(f"job_file_1 = {job_file_1}")
print(f"job_file_2 = {job_file_2}")
print(f"job_file_3 = {job_file_3}")
# 0 OKs
job_file_build = f"../unitTest-build.o{job_ids[0]}"
# 0 OKs
job_file_genTestData = f"../unitTest-genTestData.o{job_ids[1]}"
# 2 OKs
job_file_caseTests = f"../unitTest-caseTests.o{job_ids[2]}"
# 6 OKs
job_file_noncaseTests1 = f"../unitTest-noncaseTests1.o{job_ids[3]}"
# 1 OK
job_file_noncaseTests2 = f"../unitTest-noncaseTests2.o{job_ids[4]}"
# Combine the results of each test log file.
os.chdir("bin")
bigFile = []
job_files = [job_file_0, job_file_1, job_file_2, job_file_3]
job_files = [
job_file_build,
job_file_genTestData,
job_file_caseTests,
job_file_noncaseTests1,
job_file_noncaseTests2,
]
for job_file in job_files:
with open(job_file, 'r', encoding='utf-8') as f:
bigFile += f.readlines()
@@ -234,12 +241,14 @@ def main():
)
if debug:
print(f"slack_response_summary = {slack_response_summary}")
# Also write a summary file to the root folder of this test
with open(os.path.join(MAGE_TEST_SET_ROOT,'testSummary.out'), 'w', encoding='utf-8') as f:
with open(os.path.join(
MAGE_TEST_SET_ROOT, 'testSummary.out'), 'w', encoding='utf-8'
) as f:
f.write(test_report_details_string)
f.write('\n')
# ------------------------------------------------------------------------
if debug:

27
testingScripts/weeklyDash-template.pbs
View File

@@ -20,19 +20,26 @@ echo 'The currently loaded modules are:'
module list
echo 'Loading python environment.'
mage_miniconda3="${HOME}/miniconda3"
mage_conda="${mage_miniconda3}/bin/conda"
__conda_setup="$($mage_conda 'shell.bash' 'hook' 2> /dev/null)"
if [ $? -eq 0 ]; then
eval "$__conda_setup"
else
if [ -f "$mage_miniconda3/etc/profile.d/conda.sh" ]; then
. "$mage_miniconda3/etc/profile.d/conda.sh"
mage_test_root=$HOME
if [ -d "${mage_test_root}/miniconda3" ]; then
echo 'Loading local miniconda3'
mage_miniconda3="${mage_test_root}/miniconda3"
mage_conda="${mage_miniconda3}/bin/conda"
__conda_setup="$($mage_conda 'shell.bash' 'hook' 2> /dev/null)"
if [ $? -eq 0 ]; then
eval "$__conda_setup"
else
export PATH="$mage_miniconda3/bin:$PATH"
if [ -f "$mage_miniconda3/etc/profile.d/conda.sh" ]; then
. "$mage_miniconda3/etc/profile.d/conda.sh"
else
export PATH="$mage_miniconda3/bin:$PATH"
fi
fi
unset __conda_setup
else
echo 'Loading conda module'
module load conda
fi
unset __conda_setup
conda activate {{ conda_environment }}
echo "The current conda environment is ${CONDA_PREFIX}."

13
tests/genTestData-template.pbs
View File

@@ -9,6 +9,9 @@
#PBS -j oe
#PBS -m abe
# Abort script on any error.
set -e
echo "Job $PBS_JOBID started at `date` on `hostname` in directory `pwd`."
echo 'Loading modules.'
@@ -28,6 +31,16 @@ export KMP_STACKSIZE=128M
echo 'The active environment variables are:'
printenv
# Move to the directory containing the compiled code.
cd bin
echo 'Copying input files.'
test_inputs_dir="{{ mage_test_root }}/unit_test_inputs"
cp "${test_inputs_dir}/bcwind.h5" .
cp "${test_inputs_dir}/geo_mpi.xml" .
cp "${test_inputs_dir}/lfmD.h5" .
cp "${test_inputs_dir}/raijuconfig.h5" .
echo 'Generating data for testing.'
MPICOMMAND="mpiexec $KAIJUHOME/scripts/preproc/pinCpuCores.sh"
$MPICOMMAND ./voltron_mpi.x cmiD_deep_8_genRes.xml >& cmiD_deep_8_genRes.out

6
tests/runCaseTests-template.pbs
View File

@@ -9,6 +9,9 @@
#PBS -j oe
#PBS -m abe
# Abort script on any error.
set -e
echo "Job $PBS_JOBID started at `date` on `hostname` in directory `pwd`."
echo 'Loading modules.'
@@ -28,6 +31,9 @@ export KMP_STACKSIZE=128M
echo 'The active environment variables are:'
printenv
# Move to the directory containing the compiled code.
cd bin
echo 'Running non-MPI test cases.'
./caseTests >& caseTests.out
echo 'Non-MPI test cases complete.'

6
tests/runNonCaseTests1-template.pbs
View File

@@ -9,6 +9,9 @@
#PBS -j oe
#PBS -m abe
# Abort script on any error.
set -e
echo "Job $PBS_JOBID started at `date` on `hostname` in directory `pwd`."
echo 'Loading modules.'
@@ -28,6 +31,9 @@ export KMP_STACKSIZE=128M
echo 'The active environment variables are:'
printenv
# Move to the directory containing the compiled code.
cd bin
echo 'Running GAMERA tests.'
date
./gamTests >& gamTests.out

6
tests/runNonCaseTests2-template.pbs
View File

@@ -9,6 +9,9 @@
#PBS -j oe
#PBS -m abe
# Abort script on any error.
set -e
echo "Job $PBS_JOBID started at `date` on `hostname` in directory `pwd`."
echo 'Loading modules.'
@@ -28,6 +31,9 @@ export KMP_STACKSIZE=128M
echo 'The active environment variables are:'
printenv
# Move to the directory containing the compiled code.
cd bin
echo 'Running VOLTRON MPI tests.'
date
MPICOMMAND="mpiexec $KAIJUHOME/scripts/preproc/pinCpuCores.sh"

11
tests/unitTestReport-template.pbs
View File

@@ -9,6 +9,9 @@
#PBS -j oe
#PBS -m abe
# Abort script on any error.
set -e
echo "Job $PBS_JOBID started at `date` on `hostname` in directory `pwd`."
echo 'Loading modules.'
@@ -40,6 +43,7 @@ else
module load conda
fi
conda activate {{ conda_environment }}
echo "The active conda environment is ${CONDA_DEFAULT_ENV}."
echo 'Setting up MAGE environment.'
source {{ kaijuhome }}/scripts/setupEnvironment.sh
@@ -51,6 +55,13 @@ export BRANCH_OR_COMMIT={{ branch_or_commit }}
echo 'The active environment variables are:'
printenv
# Move to the directory containing the compiled code.
cd bin
if [[ $? -eq 1 ]]; then
python $KAIJUHOME/testingScripts/send_slack_message.py "Unit test build failed in `pwd`!"
exit 1
fi
echo 'Generating unit test report.'
python $KAIJUHOME/testingScripts/unitTestReport.py {{ report_options }} >& unitTestReport.out