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

Maybe working parent2child sgv interp. unit test ready but can't compile on my laptop :). Raiju getting total and corot pot from voltron

Browse Source
This commit is contained in:
Anthony M. Sciola
2025-02-11 15:09:23 -05:00
parent 5b88898e63
commit f3fb88dad5
11 changed files with 176 additions and 57 deletions
Download Patch File Download Diff File Expand all files Collapse all files

26
src/base/shellutils/shellGrid.F90
View File

@@ -13,6 +13,8 @@ module shellGrid
enum, bind(C)
enumerator :: SHGR_CC,SHGR_CORNER,SHGR_FACE_THETA,SHGR_FACE_PHI
end enum
character(len=strLen), parameter :: def_noParentName ="N/A"
type ShellGridVar_T
@@ -89,8 +91,8 @@ end type ShellGridVar_T
!> Subgrid information
!> ShellGrids that are subgrids of other shellGrids store info about their parent grid
logical :: isChild = .false.
character(len=strLen) :: parentName = "N/A"
!! Name of the parent grid that this one derives from
character(len=strLen) :: parentName = def_noParentName
!! Name of the parent grid that this one derives from
integer :: bndis,bndie,bndjs,bndje
!! Indices of parent grid that bound this grid
@@ -475,10 +477,10 @@ end type ShellGridVar_T
!! Actual bounds used
! If a bound is provided then we use that, if not we default to parent grid's bounds
is = whichInd(sub_is, pSG%is, present(sub_is))
ie = whichInd(sub_ie, pSG%ie+1, present(sub_ie))
js = whichInd(sub_js, pSG%js , present(sub_js))
je = whichInd(sub_je, pSG%je+1, present(sub_je))
is = whichInd(sub_is, pSG%is, present(sub_is))
ie = whichInd(sub_ie, pSG%ie, present(sub_ie))
js = whichInd(sub_js, pSG%js, present(sub_js))
je = whichInd(sub_je, pSG%je, present(sub_je))
! Check for valid bounds
if (is < pSG%is) then
@@ -487,10 +489,10 @@ end type ShellGridVar_T
write(*,*) "Mimumum:",pSG%is
stop
endif
if (ie > pSG%ie+1) then
if (ie > pSG%ie) then
write(*,*) "ERROR GenChildShellGrid: Invalid ie bound."
write(*,*) "Requested:",ie
write(*,*) "Maximum:",pSG%ie+1
write(*,*) "Maximum:",pSG%ie
stop
endif
if (js < pSG%js) then
@@ -499,10 +501,10 @@ end type ShellGridVar_T
write(*,*) "Mimumum:",pSG%js
stop
endif
if (je > pSG%je+1) then
if (je > pSG%je) then
write(*,*) "ERROR GenChildShellGrid: Invalid je bound."
write(*,*) "Requested:",je
write(*,*) "Maximum:",pSG%je+1
write(*,*) "Maximum:",pSG%je
stop
endif
if (is > ie) then
@@ -539,11 +541,11 @@ end type ShellGridVar_T
endif
! Otherwise, this ghost definition is okay
call GenShellGrid(cSG, pSG%th(is:ie), pSG%ph(js:je), name, nGhosts=nGhosts, radO=pSG%radius)
call GenShellGrid(cSG, pSG%th(is:ie+1), pSG%ph(js:je+1), name, nGhosts=nGhosts, radO=pSG%radius)
else
! No ghosts defined, go with default
call GenShellGrid(cSG, pSG%th(is:ie), pSG%ph(js:je), name, radO=pSG%radius)
call GenShellGrid(cSG, pSG%th(is:ie+1), pSG%ph(js:je+1), name, radO=pSG%radius)
endif

47
src/base/shellutils/shellInterp.F90
View File

