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

Realtime one-way maybe good? Still need coldStart

Browse Source
This commit is contained in:
Anthony
2024-09-27 15:44:11 -06:00
parent 1445f27548
commit 440321f9c1
15 changed files with 87 additions and 74 deletions
Download Patch File Download Diff File Expand all files Collapse all files

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

@@ -210,6 +210,8 @@ module raijutypes
real(rp) :: tiote ! Ion temp over electron temp. In the future, should be fancier
! Active domain settings
integer :: n_bndLim
!! cell dt/dp slope
real(rp) :: maxTail_buffer
!! Maximum tailward extent of the buffer region
real(rp) :: maxSun_buffer

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

@@ -110,8 +110,10 @@ module volttypes
type, extends(BaseOptions_T) :: imagOptions_T
character(len=strLen) :: swF
!! Solar wind filename
real(rp) :: lowLatBC
!! GAMERA low latitude boundary
real(rp) :: mhd_Rin
!! GAMERA near-Earth boundary
real(rp) :: mjd0
!! MJD at sim t=0
contains
end type imagOptions_T

5
src/chimp/streamline.F90
View File

@@ -317,8 +317,9 @@ module streamline
!Now accumulate into flux-tube integrals
! 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
! Only accumulate other quantities if in active domain, so 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
if (bTrc%ijk(k,XDIR) > ebGr%is+1 .and. bTrc%ijk(k+1,XDIR) > ebGr%is+1) then
bD = bD + eD*dl/bMag
bP = bP + eP*dl/bMag
bPb = bPb + ePb*dl/bMag

24
src/drivers/raijuOWDx.F90
View File

@@ -73,13 +73,6 @@ program raijuOWDx
call inpXML%Set_Val(doChmpOut,'driver/doChmpOut',.false.)
call inpXML%Set_Val(doFLOut,'driver/doFLOut',.false.)
! Init RAIJU
!call raijuInit(raiApp, inpXML)
!call raijuCpl_init(vApp, raiApp, raijuCplBase)
!if (raiApp%Model%isRestart) then
! isFirstCpl = .false.
!endif
call inpXML%Set_Val(raiApp%Model%doClockConsoleOut,'driver/doClockOut',.false.)
@@ -121,14 +114,12 @@ program raijuOWDx
if (raiApp%State%IO%doRestart(raiApp%State%t)) then
call Tic("RAIJU Restart")
!call raijuResOutput(raiApp%Model,raiApp%Grid,raiApp%State)
call raiCplApp%WriteRestart(raiApp%State%IO%nRes)
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 raiCplApp%WriteFileOutput(raiApp%State%IO%nOut)
call Toc("RAIJU Output")
@@ -166,32 +157,21 @@ program raijuOWDx
call updateRM(rmReader, rmReader%time)
call Toc("REMIX update")
! Populate RAIJU's fromV object with updated model info
! Populate RAIJU's coupling variables with updated model info
call Tic("fromV packing")
!call packFromV(raijuCplBase%fromV, vApp, rmReader, raiApp)
call packRaijuCoupler_OWD(raiCplApp, vApp, rmReader)
call Toc("fromV packing")
call Tic("fromV to State")
! Now put fomV info into RAIJU's State
!call raijuCpl_Volt2RAIJU(raijuCplBase, vApp, raiApp)
!call raijuCplApp%toIMAG(raijuCplBase, vApp, raiApp)
! Now put coupler info into RAIJU's State
call raiCplApp%toIMAG(vApp)
call Toc("fromV to State")
! Step RAIJU
call Tic("RAIJU Advance")
call raiCplApp%AdvanceModel(raiApp%Model%dt)
!call raijuAdvance(raiApp%Model,raiApp%Grid,raiApp%State, raiApp%Model%dt, isFirstCpl)
call Toc("RAIJU Advance")
!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

7
src/raiju/raijuDomain.F90
View File

