mirror of
https://github.com/JHUAPL/kaiju.git
synced 2026-01-09 18:57:53 -05:00
Starting work on new loss framework
This commit is contained in:
Anthony M. Sciola
@@ -54,6 +54,35 @@ module raijutypes
|
||||
|
||||
end type raijuSpecies_T
|
||||
|
||||
|
||||
!------
|
||||
! Losses
|
||||
!------
|
||||
|
||||
type :: baseRaijuLoss_T
|
||||
!! Base loss process type for specific loss implementations to extend
|
||||
logical :: isPrecip = .false.
|
||||
|
||||
contains
|
||||
|
||||
! Default do-nothing placeholders
|
||||
procedure baseLossInit, baseLossUpdate, baseLossCalcTau, baseLossValidSpc
|
||||
|
||||
procedure :: doInit => baseLossInit
|
||||
!! Initialize any arrays and params that persist throughout run
|
||||
procedure :: doUpdate => baseLossUpdate
|
||||
!! Update any variables based on current raijuState
|
||||
procedure :: isValidSpc => baseLossValidSpc
|
||||
!! Tell someone if this loss process applies to given species
|
||||
procedure :: calcTau => baseLossCalcTau
|
||||
!! Report instantaneous loss rate for a given energy and lat/lon
|
||||
end type baseRaijuLoss_T
|
||||
|
||||
type :: raijuLPHolder_T
|
||||
!! Container for collection of loss processes
|
||||
class(baseRaijuLoss_T), allocatable :: p
|
||||
end type raijuLPHolder_T
|
||||
|
||||
!------
|
||||
! Precipitation models
|
||||
!------
|
||||
@@ -363,6 +392,10 @@ module raijutypes
|
||||
integer , dimension(:,:), allocatable :: OCBDist
|
||||
!! (Ngi, Ngj) Cell distance from open-closed boundary
|
||||
|
||||
! Loss-related things
|
||||
type(raijuLPHolder_T), dimension(:), allocatable :: lps
|
||||
!class(baseRaijuLoss_T), dimension(:), allocatable :: lps
|
||||
!! Collection of loss processes
|
||||
! (Ngi, Ngj, Nk) Varibles coming from RAIJU
|
||||
real(rp), dimension(:,:,:), allocatable :: lossRates
|
||||
!! (Ngi, Ngj, Nk) [1/s] Loss rates for each grid and lambda point. Generally stays the same over coupling time so we store them all here
|
||||
@@ -441,9 +474,44 @@ module raijutypes
|
||||
integer, intent(in) :: i,j,k
|
||||
real(rp), dimension(2) :: lossRate2
|
||||
end function raijuELossRate_T
|
||||
|
||||
end interface
|
||||
|
||||
contains
|
||||
|
||||
! New loss stuff
|
||||
subroutine baseLossInit(this, Model, Grid, xmlInp)
|
||||
!Import :: baseRaijuLoss_T, raijuModel_T, raijuGrid_T, XML_Input_T
|
||||
class(baseRaijuLoss_T), intent(inout) :: this
|
||||
type(raijuModel_T), intent(in) :: Model
|
||||
type(raijuGrid_T) , intent(in) :: Grid
|
||||
type(XML_Input_T) , intent(in) :: xmlInp
|
||||
end subroutine baseLossInit
|
||||
|
||||
subroutine baseLossUpdate(this, Model, Grid, State)
|
||||
!Import :: baseRaijuLoss_T, raijuModel_T, raijuGrid_T, raijuState_T
|
||||
class(baseRaijuLoss_T), intent(inout) :: this
|
||||
type(raijuModel_T), intent(in) :: Model
|
||||
type(raijuGrid_T) , intent(in) :: Grid
|
||||
type(raijuState_T), intent(in) :: State
|
||||
end subroutine baseLossUpdate
|
||||
|
||||
function baseLossValidSpc(this, spc) result(isValid)
|
||||
!Import :: baseRaijuLoss_T, raijuModel_T, raijuGrid_T, raijuState_T
|
||||
class(baseRaijuLoss_T), intent(in) :: this
|
||||
type(raijuSpecies_T), intent(in) :: spc
|
||||
logical :: isValid
|
||||
end function baseLossValidSpc
|
||||
|
||||
function baseLossCalcTau(this, Model, Grid, State, i, j, k) result(tau)
|
||||
!Import :: baseRaijuLoss_T, raijuModel_T, raijuGrid_T, raijuState_T, rp
|
||||
class(baseRaijuLoss_T), intent(in) :: this
|
||||
type(raijuModel_T), intent(in) :: Model
|
||||
type(raijuGrid_T) , intent(in) :: Grid
|
||||
type(raijuState_T), intent(in) :: State
|
||||
integer, intent(in) :: i, j, k
|
||||
real(rp) :: tau
|
||||
end function baseLossCalcTau
|
||||
|
||||
|
||||
|
||||
|
||||
@@ -69,6 +69,7 @@ program raijuSAx
|
||||
|
||||
! Init RAIJU
|
||||
call raijuInit(raiApp, inpXML)
|
||||
isFirstCpl = .false.
|
||||
if (raiApp%Model%isRestart) then
|
||||
isFirstCpl = .false.
|
||||
endif
|
||||
|
||||
@@ -1,4 +1,4 @@
|
||||
file(GLOB raiju_srcs_f90 *.F90)
|
||||
file(GLOB raiju_srcs_f90 *.F90 losses/*.F90)
|
||||
file(GLOB raiju_srcs_f *.F)
|
||||
if(CMAKE_Fortran_COMPILER_ID MATCHES Intel)
|
||||
set_source_files_properties(${raiju_srcs_f} PROPERTIES COMPILE_FLAGS "-fixed -nowarn")
|
||||
|
||||
168
src/raiju/losses/raijuLoss_CX.F90
Normal file
168
src/raiju/losses/raijuLoss_CX.F90
Normal file
@@ -0,0 +1,168 @@
|
||||
module raijuLoss_CX
|
||||
|
||||
use raijudefs
|
||||
use raijutypes
|
||||
use raijuSpeciesHelper
|
||||
|
||||
implicit none
|
||||
|
||||
type, extends(baseRaijuLoss_T) :: raiLoss_CX_T
|
||||
|
||||
contains
|
||||
|
||||
!procedure :: doInit => CXLossInit
|
||||
!procedure :: doPAUpdate
|
||||
procedure :: isValidSpc => CXLossValidSpc
|
||||
procedure :: calcTau => CXLossCalcTau
|
||||
|
||||
end type raiLoss_CX_T
|
||||
|
||||
contains
|
||||
|
||||
! subroutine CXLossInit(lp, Model, Grid, xmlInp)
|
||||
! class(raiLoss_CX_T), intent(inout) :: lp
|
||||
! type(raijuModel_T), intent(in) :: Model
|
||||
! type(raijuGrid_T) , intent(in) :: Grid
|
||||
! type(XML_Input_T) , intent(in) :: xmlInp
|
||||
!
|
||||
! ! Say which species we know how to do losses for
|
||||
! !call setCXValidSpecies(lp, Model, Grid)
|
||||
! end subroutine
|
||||
|
||||
|
||||
function CXLossValidSpc(this, spc) result(isValid)
|
||||
class(raiLoss_CX_T), intent(in) :: this
|
||||
type(raijuSpecies_T), intent(in) :: spc
|
||||
logical :: isValid
|
||||
|
||||
isValid = .false.
|
||||
|
||||
if ( (spc%spcType .eq. RAIJUHPLUS) &
|
||||
.or. (spc%spcType .eq. RAIJUOPLUS) ) then
|
||||
isValid = .true.
|
||||
endif
|
||||
end function
|
||||
|
||||
|
||||
function CXLossCalcTau(this, Model, Grid, State, i, j, k) result(tau)
|
||||
class(raiLoss_CX_T), intent(in) :: this
|
||||
type(raijuModel_T), intent(in) :: Model
|
||||
type(raijuGrid_T) , intent(in) :: Grid
|
||||
type(raijuState_T), intent(in) :: State
|
||||
integer, intent(in) :: i, j, k
|
||||
real(rp) :: tau
|
||||
|
||||
real(rp) :: energy, rLoc, Ngeo, cxSig, M, V
|
||||
real(rp) :: tauEq, tauBA
|
||||
|
||||
associate(spc => Grid%spc(Grid%k2spc(k)))
|
||||
|
||||
! Neutral
|
||||
rLoc = norm2(State%xyzMincc(i,j,:)) ! [Rp]
|
||||
Ngeo = OstgaardGeocorona(rLoc) ! [#/cc]
|
||||
|
||||
! Ion
|
||||
energy = abs(Grid%alamc(k))*State%bvol_cc(i,j)**(-2./3.) * 1D-3 ! keV
|
||||
cxSig = CXSigma(energy, spc%spcType)
|
||||
M = spc%amu * dalton ! [kg]
|
||||
V = sqrt(2*(energy*kev2J)/M)*100.0 ! [m/s -> cm/s]
|
||||
|
||||
! Timescale
|
||||
tauEq = 1.0/(Ngeo*V*cxSig)
|
||||
! Bounce-averaged approximation from Smith & Bewtra 1976
|
||||
! NOTE: Technically, not derived from Ostgaard distribution. Fix me soon
|
||||
! Also, angle is mirror lat, not pitch angle. So 45 deg is probably a bad estimate
|
||||
tauBA = tauEq * cos(45*PI/180.0)**3.5
|
||||
|
||||
tau = tauBA
|
||||
|
||||
end associate
|
||||
|
||||
end function CXLossCalcTau
|
||||
|
||||
|
||||
!------
|
||||
! Helpers
|
||||
!------
|
||||
! subroutine setCXValidSpecies(this, Model, Grid)
|
||||
! class(raiLoss_CX_T), intent(inout) :: this
|
||||
! type(raijuModel_T), intent(in) :: Model
|
||||
! type(raijuGrid_T) , intent(in) :: Grid
|
||||
!
|
||||
! integer :: s
|
||||
!
|
||||
! allocate(lp%isValidSpc(Grid%nSpc))
|
||||
! this%isValidSpc = .false.
|
||||
!
|
||||
! if (trim(toUpper(Model%planet%name)) .ne. "EARTH") then
|
||||
! return
|
||||
! endif
|
||||
!
|
||||
! do s=1,Grid%nSpc
|
||||
! if ( (SpcType(Grid%spc(s)) .eq. RAIJUHPLUS) &
|
||||
! .or. (SpcType(Grid%spc(s)) .eq. RAIJUOPLUS) ) then
|
||||
! this%isValidSpc(s) = .true.
|
||||
! endif
|
||||
! enddo
|
||||
!
|
||||
! end subroutine setCXValidSpecies
|
||||
|
||||
|
||||
function OstgaardGeocorona(L) result(Ngeo)
|
||||
!! Geocoronal density afa L [#/cc], Taken from Ostgaard 2003
|
||||
real(rp), intent(in) :: L
|
||||
!! L shell (radial dist?) [Rp]
|
||||
real(rp) :: Ngeo
|
||||
!! Density [#/cc]
|
||||
|
||||
Ngeo = 10000.0*exp(-L/1.02) + 70.0*exp(-L/8.2)
|
||||
end function OstgaardGeocorona
|
||||
|
||||
|
||||
!Charge exchange cross-section for K [keV] and species ispc
|
||||
!Sig in cm2
|
||||
!Using Lindsay & Stebbings 2005
|
||||
function CXSigma(K,ispc) result(Sig)
|
||||
real(rp), intent(in) :: K
|
||||
!! Energy [keV]
|
||||
integer, intent(in) :: ispc
|
||||
!! RAIJU enum for species type
|
||||
real(rp) :: Sig
|
||||
!! Cross-section [cm^2]
|
||||
real(rp) :: Sig0, KSig,a1,a2,a3,B1,B2
|
||||
|
||||
Sig0 = 1.0e-16
|
||||
|
||||
select case(ispc)
|
||||
case(RAIJUHPLUS)
|
||||
!Charge exchange cross-section for H+/H
|
||||
!Cap for validity of CX cross-section
|
||||
KSig = K
|
||||
call ClampValue(KSig,0.005_rp,250.0_rp)
|
||||
|
||||
a1 = 4.15
|
||||
a2 = 0.531
|
||||
a3 = 67.3
|
||||
|
||||
B1 = (a1-a2*log(KSig))**2.0
|
||||
B2 = 1.0-exp(-a3/KSig)
|
||||
Sig = Sig0*B1*(B2**(4.5))
|
||||
case(RAIJUOPLUS)
|
||||
!Charge exchange cross-section for O+/H
|
||||
!Cap for validity of CX cross-section
|
||||
KSig = K
|
||||
call ClampValue(KSig,0.025_rp,600.0_rp)
|
||||
a1 = 3.13
|
||||
a2 = 0.170
|
||||
a3 = 87.5
|
||||
|
||||
B1 = (a1-a2*log(KSig))**2.0
|
||||
B2 = 1.0-exp(-a3/KSig)
|
||||
Sig = Sig0*B1*(B2**(0.8))
|
||||
case default
|
||||
Sig = 0.0
|
||||
end select
|
||||
|
||||
END FUNCTION CXSigma
|
||||
|
||||
end module raijuLoss_CX
|
||||
@@ -83,16 +83,22 @@ module raijuPreAdvancer
|
||||
call EvalMoments(Grid, State)
|
||||
call Toc("Moments Eval PreAdvance")
|
||||
|
||||
call Tic("Calc loss rates")
|
||||
call Tic("Losses")
|
||||
if (Model%doLosses) then
|
||||
call Tic("Update loss states")
|
||||
call updateRaiLosses(Model, Grid, State)
|
||||
call Toc("Update loss states")
|
||||
|
||||
call Tic("Calc loss rates")
|
||||
!$OMP PARALLEL DO default(shared) &
|
||||
!$OMP schedule(dynamic) &
|
||||
!$OMP private(k)
|
||||
do k=1,Grid%Nk
|
||||
call calcChannelLossRates(Model, Grid, State, k)
|
||||
enddo
|
||||
call Toc("Calc loss rates")
|
||||
endif
|
||||
call Toc("Calc loss rates")
|
||||
call Toc("Losses")
|
||||
|
||||
end subroutine raijuPreAdvance
|
||||
|
||||
|
||||
@@ -15,6 +15,7 @@ module raijustarter
|
||||
use raijuout
|
||||
use raijuICHelpers
|
||||
use raijuELossWM
|
||||
use raijulosses, only : initRaiLosses
|
||||
|
||||
! Cmake points to this
|
||||
use raijuuseric
|
||||
@@ -196,19 +197,20 @@ module raijustarter
|
||||
call iXML%Set_Val(Model%doCC , "losses/doCC" ,.true. )
|
||||
call iXML%Set_Val(Model%doCX , "losses/doCX" ,.false.)
|
||||
call iXML%Set_Val(Model%doFLC , "losses/doFLC",.false.)
|
||||
!call iXML%Set_Val(Model%doEWM , "losses/doELoss",.false.)
|
||||
|
||||
! Electron loss model
|
||||
call iXML%Set_Val(tmpStr, "losses/eLossModel","WM")
|
||||
select case (tmpStr)
|
||||
case ("WM")
|
||||
write(*,*) "(RAIJU) Using Wang-Bao electron wave model"
|
||||
Model%eLossModel = RaiELOSS_WM
|
||||
Model%eLossRateFn => calcELossRate_WM
|
||||
! We init later now, up in raijuInit, since initEWM needs Grid info for diagnostic output
|
||||
case default
|
||||
write(*,*) "(RAIJU) Did not get a valid electron loss model, goodbye"
|
||||
stop
|
||||
end select
|
||||
!call iXML%Set_Val(tmpStr, "losses/eLossModel","WM")
|
||||
!select case (tmpStr)
|
||||
! case ("WM")
|
||||
! write(*,*) "(RAIJU) Using Wang-Bao electron wave model"
|
||||
! Model%eLossModel = RaiELOSS_WM
|
||||
! Model%eLossRateFn => calcELossRate_WM
|
||||
! ! We init later now, up in raijuInit, since initEWM needs Grid info for diagnostic output
|
||||
! case default
|
||||
! write(*,*) "(RAIJU) Did not get a valid electron loss model, goodbye"
|
||||
! stop
|
||||
!end select
|
||||
|
||||
call iXML%Set_Val(Model%doFatOutput, "output/doFat",.false.)
|
||||
! TODO: Add flags to output certain data, like coupling information
|
||||
@@ -438,6 +440,9 @@ module raijustarter
|
||||
|
||||
call Model%initState(Model, Grid, State, iXML)
|
||||
|
||||
! Do losses after everything else has been set, just in case they need something from it
|
||||
call initRaiLosses(Model, Grid, State, iXML)
|
||||
|
||||
end subroutine raijuInitState
|
||||
|
||||
|
||||
|
||||
@@ -1,36 +1,166 @@
|
||||
module raijulosses
|
||||
|
||||
use kdefs
|
||||
use XML_Input
|
||||
|
||||
use raijudefs
|
||||
use raijutypes
|
||||
use raijuspecieshelper, only : spcIdx
|
||||
|
||||
use raijuLoss_CX
|
||||
|
||||
implicit none
|
||||
|
||||
contains
|
||||
|
||||
|
||||
subroutine initRaiLosses(Model, Grid, State, xmlInp)
|
||||
type(raijuModel_T), intent(in) :: Model
|
||||
type(raijuGrid_T), intent(in) :: Grid
|
||||
type(raijuState_T), intent(inout) :: State
|
||||
type(XML_Input_T), intent(in) :: xmlInp
|
||||
|
||||
integer :: iLP
|
||||
!! iterator
|
||||
integer :: numLPs
|
||||
!! number of loss proccces we're gonna have
|
||||
integer :: initCX, initCC, initSS, initFLC = 0
|
||||
!! Flag for if we need to allocate and init this LP
|
||||
|
||||
State%lossRates = 0.0
|
||||
State%lossRatesPrecip = 0.0
|
||||
|
||||
numLPs = 0
|
||||
|
||||
if (Model%doCX ) initCX = 1
|
||||
if (Model%doCC ) initCC = 1
|
||||
if (Model%doSS ) initSS = 1
|
||||
if (Model%doFLC) initFLC = 1
|
||||
numLPs = initCX + initCC + initSS + initFLC
|
||||
|
||||
allocate(State%lps(numLPs))
|
||||
|
||||
do iLP=1,numLPs
|
||||
! Determine which loss process is next in line for initting
|
||||
if (initCX==1) then
|
||||
if (allocated(State%lps(iLP)%p)) deallocate(State%lps(iLP)%p)
|
||||
allocate( raiLoss_CX_T :: State%lps(iLP)%p )
|
||||
initCX = 0
|
||||
elseif(.false.) then
|
||||
continue
|
||||
endif
|
||||
|
||||
! Init newly-allocated LP
|
||||
call State%lps(iLP)%p%doInit(Model, Grid, xmlInp)
|
||||
enddo
|
||||
|
||||
end subroutine initRaiLosses
|
||||
|
||||
|
||||
subroutine updateRaiLosses(Model, Grid, State)
|
||||
!! Update loss processes with step-specific information
|
||||
type(raijuModel_T), intent(in) :: Model
|
||||
type(raijuGrid_T) , intent(in) :: Grid
|
||||
type(raijuState_T), intent(inout) :: State
|
||||
|
||||
integer :: nLP, iLP
|
||||
|
||||
if (allocated(State%lps)) then
|
||||
nLP = size(State%lps)
|
||||
else
|
||||
nLP = 0.0
|
||||
return
|
||||
endif
|
||||
|
||||
do iLP=1,nLP
|
||||
call State%lps(iLP)%p%doUpdate(Model, Grid, State)
|
||||
enddo
|
||||
end subroutine updateRaiLosses
|
||||
|
||||
|
||||
subroutine calcChannelLossRates(Model, Grid, State, k)
|
||||
!! Calculate 2D loss rates for channel k
|
||||
!! Usually this will stay constant over a coupling period, so it can be called during pre-advance and not touched afterwards
|
||||
type(raijuModel_T), intent(in) :: Model
|
||||
type(raijuGrid_T) , intent(in) :: Grid
|
||||
type(raijuState_T), intent(inout) :: State
|
||||
integer, intent(in) :: k
|
||||
|
||||
integer :: i,j,l
|
||||
real(rp) :: rate, rateSS, ratePrecip
|
||||
logical, dimension(Grid%shGrid%isg:Grid%shGrid%ieg, &
|
||||
Grid%shGrid%jsg:Grid%shGrid%jeg) :: isG
|
||||
|
||||
|
||||
! Electrons are special, handle them on their own
|
||||
if (Grid%spc(Grid%k2spc(k))%spcType .eq. RAIJUELE) then
|
||||
!call calcElectronLossRate(Model, Grid, State, k)
|
||||
return
|
||||
endif
|
||||
|
||||
! Otherwise, do default loss behavior
|
||||
|
||||
State%lossRates(:,:,k) = 0.0 ! 1/s, so 0 means we lose nothing no matter the dt
|
||||
State%lossRatesPrecip(:,:,k) = 0.0
|
||||
|
||||
! Mask inactive, go ahead and calc losses in buffer just in case anyone wants it
|
||||
where (State%active .ne. RAIJUINACTIVE)
|
||||
isG = .true.
|
||||
elsewhere
|
||||
isG = .false.
|
||||
end where
|
||||
|
||||
associate(spc=>Grid%spc(Grid%k2spc(k)), lps=>State%lps)
|
||||
|
||||
do j=Grid%shGrid%jsg,Grid%shGrid%jeg
|
||||
do i=Grid%shGrid%isg,Grid%shGrid%ieg
|
||||
do l=1,size(lps)
|
||||
if (.not. isG(i,j)) then
|
||||
cycle
|
||||
endif
|
||||
if ( .not. lps(l)%p%isValidSpc(spc) ) then
|
||||
cycle
|
||||
endif
|
||||
|
||||
rate = 1.0_rp/lps(l)%p%calcTau(Model, Grid, State, i,j,k)
|
||||
|
||||
State%lossRates(i,j,k) = State%lossRates(i,j,k) + rate
|
||||
if (lps(l)%p%isPrecip) then
|
||||
State%lossRatesPrecip(i,j,k) = State%lossRatesPrecip(i,j,k) + rate
|
||||
endif
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
|
||||
end associate
|
||||
|
||||
|
||||
end subroutine calcChannelLossRates
|
||||
|
||||
! ----- OLD -----
|
||||
|
||||
|
||||
!------
|
||||
! High-level calc
|
||||
!------
|
||||
|
||||
subroutine calcChannelLossRates(Model, Grid, State, k)
|
||||
!! Calculate 2D loss rates for channel k
|
||||
!! Usually this will stay constant over a coupling period, so it can be called during pre-advance and not touched afterwards
|
||||
type(raijuModel_T), intent(inout) :: Model
|
||||
type(raijuGrid_T), intent(in) :: Grid
|
||||
type(raijuState_T), intent(inout) :: State
|
||||
integer, intent(in) :: k
|
||||
|
||||
if (Grid%spc(Grid%k2spc(k))%spcType .eq. RAIJUHPLUS) then
|
||||
call calcProtonLossRate(Model, Grid, State, k)
|
||||
endif
|
||||
|
||||
if (Grid%spc(Grid%k2spc(k))%spcType .eq. RAIJUELE) then
|
||||
call calcElectronLossRate(Model, Grid, State, k)
|
||||
endif
|
||||
|
||||
end subroutine calcChannelLossRates
|
||||
!subroutine calcChannelLossRates(Model, Grid, State, k)
|
||||
! !! Calculate 2D loss rates for channel k
|
||||
! !! Usually this will stay constant over a coupling period, so it can be called during pre-advance and not touched afterwards
|
||||
! type(raijuModel_T), intent(inout) :: Model
|
||||
! type(raijuGrid_T), intent(in) :: Grid
|
||||
! type(raijuState_T), intent(inout) :: State
|
||||
! integer, intent(in) :: k
|
||||
!
|
||||
! if (Grid%spc(Grid%k2spc(k))%spcType .eq. RAIJUHPLUS) then
|
||||
! call calcProtonLossRate(Model, Grid, State, k)
|
||||
! endif
|
||||
!
|
||||
! if (Grid%spc(Grid%k2spc(k))%spcType .eq. RAIJUELE) then
|
||||
! call calcElectronLossRate(Model, Grid, State, k)
|
||||
! endif
|
||||
!
|
||||
!end subroutine calcChannelLossRates
|
||||
|
||||
|
||||
subroutine calcProtonLossRate(Model, Grid, State, k)
|
||||
|
||||
Reference in New Issue
Block a user