@@ -17,8 +17,6 @@ module shellInterp
! InterpShellVar - TODO
! InterpShellVar_ChildToParent - TODO
! InterpShellVar_ParentToChild - TODO
! InterpShellVar_TSC_SG
! InterpShellVar_TSC_pnt
subroutine InterpShellVar(sgSource, sgVar, sgDest, varOut)
!! This is meant to be the highest-abstraction option for interpolation
@@ -184,7 +182,7 @@ module shellInterp
end subroutine InterpShellVar_TSC_SG
subroutine InterpShellVar_TSC_pnt(sgsource, sgVar, th, pin, Qinterp, dThetaO, dPhiO, srcMaskO, goodInterpO)
subroutine InterpShellVar_TSC_pnt(sgSource, sgVar, th, pin, Qinterp, dThetaO, dPhiO, srcMaskO, goodInterpO)
!! Given the source information, interpolate sgVar to point (t,pin) and return as Qout
type(ShellGrid_T ), intent(in) :: sgSource
!! Source ShellGrid
@@ -265,7 +263,7 @@ module shellInterp
if (present(dThetaO)) then
! First, make sure dTheta and dPhi are defined at sgVar locations
if ( size(dThetaO) .ne. sgVar%Ni ) then
write(*,*)"ERROR in InterpShellVar_TSC_pnt: dTheta != sgVar%Ni"
write(*,*)"ERROR in InterpShellVar_TSC_pnt: size(dTheta) != sgVar%Ni"
write(*,*) size(dThetaO), sgVar%Ni
stop
endif
@@ -281,7 +279,7 @@ module shellInterp
if (present(dPhiO)) then
if ( size(dPhiO) .ne. sgVar%Nj ) then
write(*,*)"ERROR in InterpShellVar_TSC_pnt: dPhi != sgVar%Nj"
write(*,*)"ERROR in InterpShellVar_TSC_pnt: size(dPhi) != sgVar%Nj"
write(*,*) size(dPhiO), sgVar%Nj
stop
endif
@@ -397,6 +395,45 @@ module shellInterp
end subroutine InterpShellVar_TSC_pnt
subroutine InterpShellVar_ParentToChild(srcGrid, srcVar, destGrid, destVar)
type(ShellGrid_T ), intent(in) :: srcGrid
type(ShellGridVar_T), intent(in) :: srcVar
type(ShellGrid_T ), intent(in) :: destGrid
type(ShellGridVar_T), intent(inout) :: destVar
integer :: subis, subie, subjs, subje, iExtra, jExtra
type(ShellGridVar_T) :: src2destVar
!! srcVar chopped down to dest grid's domain
! Are you actually my parent or do I need to call the police?
if (trim(destGrid%parentName) == def_noParentName .or. trim(srcGrid%name) .ne. trim(destGrid%parentName)) then
write(*,*) " ERROR in InterpShellVar_ParentToChild: Src grid not parent of Dest grid"
write(*,*) " Src: ", trim(srcGrid%name), ", Dest: ",trim(destGrid%parentName)
stop
endif
! Determine SGV bounds if source var that map to destination grid's bounds
call ijExtra_SGV(srcVar%loc, iExtra, jExtra)
subis = srcVar%isv + destGrid%bndis-1 - destGrid%Ngn
subie = srcVar%isv + destGrid%bndie-1 + destGrid%Ngs + iExtra
subjs = srcVar%jsv + destGrid%bndjs-1 - destGrid%Ngw
subje = srcVar%jsv + destGrid%bndje-1 + destGrid%Nge + jExtra
! Now fill in destination variable
if (srcVar%loc == destVar%loc) then
destVar%data = srcVar%data(subis:subie, subjs:subje)
destVar%mask = srcVar%mask(subis:subie, subjs:subje)
else ! Need to interpolate
! Chop srcVar into destination grid size
call initShellVar(destGrid, srcVar%loc, src2destVar)
src2destVar%data = srcVar%data(subis:subie, subjs:subje)
src2destVar%mask = srcVar%mask(subis:subie, subjs:subje)
call InterpShellVar_TSC_SG(destGrid, src2destVar, destGrid, destVar)
endif
end subroutine InterpShellVar_ParentToChild
subroutine interpPole(shGr,Qin,tin,pin,Qinterp)
type(ShellGrid_T), intent(in) :: shGr ! source grid
type(ShellGridVar_T), intent(in) :: Qin ! source grid variable