@@ -109,11 +109,8 @@ module raijuDomain
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"
n_bndLim = Model%n_bndLim
! Right sweep
do j=sh%js,sh%je
bndR(j) = max(bndR(j), bndR(j-1)-n_bndLim)

6
src/raiju/raijuICHelpers.F90
View File

@@ -173,9 +173,9 @@ module raijuICHelpers
call iXML%Set_Val(L_peak,"prob/L" , 4.0 ) ! [Re]
call iXML%Set_Val(dL,"prob/dL", 0.625 ) ! [Re]
!$OMP PARALLEL DO default(shared) &
!$OMP schedule(dynamic) &
!$OMP private(i,j,L,bVol,vm,D)
! ! $OMP PARALLEL DO default(shared) &
! ! $OMP schedule(dynamic) &
! ! $OMP private(i,j,L,bVol,vm,D)
do i=Grid%shGrid%isg,Grid%shGrid%ieg
! Fully symmetric, so only need to set certain things per i
L = DipColat2L(Grid%shGrid%thc(i))

12
src/raiju/raijuIO.F90
View File

@@ -244,7 +244,6 @@ module raijuIO
enddo
call AddOutVar(IOVars,"intensity",outTmp3D(is:ie,js:je, :),uStr="1/(s*sr*keV*cm^2)")
deallocate(outTmp3D)
! Coupling things
call AddOutVar(IOVars,"dtCpl" ,State%dt) ! Attr
@@ -278,17 +277,16 @@ module raijuIO
! Moments
call AddOutVar(IOVars,"Pressure",State%Press(is:ie,js:je,:),uStr="nPa")
do s=0,Grid%nFluidIn
write(*,*)"Davg_in ",s,maxval(State%Davg(is:ie,js:je,s))
enddo
!do s=0,Grid%nFluidIn
! write(*,*)"Davg_in ",s,maxval(State%Davg(is:ie,js:je,s))
!enddo
! Add density moment as #/cc instead of amu/cc
!allocate(outDen(is:ie,js:je,Grid%nSpc+1))
allocate(outTmp3D(is:ie,js:je,Grid%nSpc+1))
outTmp3D = 0.0
do s=1, Grid%nSpc
! Convert amu/cc to #/cc
outTmp3D(:,:,s+1) = State%Den(is:ie,js:je,s+1)/Grid%spc(s)%amu
write(*,*)"Den out ",s,maxval(outTmp3D(:,:,s+1))
!write(*,*)"Den out ",s,maxval(outTmp3D(:,:,s+1))
! Don't include electrons to total number density
if(Grid%spc(s)%spcType .ne. RAIJUELE) then
outTmp3D(:,:,1) = outTmp3D(:,:,1) + outTmp3D(:,:,s+1)
@@ -298,7 +296,6 @@ module raijuIO
deallocate(outTmp3D)
! Calculate flux tube entropy using bulk pressure
!allocate(outEnt(is:ie,js:je))
allocate(outTmp2D(is:ie,js:je))
outTmp2D = -1.0
!$OMP PARALLEL DO default(shared) &
@@ -422,7 +419,6 @@ module raijuIO
allocate(outTmp2D(Grid%shGrid%isg:Grid%shGrid%ieg ,Grid%shGrid%jsg:Grid%shGrid%jeg))
call calcMapJacNorm(Grid, State%xyzMin, outTmp2D)
call AddOutVar(IOVars, "mapJacNorm", outTmp2D(is:ie,js:je), dStr="L_(2,1) norm of lat/lon => xyzMin Jacobian")
endif
call WriteVars(IOVars,.true.,Model%raijuH5, gStr)

24
src/raiju/raijuOut.F90
View File