25
src/base/shellutils/shellUtils.F90
View File

@@ -8,6 +8,31 @@ module shellUtils
contains
subroutine ijExtra_SGV(loc, iExtra, jExtra)
!! Determine which dimensions have extra index relative to # cells based on variable's location on grid
integer, intent(in) :: loc
integer, intent(out) :: iExtra
integer, intent(out) :: jExtra
select case(loc)
case(SHGR_CC)
iExtra = 0
jExtra = 0
case(SHGR_CORNER)
iExtra = 1
jExtra = 1
case(SHGR_FACE_THETA)
iExtra = 1
jExtra = 0
case(SHGR_FACE_PHI)
iExtra = 0
jExtra = 1
case default
write(*,*) "initShellGridVar got an invalid data location:",loc
stop
end select
end subroutine ijExtra_SGV
subroutine wrapJ_SGV(sh, sgVar)
!! Wrap a ShellGridVar around the periodic j/phi boundary
!! Assumes all vars within active domain are valid to wrap

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

@@ -408,6 +408,8 @@ module raijutypes
! -- Coming from ionosphere solve -- !
real(rp), dimension(:,:), allocatable :: espot
!! (Ngi+1, Ngj+1) [kV] electro-static potential
real(rp), dimension(:,:), allocatable :: pot_corot
!! (Ngi+1, Ngj+1) [kV] corotation potential
integer, dimension(:), allocatable :: bndLoc
!! (Ngi) i value of boundary between inactive and buffer domain

6
src/base/types/volttypes.F90
View File

@@ -204,8 +204,10 @@ module volttypes
!! (Ngi, Ngj) Bounce timesale
type(ShellGridVar_T) :: pot
!! electrostatic potential from ionosphere [kV]
type(ShellGridVar_T) :: pot_total
!! Total electrostatic potential from (ionosphere + corot) [kV]
type(ShellGridVar_T) :: pot_corot
!! Just corotation potential, just for output purposes [kV]
contains

6
src/raiju/raijuGrids.F90
View File

@@ -254,8 +254,8 @@ module raijugrids
end subroutine raijuGenWarpSphGrid_Shafee2008
subroutine raijuGenGridFromShGrid(Grid, shGrid, iXML)
type(raijuGrid_T) , intent(inout) :: Grid
subroutine raijuGenGridFromShGrid(raijuGrid, shGrid, iXML)
type(ShellGrid_T) , intent(inout) :: raijuGrid
type(ShellGrid_T), intent(in) :: shGrid
type(XML_Input_T), intent(in) :: iXML
@@ -304,7 +304,7 @@ module raijugrids
endif
! Now we can make our grid
call GenChildShellGrid(shGrid, Grid%shGrid, RAI_SG_NAME, nGhosts, sub_is=iStart, sub_ie=iEnd)
call GenChildShellGrid(shGrid, raijuGrid, RAI_SG_NAME, nGhosts, sub_is=iStart, sub_ie=iEnd)
end associate

10
src/raiju/raijuIO.F90
View File

@@ -255,7 +255,6 @@ module raijuIO
call AddOutVar(IOVars,"lonc" ,State%phcon (is:ie+1,js:je+1) ,uStr="radians",dStr="(corners) Congugate longitude")
call AddOutVar(IOVars,"active" ,State%active (is:ie,js:je)*1.0_rp ,uStr="-1=Inactive, 0=Buffer, 1=Active")
call AddOutVar(IOVars,"OCBDist",State%OCBDist(is:ie,js:je)*1.0_rp ,dStr="Cell distance from an open closed boundary")
call AddOutVar(IOVars,"espot" ,State%espot (is:ie+1,js:je+1) ,uStr="kV",dStr="(corners) Electrostatic potential")
call AddOutVar(IOVars,"bVol" ,State%bvol (is:ie+1,js:je+1) ,uStr="Rx/nT",dStr="(corners) Flux Tube Volume")
call AddOutVar(IOVars,"bVol_cc",State%bvol_cc(is:ie+1,js:je+1) ,uStr="Rx/nT",dStr="(corners) Flux Tube Volume")
call AddOutVar(IOVars,"vaFrac" ,State%vaFrac (is:ie+1,js:je+1) ,uStr="fraction",dStr="Fraction of Alfven speed over magnetofast + ExB speed")
@@ -265,6 +264,15 @@ module raijuIO
call AddOutVar(IOVars,"Dstd_in",State%Dstd (is:ie,js:je, :) ,uStr="normalized" ,dStr="Std. dev. of species density from imagtubes")
call AddOutSGV(IOVars,"Tbounce",State%Tb, outBndsO=outBnds2D, uStr="[s]", dStr="Bounce timescale along field line (Alfven crossing time)", doWriteMaskO=.false.)
if (Model%doOwnCorot) then
! If we are addng our own corotation potential, then its not included in espot
call AddOutVar(IOVars,"pot_iono",State%espot(is:ie+1,js:je+1) ,uStr="kV",dStr="(corners) Ionospheric potential")
else
! Otherwise, its in espot, so substract it out
call AddOutVar(IOVars,"pot_iono",State%espot(is:ie+1,js:je+1)-State%pot_corot(is:ie+1,js:je+1) ,uStr="kV",dStr="(corners) Ionospheric potential")
endif
call AddOutVar(IOVars,"pot_corot",State%pot_corot(is:ie+1,js:je+1),uStr="kV",dStr="(corners) Corotation potential")
! Idk about you but I did not expect true to equal -1
allocate(outTmp2D(is:ie, Grid%Nk))
where (State%activeShells)

13
src/raiju/raijuStarter.F90
View File

@@ -322,7 +322,7 @@ module raijustarter
! Then we should be receiving a predefined ShellGrid that Voltron has set up
if(present(shGridO)) then
shGrid = shGridO
call raijuGenGridFromShGrid(Grid, shGrid, iXML)
call raijuGenGridFromShGrid(Grid%shGrid, shGrid, iXML)
else
write(*,*) "RAIJU expecting a ShellGrid_T but didn't receive one. Dying."
endif
@@ -409,11 +409,12 @@ module raijustarter
allocate( State%thcon (sh%isg:sh%ieg+1, sh%jsg:sh%jeg+1 ) )
allocate( State%phcon (sh%isg:sh%ieg+1, sh%jsg:sh%jeg+1 ) )
! 2D corner quantities
allocate( State%topo (sh%isg:sh%ieg+1, sh%jsg:sh%jeg+1) )
allocate( State%espot (sh%isg:sh%ieg+1, sh%jsg:sh%jeg+1) )
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) )
allocate( State%topo (sh%isg:sh%ieg+1, sh%jsg:sh%jeg+1) )
allocate( State%espot (sh%isg:sh%ieg+1, sh%jsg:sh%jeg+1) )
allocate( State%pot_corot(sh%isg:sh%ieg+1, sh%jsg:sh%jeg+1) )
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

2
src/voltron/modelInterfaces/raiCplTypesSub.F90
View File

@@ -31,7 +31,7 @@ submodule (volttypes) raijuCplTypesSub
App%raiApp%State%isFirstCpl = .false.
endif
! Then allocate and initialize coupling variables based on raiju app
call raijuCpl_init(App)
call raijuCpl_init(App, xml)
end subroutine raiCplInitModel

64
src/voltron/modelInterfaces/raijuCplHelper.F90
View File