@@ -20,7 +20,6 @@ module raijuOut
call timeStrFmt(State%t, tStr)
write (*, '(a,a,a,a,a)') ANSIGREEN, '<Writing RAIJU HDF5 DATA @ t = ', trim(tStr), ' >', ANSIRESET
endif
call WriteRAIJU(Model, Grid, State, gStr, Model%writeGhosts)
!Setup for next output
@@ -91,6 +90,9 @@ module raijuOut
character(len=strLen) :: utStr, tStr, tStr2
integer :: minDtLoc, maxDtLoc
real(rp) :: minDt, maxDt
integer :: s, sIdx
integer, dimension(2) ::mpLoc
real(rp) :: maxPress, maxDen, maxP_Xmin, maxP_Ymin
call mjd2utstr(State%mjd,utStr)
minDtLoc = minloc(State%dtk,dim=1)
@@ -107,8 +109,24 @@ module raijuOut
call timeStrFmt(State%dtk(maxDtLoc), tStr )
call timeStrFmt(State%dtk(minDtLoc), tStr2)
write(*,'(a)' ) ' Max/Min dt @ k:'
write(*,'(a,a,a,I0)') ' Max', trim(tStr2), ' @ ', maxDtLoc
write(*,'(a,a,a,I0)') ' Min', trim(tStr) , ' @ ', minDtLoc
write(*,'(a,a,a,I0)') ' Max', trim(tStr), ' @ ', maxDtLoc
write(*,'(a,a,a,I0)') ' Min', trim(tStr2) , ' @ ', minDtLoc
write(*,'(a)' ) ' Flav : max Press/Den @ Xmin,Ymin:'
do s=1, Model%nSpc
sIdx = spcIdx(Grid, Grid%spc(s)%flav)
if (Grid%spc(s)%flav == F_PSPH) then
mPLoc = maxloc(State%Den(:,:,1+sIdx))
else
mPLoc = maxloc(State%Press(:,:,1+sIdx))
endif
maxPress = State%Press(mPLoc(1), mpLoc(2), 1+sIdx)
maxDen = State%Den (mPLoc(1), mpLoc(2), 1+sIdx)
maxP_Xmin = State%xyzMincc(mPLoc(1), mpLoc(2),XDIR)
maxP_Ymin = State%xyzMincc(mPLoc(1), mpLoc(2),YDIR)
write(*,'(a,I0,a,f5.2,a,f5.2,a,f5.2,a,f5.2,a)') ' ', Grid%spc(s)%flav, ': P=', maxPress,',D=',maxDen,' @',maxP_Xmin,',',maxP_Ymin,' [Rp]'
enddo
write(*,'(a)',advance="no") ANSIRESET
State%IO%tCon = State%IO%tCon + State%IO%dtCon

18
src/raiju/raijuPreAdvancer.F90
View File

@@ -120,7 +120,7 @@ module raijuPreAdvancer
endif
!!TEST
write(*,*)"Setting last state to current indiscriminately"
!write(*,*)"Setting last state to current indiscriminately"
State%active_last = State%active
State%eta_last = State%eta
return
@@ -129,14 +129,14 @@ module raijuPreAdvancer
! we need to set its eta_last to something usable in the first eta_half calculation
! We do this in the buffer region as well, because this was freshly populated with MHD moments and any existing eta_last is stale
! After the first dt, eta_last will be set with valid information and no further adjustments are necessay
do j=sh%jsg,sh%jeg
do i=sh%isg,sh%ieg
if ( (State%active(i,j) .eq. RAIJUACTIVE .and. State%active_last(i,j) .ne. RAIJUACTIVE) &
.or. (State%active(i,j) .eq. RAIJUBUFFER) ) then
State%eta_last(i,j,:) = State%eta(i,j,:)
endif
enddo
enddo
!do j=sh%jsg,sh%jeg
! do i=sh%isg,sh%ieg
! if ( (State%active(i,j) .eq. RAIJUACTIVE .and. State%active_last(i,j) .ne. RAIJUACTIVE) &
! .or. (State%active(i,j) .eq. RAIJUBUFFER) ) then
! State%eta_last(i,j,:) = State%eta(i,j,:)
! endif
! enddo
!enddo
! NOTE: Any case not addressed here should mean that eta_last and active_last are valid in those cases

5
src/raiju/raijuStarter.F90
View File

@@ -106,14 +106,14 @@ module raijustarter
if (Model%isSA) then
tmpStr = "restart/doRes"
else
tmpStr = "Kaiju/Gamera/restart/doRes"
tmpStr = "/Kaiju/Gamera/restart/doRes"
endif
call iXML%Set_Val(Model%isRestart, trim(tmpStr),.false.)
if (Model%isRestart) then
if (Model%isSA) then
tmpStr = "restart/nRes"
else
tmpStr = "Kaiju/Gamera/restart/nRes"
tmpStr = "/Kaiju/Gamera/restart/nRes"
endif
call iXML%Set_Val(Model%nResIn, trim(tmpStr), Model%nResIn)
call genResInFname(Model, Model%ResF) ! Determine filename to read from
@@ -141,6 +141,7 @@ module raijustarter
call iXML%Set_Val(Model%nSpc, "prob/nSpc",Model%nSpc)
! Domain constraints
call iXML%Set_Val(Model%n_bndLim , "domain/bndLim" , 3)
call iXML%Set_Val(Model%maxTail_buffer, "domain/tail_buffer", def_maxTail_buffer)
call iXML%Set_Val(Model%maxSun_buffer , "domain/sun_buffer" , def_maxSun_buffer)
call iXML%Set_Val(Model%maxTail_active, "domain/tail_active", def_maxTail_active)

4
src/remix/mixmain.F90
View File

@@ -48,10 +48,10 @@ module mixmain
if (present(mixIOobj)) then
call init_grid(I(h),mixIOobj)
call init_state(I(h)%G,I(h)%St,mixIOobj%vars(:,:,:,h)) ! passing only the hemisphere variables
call init_state(I(h)%G,I(h)%St,mixIOobj%vars(:,:,:,h), shGrO=I(h)%shGr) ! passing only the hemisphere variables
else
call init_grid(I(h))
call init_state(I(h)%G,I(h)%St)
call init_state(I(h)%G,I(h)%St, shGrO=I(h)%shGr)
end if
call conductance_init(I(h)%conductance,I(h)%P,I(h)%G)

11
src/voltron/innermagsphere.F90
View File

@@ -26,6 +26,8 @@ module innermagsphere
type(XML_Input_T), intent(inout) :: iXML
character(len=strLen) :: imStr
real(rp) :: mhd_Rin = 2.0
if (.not. vApp%doDeep) return !Why are you even here?
call iXML%Set_Val(imStr,"coupling/imType","SST")
@@ -60,7 +62,14 @@ module innermagsphere
allocate(raijuCoupler_T :: vApp%imagApp)
allocate(imagOptions_T :: vApp%imagApp%opt)
vApp%imagApp%opt%swF = vApp%symh%wID
vapp%imagApp%opt%lowLatBC = vApp%mhd2chmp%lowlatBC
!vapp%imagApp%opt%lowLatBC = vApp%mhd2chmp%lowlatBC
vApp%imagApp%opt%mjd0 = vApp%MJD
! Calc lowlat BC from Gamera
associate(Gr=>gApp%Grid)
vApp%imagApp%opt%mhd_Rin = norm2(Gr%xyz(Gr%is,Gr%js,Gr%ks,:))
end associate
case DEFAULT
write(*,*) 'Unkown imType, bailing ...'
stop

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