@@ -5,17 +5,19 @@ module raijuCplHelper
use remixReader
use shellinterp
use ebtypes
use raijugrids
use imagtubes
use mixdefs
implicit none
contains
subroutine raijuCpl_init(raiCpl)
subroutine raijuCpl_init(raiCpl, iXML)
class(raijuCoupler_T), intent(inout) :: raiCpl
type(XML_Input_T), intent(inout) :: iXML
integer, dimension(4) :: shGhosts
integer :: i
@@ -32,10 +34,13 @@ module raijuCplHelper
shGhosts(SOUTH) = sh%Ngs
shGhosts(EAST) = sh%Nge
shGhosts(WEST) = sh%Ngw
call GenChildShellGrid(sh, raiCpl%shGr, "raijuCpl", nGhosts=shGhosts)
call initShellVar(raiCpl%shGr, SHGR_CORNER, raiCpl%pot)
call initShellVar(sh, SHGR_CC, raiCpl%bvol_cc)
associate(mask_corner=>raiCpl%pot%mask, mask_cc=>raiCpl%bvol_cc%mask)
!call GenChildShellGrid(sh, raiCpl%shGr, "raijuCpl", nGhosts=shGhosts)
!call GenChildShellGrid(sh, raiCpl%opt%voltGrid, "raijuCpl", nGhosts=shGhosts)
call raijuGenGridFromShGrid(raiCpl%shGr, raiCpl%opt%voltGrid, iXML)
call initShellVar(raiCpl%shGr, SHGR_CORNER, raiCpl%pot_total)
call initShellVar(raiCpl%shGr, SHGR_CORNER, raiCpl%pot_corot)
call initShellVar(raiCpl%shGr, SHGR_CC, raiCpl%bvol_cc)
associate(mask_corner=>raiCpl%pot_total%mask, mask_cc=>raiCpl%bvol_cc%mask)
mask_corner = .true.
mask_cc = .true.
@@ -66,8 +71,10 @@ module raijuCplHelper
! Initial values
raiCpl%tLastUpdate = -1.0*HUGE
raiCpl%pot%data = 0.0
raiCpl%pot%mask = .true.
raiCpl%pot_total%data = 0.0
raiCpl%pot_total%mask = .true.
raiCpl%pot_corot%data = 0.0
raiCpl%pot_corot%mask = .true.
end subroutine raijuCpl_init
@@ -245,7 +252,8 @@ module raijuCplHelper
associate(Model=>raiCpl%raiApp%Model, State=>raiCpl%raiApp%State, shGr=>raiCpl%shGr)
!--- Ionosphere data ---!
State%espot(:,:) = raiCpl%pot%data(:,:) ! They live on the same grid so this is okay
State%espot(:,:) = raiCpl%pot_total%data(:,:) ! They live on the same grid so this is okay
State%pot_corot(:,:) = raiCpl%pot_corot%data(:,:)
!--- Mag data ---!
@@ -307,7 +315,7 @@ module raijuCplHelper
! Using chimp, populate imagTubes
call genImagTubes(raiCpl, vApp)
! Set potential
call InterpShellVar_TSC_SG(rmReader%shGr, rmReader%nsPot(1), raiCpl%shGr, raiCpl%pot)
call InterpShellVar_TSC_SG(rmReader%shGr, rmReader%nsPot(1), raiCpl%shGr, raiCpl%pot_total)
end subroutine packRaijuCoupler_OWD
@@ -367,25 +375,27 @@ module raijuCplHelper
!call genImagTubes(raiCpl, vApp)
call tubeShell2RaiCpl(vApp%shGrid, vApp%State%tubeShell, raiCpl)
!call mixPot2Raiju_RT(raiCpl, vApp%remixApp)
call InterpShellVar_TSC_SG(vApp%shGrid, vApp%State%potential_total, raiCpl%shGr, raiCpl%pot)
call InterpShellVar_ParentToChild(vApp%shGrid, vApp%State%potential_total, raiCpl%shGr, raiCpl%pot_total)
call InterpShellVar_ParentToChild(vApp%shGrid, vApp%State%potential_corot, raiCpl%shGr, raiCpl%pot_corot)
end subroutine
subroutine mixPot2Raiju_RT(raiCpl, rmApp)
!! Take remix's potential, shove it into a remix ShellGrid, use InterpShellVar to get it onto raiju's ShellGrid
class(raijuCoupler_T), intent(inout) :: raiCpl
class(mixApp_T), intent(inout) :: rmApp
real(rp), dimension(rmApp%ion(NORTH)%shGr%Nt,rmApp%ion(NORTH)%shGr%Np) :: tmpPot
associate(rmHemi=>rmApp%ion(NORTH), Nt=>rmApp%ion(NORTH)%shGr%Nt, Np=>rmApp%ion(NORTH)%shGr%Np)
rmHemi%St%pot_shGr%data(:,1:Np) = transpose(rmHemi%St%Vars(:,:,POT))
rmHemi%St%pot_shGr%data(:,Np+1) = rmHemi%St%pot_shGr%data(:,1)
rmHemi%St%pot_shGr%mask = .true.
call InterpShellVar_TSC_SG(rmHemi%shGr, rmHemi%St%pot_shGr, raiCpl%shGr, raiCpl%pot)
end associate
end subroutine mixPot2Raiju_RT
! subroutine mixPot2Raiju_RT(raiCpl, rmApp)
! !! Take remix's potential, shove it into a remix ShellGrid, use InterpShellVar to get it onto raiju's ShellGrid
! class(raijuCoupler_T), intent(inout) :: raiCpl
! class(mixApp_T), intent(inout) :: rmApp
!
! real(rp), dimension(rmApp%ion(NORTH)%shGr%Nt,rmApp%ion(NORTH)%shGr%Np) :: tmpPot
!
! associate(rmHemi=>rmApp%ion(NORTH), Nt=>rmApp%ion(NORTH)%shGr%Nt, Np=>rmApp%ion(NORTH)%shGr%Np)
! rmHemi%St%pot_shGr%data(:,1:Np) = transpose(rmHemi%St%Vars(:,:,POT))
! rmHemi%St%pot_shGr%data(:,Np+1) = rmHemi%St%pot_shGr%data(:,1)
! rmHemi%St%pot_shGr%mask = .true.
! call InterpShellVar_TSC_SG(rmHemi%shGr, rmHemi%St%pot_shGr, raiCpl%shGr, raiCpl%pot)
! end associate
!
! end subroutine mixPot2Raiju_RT
end module raijuCplHelper

32
tests/shellgrid/testSGInterp.pf
View File

@@ -9,9 +9,18 @@ module testSGInterp
implicit none
integer, parameter :: Npnt = 50
integer, parameter :: sub_is = 13
integer, parameter :: sub_ie = 17
!integer, parameter :: sub_ie = 37
integer, parameter :: sub_js = 17
integer, parameter :: sub_je = 21
!integer, parameter :: sub_je = 41
type(ShellGrid_T) :: sh
type(ShellGridVar_T) :: shVar_rand_corner
character(len=strLen) :: checkMessage
integer, parameter :: rand_seed = 86456
contains
@before
@@ -33,6 +42,12 @@ module testSGInterp
enddo
call GenShellGrid(sh, thetas, phis,'shgrTest_loc', radO=1.1_rp)
call initShellVar(sh, SHGR_CORNER, shVar_rand_corner)
do j=shVar_rand_corner%jsv:shVar_rand_corner%jev
do i=shVar_rand_corner%isv:shVar_rand_corner%iev
shVar_rand_corner%data(i,j) = rand(rand_seed)
enddo
enddo
end subroutine setup
@@ -143,4 +158,21 @@ module testSGInterp
continue
end subroutine testEquivDCell_Corner
@test
subroutine testParentToChild_sameLoc()
type(ShellGrid_T ) :: destGrid
type(ShellGridVar_T) :: destVar ! Contains source data
call GenChildShellGrid(sh, destGrid, "testChild", sub_is=sub_is,sub_ie=sub_ie,sub_js=sub_js,sub_je=sub_je)
call initShellVar(destGrid, shVar_rand_corner%loc, destVar)
call InterpShellVar_ParentToChild(sh, shVar_rand_corner, destGrid, destVar)
write(*,*)shVar_rand_corner%data(sub_is:sub_ie,sub_js:sub_je)
write(*,*)'-----'
write(*,*)destVar%data
end subroutine testParentToChild_sameLoc
end module testSGInterp