@@ -15,6 +15,10 @@ submodule (volttypes) raijuCplTypesSub
allocate(raijuApp_T :: App%raiApp)
! Init raiju app itself
call App%raiApp%InitModel(xml)
! Update MJD with whatever voltron handed us
! If we are restarting, this will get replaced with whatever's in file later
App%raiApp%State%mjd = App%opt%mjd0
write(*,*)"MJD0=",App%opt%mjd0
! Then allocate and initialize coupling variables based on raiju app
call raijuCpl_init(App)
@@ -33,7 +37,7 @@ submodule (volttypes) raijuCplTypesSub
endif
call imagTubes2RAIJU(raiApp%Model, raiApp%Grid, raiApp%State, App%ijTubes)
! Potential
! Our App%pot should have been updated by someone. We can plug it straight into our State's copy because they live on the same grid
raiApp%State%espot(:,:) = App%pot%data(:,:)
end associate
end subroutine volt2RAIJU
@@ -110,6 +114,9 @@ submodule (volttypes) raijuCplTypesSub
real(rp), intent(in) :: dt
call App%raiApp%AdvanceModel(dt)
App%raiApp%State%t = App%raiApp%State%t + dt
App%raiApp%State%ts = App%raiApp%State%ts + 1
App%raiApp%State%mjd = T2MJD(dt,App%raiApp%State%mjd)
end subroutine

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

@@ -190,7 +190,7 @@ module raijuCplHelper
type(voltApp_T), intent(in ) :: vApp
integer :: i,j
real(rp) :: seedR
real(rp) :: seedR, eqR
type(magLine_T) :: magLine
! Get field line info and potential from voltron
! And put the data into RAIJU's fromV coupling object
@@ -199,16 +199,17 @@ module raijuCplHelper
planet=>raiCpl%raiApp%Model%planet, &
ebApp =>vApp%ebTrcApp)
seedR = planet%ri_m/planet%rp_m + TINY
seedR = planet%ri_m/planet%rp_m
! Do field line tracing, populate fromV%ijTubes
!$OMP PARALLEL DO default(shared) collapse(2) &
!$OMP PARALLEL DO default(shared) &
!$OMP schedule(dynamic) &
!$OMP private(i,j)
!$OMP private(i,j,eqR)
do i=sh%isg,sh%ieg+1
do j=sh%jsg,sh%jeg+1
call CleanLine(raiCpl%magLines(i,j))
if ((PI/2. - sh%th(i)) < raiCpl%opt%lowLatBC) then
eqR = DipColat2L(raiCpl%raiApp%Grid%thRp(i)) ! Function assumes colat coming from 1 Rp, make sure we use the right theta value
if (eqR < 0.5*raiCpl%opt%mhd_Rin) then
call DipoleTube(vApp, sh%th(i), sh%ph(j), raiCpl%ijTubes(i,j), raiCpl%magLines(i,j))
else
call MHDTube(ebApp, planet, & !ebTrcApp, planet
@@ -229,8 +230,9 @@ module raijuCplHelper
! Placeholder real-time coupling stuff
!------
! Temporary imagTube generator for realtime coupling
subroutine packRaijuCoupler_RT(raiCpl, vApp)
!! Temporary imagTube generator for realtime coupling
!! Eventually, someone else should probably be packing raiCpl objects for us
class(raijuCoupler_T), intent(inout) :: raiCpl
class(voltApp_T), intent(inout) :: vApp
@@ -243,19 +245,17 @@ module raijuCplHelper
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)
write(*,*)"asdfasdfasfasflkashd as-------"
write(*,*)shape(rmHemi%St%Vars(:,:,POT))
write(*,*)shape(tmpPot)
!rmHemi%St%pot_shGr%data = rmHemi%St%Vars(:,:,POT)
tmpPot(:,1:Np-1) = rmHemi%St%Vars(:,:,POT)
tmpPot(:,Np) = tmpPot(:,1)
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

2
src/voltron/voltapp.F90
View File

@@ -227,7 +227,7 @@ module voltapp
call InitInnerMag(vApp,gApp,xmlInp)
if(gApp%Model%isRestart) then
call vApp%imagApp%ReadRestart(gApp%Model%RunID, vApp%IO%nRes)
call vApp%imagApp%ReadRestart(gApp%Model%RunID, nRes)
!select type(rcmApp=>vApp%imagApp)
! type is (rcmIMAG_T)
! !Check if Voltron and RCM have the same restart number