From db1a5c0f66f2bd263e02320f7c00c2031a2b98b8 Mon Sep 17 00:00:00 2001 From: Anthony Date: Fri, 5 Sep 2025 12:20:00 -0600 Subject: [PATCH 01/15] Removing allocatable from fillArrays so we can pass it static arrays too --- src/base/arrayutil.F90 | 20 ++++++++++---------- 1 file changed, 10 insertions(+), 10 deletions(-) diff --git a/src/base/arrayutil.F90 b/src/base/arrayutil.F90 index 921e61b8..cda9ccb6 100644 --- a/src/base/arrayutil.F90 +++ b/src/base/arrayutil.F90 @@ -15,7 +15,7 @@ module arrayutil subroutine fillArray1D_R(array, val) - real(rp), dimension(:), allocatable, intent(inout) :: array + real(rp), dimension(:), intent(inout) :: array real(rp), intent(in) :: val integer :: lbds(1),ubds(1),i @@ -37,7 +37,7 @@ module arrayutil end subroutine fillArray1D_R subroutine fillArray1D_I(array, val) - integer, dimension(:), allocatable, intent(inout) :: array + integer, dimension(:), intent(inout) :: array integer, intent(in) :: val integer :: lbds(1),ubds(1),i @@ -59,7 +59,7 @@ module arrayutil end subroutine fillArray1D_I subroutine fillArray2D_R(array, val) - real(rp), dimension(:,:), allocatable, intent(inout) :: array + real(rp), dimension(:,:), intent(inout) :: array real(rp), intent(in) :: val integer :: lbds(2),ubds(2),i,j @@ -83,7 +83,7 @@ module arrayutil end subroutine fillArray2D_R subroutine fillArray2D_I(array, val) - integer, dimension(:,:), allocatable, intent(inout) :: array + integer, dimension(:,:), intent(inout) :: array integer, intent(in) :: val integer :: lbds(2),ubds(2),i,j @@ -107,7 +107,7 @@ module arrayutil end subroutine fillArray2D_I subroutine fillArray3D_R(array, val) - real(rp), dimension(:,:,:), allocatable, intent(inout) :: array + real(rp), dimension(:,:,:), intent(inout) :: array real(rp), intent(in) :: val integer :: lbds(3),ubds(3),i,j,k @@ -133,7 +133,7 @@ module arrayutil end subroutine fillArray3D_R subroutine fillArray3D_I(array, val) - integer, dimension(:,:,:), allocatable, intent(inout) :: array + integer, dimension(:,:,:), intent(inout) :: array integer, intent(in) :: val integer :: lbds(3),ubds(3),i,j,k @@ -159,7 +159,7 @@ module arrayutil end subroutine fillArray3D_I subroutine fillArray4D_R(array, val) - real(rp), dimension(:,:,:,:), allocatable, intent(inout) :: array + real(rp), dimension(:,:,:,:), intent(inout) :: array real(rp), intent(in) :: val integer :: lbds(4),ubds(4),i,j,k @@ -185,7 +185,7 @@ module arrayutil end subroutine fillArray4D_R subroutine fillArray4D_I(array, val) - integer, dimension(:,:,:,:), allocatable, intent(inout) :: array + integer, dimension(:,:,:,:), intent(inout) :: array integer, intent(in) :: val integer :: lbds(4),ubds(4),i,j,k @@ -211,7 +211,7 @@ module arrayutil end subroutine fillArray4D_I subroutine fillArray5D_R(array, val) - real(rp), dimension(:,:,:,:,:), allocatable, intent(inout) :: array + real(rp), dimension(:,:,:,:,:), intent(inout) :: array real(rp), intent(in) :: val integer :: lbds(5),ubds(5),i,j,k @@ -237,7 +237,7 @@ module arrayutil end subroutine fillArray5D_R subroutine fillArray5D_I(array, val) - integer, dimension(:,:,:,:,:), allocatable, intent(inout) :: array + integer, dimension(:,:,:,:,:), intent(inout) :: array integer, intent(in) :: val integer :: lbds(5),ubds(5),i,j,k From 207bb848107c8d0d70bca4121923e11f69dd0a82 Mon Sep 17 00:00:00 2001 From: Anthony Date: Fri, 5 Sep 2025 12:27:42 -0600 Subject: [PATCH 02/15] Re-working raiju coldStart to happen inside raiju Proper, so we can apply moments-->eta and then optionally add coldstart plasma to it. Compiles, not tested --- src/base/types/raijuTypes.F90 | 5 + src/raiju/raijuBCs.F90 | 9 +- src/raiju/raijuColdStartHelper.F90 | 118 +++++++++--------- src/raiju/raijuEtautils.F90 | 8 +- src/raiju/raijuIO.F90 | 11 ++ src/raiju/raijuOut.F90 | 2 +- src/raiju/raijuPreAdvancer.F90 | 3 + .../modelInterfaces/raiCplTypesSub.F90 | 49 +------- 8 files changed, 98 insertions(+), 107 deletions(-) diff --git a/src/base/types/raijuTypes.F90 b/src/base/types/raijuTypes.F90 index 3c13a82d..9a35b0fc 100644 --- a/src/base/types/raijuTypes.F90 +++ b/src/base/types/raijuTypes.F90 @@ -147,6 +147,11 @@ module raijutypes !--- State ---! + logical :: doCS_next_preAdv = .false. + !! Signal to run coldstart next time raiju preAdvances + real(rp) :: modelDst_next_preAdv = 0.0_rp + !! Target Dst [nT] when we run coldstart next + logical :: doneFirstCS = .false. !! Have we executed once already? real(rp) :: lastEval = -1*HUGE diff --git a/src/raiju/raijuBCs.F90 b/src/raiju/raijuBCs.F90 index f527e46a..f4bd3acf 100644 --- a/src/raiju/raijuBCs.F90 +++ b/src/raiju/raijuBCs.F90 @@ -5,6 +5,7 @@ module raijuBCs use raijutypes use raijuetautils use raijudomain + use raijuColdStartHelper implicit none @@ -29,12 +30,14 @@ module raijuBCs doWholeDomain = .false. endif - ! Now that topo is set, we can calculate active domain - call setActiveDomain(Model, Grid, State) - call calcMomentIngestionLocs(Model, Grid, State, doWholeDomain, doMomentIngest) call applyMomentIngestion(Model, Grid, State, doMomentIngest) + if (State%coldStarter%doCS_next_preAdv) then + call raijuGeoColdStart(Model, Grid, State, State%t, State%coldStarter%modelDst_next_preAdv, doAccumulateO=.true.) + State%coldStarter%doCS_next_preAdv = .false. + endif + if (Model%doActiveShell ) then ! Do first round of determining active shells for each k diff --git a/src/raiju/raijuColdStartHelper.F90 b/src/raiju/raijuColdStartHelper.F90 index f14c261d..43d28ccf 100644 --- a/src/raiju/raijuColdStartHelper.F90 +++ b/src/raiju/raijuColdStartHelper.F90 @@ -4,6 +4,7 @@ module raijuColdStartHelper use raijutypes use imaghelper use earthhelper + use arrayutil use raijuetautils use raijuloss_CX @@ -41,7 +42,7 @@ module raijuColdStartHelper ! Worker routines !------ - subroutine raijuGeoColdStart(Model, Grid, State, t0, dstModel) + subroutine raijuGeoColdStart(Model, Grid, State, t0, dstModel, doAccumulateO) !! Cold start RAIJU assuming we are at Earth sometime around 21st century type(raijuModel_T), intent(in) :: Model type(raijuGrid_T), intent(in) :: Grid @@ -50,24 +51,28 @@ module raijuColdStartHelper !! Target time to pull SW values from real(rp), intent(in) :: dstModel !! Current dst of global model + logical, optional, intent(in) :: doAccumulateO + !! If true, keep State%eta and add coldstart to it. + !! If false, replace State%eta with coldStart info + !! Default: false - logical :: isFirstCS + logical :: isFirstCS, doAccumulate integer :: s, sIdx_p, sIdx_e real(rp) :: dstReal, dstTarget real(rp) :: dps_current, dps_preCX, dps_postCX, dps_rescale, dps_ele real(rp) :: etaScale - logical, dimension(Grid%shGrid%isg:Grid%shGrid%ieg, Grid%shGrid%jsg:Grid%shGrid%jeg) :: isGood + logical , dimension(Grid%shGrid%isg:Grid%shGrid%ieg, Grid%shGrid%jsg:Grid%shGrid%jeg) :: isGood + real(rp), dimension(Grid%shGrid%isg:Grid%shGrid%ieg, Grid%shGrid%jsg:Grid%shGrid%jeg, Grid%Nk) :: etaCS - associate(cs=>State%coldStarter) - !write(*,*)"Coldstart running..." - - - isFirstCS = .not. State%coldStarter%doneFirstCS - - if (.not. isFirstCS .and. .not. cs%doUpdate) then - return + if (present(doAccumulateO)) then + doAccumulate = doAccumulateO + else + doAccumulate = .false. endif + call fillArray(etaCS, 0.0_rp) + + associate(cs=>State%coldStarter) where (State%active .eq. RAIJUACTIVE) isGood = .true. @@ -77,34 +82,17 @@ module raijuColdStartHelper sIdx_p = spcIdx(Grid, F_HOTP) sIdx_e = spcIdx(Grid, F_HOTE) - - - if (isFirstCS) then - ! Start by nuking all etas we will set up ourselves - do s=1,Grid%nSpc - !! Skip plasmashere, let that be handled on its own - !if ( Grid%spc(s)%flav == F_PSPH) then - ! continue - !endif - State%eta(:,:,Grid%spc(s)%kStart:Grid%spc(s)%kEnd) = 0.0 - enddo - endif - - !! Init psphere - !if (isFirstCS .or. cs%doPsphUpdate) then - ! call setRaijuInitPsphere(Model, Grid, State, Model%psphInitKp) - !endif ! Update Dst target dstReal = GetSWVal('symh', Model%tsF, t0) - if (isFirstCS) then + if (.not. cs%doneFirstCS) then ! On first try, we assume there is no existing ring current, and its our job to make up the entire difference dstTarget = dstReal - dstModel else if (t0 > (cs%lastEval + cs%evalCadence)) then ! If we are updating, there should already be some ring current ! If dstReal - dstModel is still < 0, we need to add ADDITIONAL pressure to get them to match - dps_current = spcEta2DPS(Model, Grid, State, Grid%spc(sIdx_p), isGood) + spcEta2DPS(Model, Grid, State, Grid%spc(sIdx_e), isGood) + dps_current = spcEta2DPS(Model, Grid, State%bvol_cc, State%eta, Grid%spc(sIdx_p), isGood) + spcEta2DPS(Model, Grid, State%bvol_cc, State%eta, Grid%spc(sIdx_e), isGood) dstTarget = dstReal - (dstModel - dps_current) else ! Otherwise we have nothing to do, just chill til next update time @@ -118,28 +106,26 @@ module raijuColdStartHelper return endif - if (isFirstCS) then - ! Init psphere + if (.not. cs%doneFirstCS) then + ! Init psphere (doesn't care about accumulation, always hard resets) call setRaijuInitPsphere(Model, Grid, State, Model%psphInitKp) ! Init hot protons - call raiColdStart_initHOTP(Model, Grid, State, t0, dstTarget) - dps_preCX = spcEta2DPS(Model, Grid, State, Grid%spc(sIdx_p), isGood) + call raiColdStart_initHOTP(Model, Grid, State, t0, dstTarget, etaCS) + dps_preCX = spcEta2DPS(Model, Grid, State%bvol_cc, etaCS, Grid%spc(sIdx_p), isGood) ! Hit it with some charge exchange if (cs%doCX) then - call raiColdStart_applyCX(Model, Grid, State, Grid%spc(sIdx_p)) + call raiColdStart_applyCX(Model, Grid, State, Grid%spc(sIdx_p), etaCS) endif - dps_postCX = spcEta2DPS(Model, Grid, State, Grid%spc(sIdx_p), isGood) - ! Calc moments to update pressure and density - call EvalMoments(Grid, State) + dps_postCX = spcEta2DPS(Model, Grid, State%bvol_cc, etaCS, Grid%spc(sIdx_p), isGood) ! Use HOTP moments to set electrons - call raiColdStart_initHOTE(Model, Grid, State) - dps_ele = spcEta2DPS(Model, Grid, State, Grid%spc(sIdx_e), isGood) + call raiColdStart_initHOTE(Model, Grid, State, etaCS) + dps_ele = spcEta2DPS(Model, Grid, State%bvol_cc, etaCS, Grid%spc(sIdx_e), isGood) dps_current = dps_postCX ! Note: if using fudge we're gonna lose electrons immediately, don't include them in current dst for now endif etaScale = abs(dstTarget / dps_current) - State%eta(:,:,Grid%spc(sIdx_p)%kStart:Grid%spc(sIdx_p)%kEnd) = etaScale*State%eta(:,:,Grid%spc(sIdx_p)%kStart:Grid%spc(sIdx_p)%kEnd) - dps_rescale = spcEta2DPS(Model, Grid, State, Grid%spc(sIdx_p), isGood) + etaCS(:,:,Grid%spc(sIdx_p)%kStart:Grid%spc(sIdx_p)%kEnd) = etaScale*etaCS(:,:,Grid%spc(sIdx_p)%kStart:Grid%spc(sIdx_p)%kEnd) + dps_rescale = spcEta2DPS(Model, Grid, State%bvol_cc, etaCS, Grid%spc(sIdx_p), isGood) if (isfirstCS) then write(*,*) "RAIJU Cold starting..." @@ -161,18 +147,26 @@ module raijuColdStartHelper end associate + ! finally, put it into raiju state + if(doAccumulate) then + State%eta = State%eta + etaCS + else + State%eta = etaCS + endif + State%coldStarter%doneFirstCS = .true. end subroutine raijuGeoColdStart - subroutine raiColdStart_initHOTP(Model, Grid, State, t0, dstTarget) + subroutine raiColdStart_initHOTP(Model, Grid, State, t0, dstTarget, etaCS) type(raijuModel_T), intent(in) :: Model type(raijuGrid_T), intent(in) :: Grid - type(raijuState_T), intent(inout) :: State + type(raijuState_T), intent(in) :: State real(rp), intent(in) :: t0 !! Target time to pull SW values from real(rp), intent(in) :: dstTarget + real(rp), dimension(Grid%shGrid%isg:Grid%shGrid%ieg, Grid%shGrid%jsg:Grid%shGrid%jeg, Grid%Nk), intent(inout) :: etaCS real(rp) :: dstTarget_p logical :: isInTM03 @@ -218,7 +212,7 @@ module raijuColdStartHelper endif ! Set everything to zero to start - State%eta(:,:,spc%kStart:spc%kEnd) = 0.0_rp + etaCS(:,:,spc%kStart:spc%kEnd) = 0.0_rp ! Now set our initial density and pressure profile do j=sh%jsg,sh%jeg @@ -250,7 +244,7 @@ module raijuColdStartHelper endif ! Finally map it to HOTP etas - call DkT2SpcEta(Model, spc, State%eta(i,j,spc%kStart:spc%kEnd), D_final, kt_rc, vm) + call DkT2SpcEta(Model, spc, etaCS(i,j,spc%kStart:spc%kEnd), D_final, kt_rc, vm) enddo enddo @@ -260,11 +254,12 @@ module raijuColdStartHelper end subroutine raiColdStart_initHOTP - subroutine raiColdStart_applyCX(Model, Grid, State, spc) + subroutine raiColdStart_applyCX(Model, Grid, State, spc, etaCS) type(raijuModel_T), intent(in) :: Model type(raijuGrid_T), intent(in) :: Grid - type(raijuState_T), intent(inout) :: State + type(raijuState_T), intent(in) :: State type(raijuSpecies_T), intent(in) :: spc + real(rp), dimension(Grid%shGrid%isg:Grid%shGrid%ieg, Grid%shGrid%jsg:Grid%shGrid%jeg, Grid%Nk), intent(inout) :: etaCS integer :: i,j,k type(raiLoss_CX_T) :: lossCX @@ -281,7 +276,7 @@ module raijuColdStartHelper do i=Grid%shGrid%isg,Grid%shGrid%ieg do k = spc%kStart,spc%kEnd tau = lossCX%calcTau(Model, Grid, State, i, j, k) - State%eta(i,j,k) = State%eta(i,j,k)*exp(-tCX/tau) + etaCS(i,j,k) = etaCS(i,j,k)*exp(-tCX/tau) enddo enddo enddo @@ -289,20 +284,24 @@ module raijuColdStartHelper end subroutine raiColdStart_applyCX - subroutine raiColdStart_initHOTE(Model, Grid, State) + subroutine raiColdStart_initHOTE(Model, Grid, State, etaCS) type(raijuModel_T), intent(in) :: Model type(raijuGrid_T), intent(in) :: Grid - type(raijuState_T), intent(inout) :: State + type(raijuState_T), intent(in) :: State + real(rp), dimension(Grid%shGrid%isg:Grid%shGrid%ieg, Grid%shGrid%jsg:Grid%shGrid%jeg, Grid%Nk), intent(inout) :: etaCS integer :: sIdx_e, sIdx_p integer :: i,j + real(rp) :: press_p, den_p real(rp) :: kt_p, kt_e, den, vm sIdx_p = spcIdx(Grid, F_HOTP) sIdx_e = spcIdx(Grid, F_HOTE) + + associate(spc_p=>Grid%spc(sIdx_p), spc_e=>Grid%spc(sIdx_e)) ! Set everything to zero to start - State%eta(:,:,Grid%spc(sIdx_e)%kStart:Grid%spc(sIdx_e)%kEnd) = 0.0_rp + etaCS(:,:,spc_e%kStart:spc_e%kEnd) = 0.0_rp !$OMP PARALLEL DO default(shared) & !$OMP schedule(dynamic) & @@ -312,15 +311,22 @@ module raijuColdStartHelper if (State%active(i,j) .eq. RAIJUINACTIVE) cycle vm = State%bvol_cc(i,j)**(-2./3.) - den = State%Den(sIdx_p)%data(i,j) - kt_p = DP2kT(den, State%Press(sIdx_p)%data(i,j)) + !den = State%Den(sIdx_p)%data(i,j) + !kt_p = DP2kT(den, State%Press(sIdx_p)%data(i,j)) + + den_p = SpcEta2Den (spc_p, etaCS(i,j,spc_p%kStart:spc_p%kEnd), State%bvol_cc(i,j)) + press_p = SpcEta2Press(spc_p, etaCS(i,j,spc_p%kStart:spc_p%kEnd), State%bvol_cc(i,j)) + kt_p = DP2kT(den_p, press_p) + kt_e = kt_p / Model%tiote call DkT2SpcEta(Model, Grid%spc(sIdx_e), & - State%eta(i,j,Grid%spc(sIdx_e)%kStart:Grid%spc(sIdx_e)%kEnd), & - den, kt_e, vm) + etaCS(i,j,Grid%spc(sIdx_e)%kStart:Grid%spc(sIdx_e)%kEnd), & + den_p, kt_e, vm) enddo enddo + end associate + end subroutine raiColdStart_initHOTE diff --git a/src/raiju/raijuEtautils.F90 b/src/raiju/raijuEtautils.F90 index ff368587..3c6d248e 100644 --- a/src/raiju/raijuEtautils.F90 +++ b/src/raiju/raijuEtautils.F90 @@ -178,11 +178,12 @@ module raijuetautils end function SpcEta2Press - function spcEta2DPS(Model, Grid, State, spc, isGood) result(dpsdst) + function spcEta2DPS(Model, Grid, bVol_cc, eta, spc, isGood) result(dpsdst) !! Calculate total DPS-Dst for given species within the defined isGood domain 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, Grid%shGrid%jsg:Grid%shGrid%jeg), intent(in) :: bVol_cc + real(rp), dimension(Grid%shGrid%isg:Grid%shGrid%ieg, Grid%shGrid%jsg:Grid%shGrid%jeg, Grid%Nk), intent(in) :: eta type(raijuSpecies_T), intent(in) :: spc logical, dimension(Grid%shGrid%isg:Grid%shGrid%ieg, Grid%shGrid%jsg:Grid%shGrid%jeg), intent(in) :: isGood !! Eval mask, true = point is included in calculation @@ -197,8 +198,7 @@ module raijuetautils do j=Grid%shGrid%jsg,Grid%shGrid%jeg do i=Grid%shGrid%isg,Grid%shGrid%ieg if (.not. isGood(i,j)) cycle - bVol = State%bvol_cc(i,j) - press = SpcEta2Press(spc, State%eta(i,j,spc%kStart:spc%kEnd), bVol) ! [nPa] + press = SpcEta2Press(spc, eta(i,j,spc%kStart:spc%kEnd), bvol_cc(i,j)) ! [nPa] energyDen = (press*1.0D-9) * (bVol*Model%planet%rp_m*1.0D9) * (Grid%Brcc(i,j)*1.0D-9)/kev2J ! p[J/m^3] * bVol[m/T] * B[T] = [J/m^2] * keV/J = [keV/m^2] energy = energyDen*(Grid%areaCC(i,j)*Model%planet%ri_m**2) ! [keV/m^2]* Re^2[m^2] = [keV] dpsdst = dpsdst - 4.2*(1.0D-30)*energy ! [nT] diff --git a/src/raiju/raijuIO.F90 b/src/raiju/raijuIO.F90 index 889f8ea3..3fe18b36 100644 --- a/src/raiju/raijuIO.F90 +++ b/src/raiju/raijuIO.F90 @@ -503,6 +503,12 @@ module raijuIO else call AddOutVar(IOVars,"cs_doneFirstCS", 0) endif + if (State%coldStarter%doCS_next_preAdv) then + call AddOutVar(IOVars,"cs_doCS_next_preAdv", 1) + else + call AddOutVar(IOVars,"cs_doCS_next_preAdv", 0) + endif + call AddOutVar(IOVars, "cs_modelDst_next_preAdv", State%coldStarter%modelDst_next_preAdv) call AddOutVar(IOVars, "cs_lastEval", State%coldStarter%lastEval) call AddOutVar(IOVars, "cs_lastTarget", State%coldStarter%lastTarget) @@ -598,6 +604,8 @@ module raijuIO call AddInVar(IOVars,"nRes") call AddInVar(IOVars,"tRes") call AddInVar(IOVars,"isFirstCpl") + call AddInVar(IOVars,"cs_doCS_next_preAdv") + call AddInVar(IOVars,"cs_modelDst_next_preAdv") call AddInVar(IOVars,"cs_doneFirstCS") call AddInVar(IOVars,"cs_lastEval") call AddInVar(IOVars,"cs_lastTarget") @@ -651,12 +659,15 @@ module raijuIO ! Coldstarter State%coldStarter%lastEval = GetIOReal(IOVars, "cs_lastEval") State%coldStarter%lastTarget = GetIOReal(IOVars, "cs_lastTarget") + State%coldStarter%modelDst_next_preAdv = GetIOReal(IOVars, "cs_modelDst_next_preAdv") ! Handle boolean attributes tmpInt = GetIOInt(IOVars, "isFirstCpl") State%isFirstCpl = tmpInt .eq. 1 tmpInt = GetIOInt(IOVars, "cs_doneFirstCS") State%coldStarter%doneFirstCS = tmpInt .eq. 1 + tmpInt = GetIOInt(IOVars, "cs_doCS_next_preAdv") + State%coldStarter%doCS_next_preAdv = tmpInt .eq. 1 call IOArray2DFill(IOVars, "xmin" , State%xyzMin(:,:,XDIR)) call IOArray2DFill(IOVars, "ymin" , State%xyzMin(:,:,YDIR)) diff --git a/src/raiju/raijuOut.F90 b/src/raiju/raijuOut.F90 index 3133ac15..ebdfa689 100644 --- a/src/raiju/raijuOut.F90 +++ b/src/raiju/raijuOut.F90 @@ -131,7 +131,7 @@ module raijuOut if (maxP_MLT > 24) maxP_MLT = maxP_MLT - 24D0 write(*,'(a,I0,a,f7.2,a,f7.2,a,f5.2,a,f5.2,a,f7.2)') ' ', & Grid%spc(s)%flav, ': P =', maxPress,', D =',maxDen,' @ ',maxP_L,' Rp,',maxP_MLT, & - " MLT; DPS:",spcEta2DPS(Model, Grid, State, Grid%spc(sIdx), State%active .eq. RAIJUACTIVE) + " MLT; DPS:",spcEta2DPS(Model, Grid, State%bvol_cc, State%eta_avg, Grid%spc(sIdx), State%active .eq. RAIJUACTIVE) enddo write(*,'(a)',advance="no") ANSIRESET diff --git a/src/raiju/raijuPreAdvancer.F90 b/src/raiju/raijuPreAdvancer.F90 index 3d1ff0a8..44bb54a1 100644 --- a/src/raiju/raijuPreAdvancer.F90 +++ b/src/raiju/raijuPreAdvancer.F90 @@ -36,6 +36,9 @@ module raijuPreAdvancer call fillArray(State%eta_avg, 0.0_rp) ! (losses handled in updateRaiLosses) + ! Now that topo is set, we can calculate active domain + call setActiveDomain(Model, Grid, State) + ! Moments to etas, initial active shell calculation call Tic("BCs") call applyRaijuBCs(Model, Grid, State, doWholeDomainO=State%isFirstCpl) ! If fullEtaMap=True, mom2eta map is applied to the whole domain diff --git a/src/voltron/modelInterfaces/raiCplTypesSub.F90 b/src/voltron/modelInterfaces/raiCplTypesSub.F90 index 4a3f870b..fbda3c80 100644 --- a/src/voltron/modelInterfaces/raiCplTypesSub.F90 +++ b/src/voltron/modelInterfaces/raiCplTypesSub.F90 @@ -43,11 +43,7 @@ submodule (volttypes) raijuCplTypesSub ! 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 - if (App%opt%doColdStart) then - ! We are gonna cold start, so ignore plasma ingestion rules for first coupling - App%raiApp%State%isFirstCpl = .false. - endif + ! Then allocate and initialize coupling variables based on raiju app call raijuCpl_init(App, xml) @@ -58,30 +54,12 @@ submodule (volttypes) raijuCplTypesSub class(raijuCoupler_T), intent(inout) :: App class(voltApp_T), intent(inout) :: vApp - logical :: doFirstColdStart - logical :: doUpdateColdStart real(rp) :: BSDst - doFirstColdStart = .false. - doUpdateColdStart = .false. - associate(raiApp=>App%raiApp) ! If we are running realtime, its our job to get everything we need from vApp into raiCpl if (.not. App%raiApp%Model%isSA) then - ! First, determine if we should cold start, i.e. Completely reset raiju's eta's to match some target conditions - ! Determine if we should cold start before packing coupler because it will set tLastUpdate to vApp%time and then we can't do the checks we want - ! But actually do cold start after coupler packing completes so we can use real field line info - - ! Do we do our very first coldstart ever - if (App%opt%doColdStart .and. App%tLastUpdate < 0.0 .and. vApp%time >= 0.0) then - doFirstColdStart = .true. - endif - ! Do we do "updates" to our coldstart during pre-conditioning period - if(App%opt%doColdStart .and. App%tLastUpdate > 0.0 .and. vApp%time < App%startup_blendTscl) then - doUpdateColdStart = .true. - endif - call packRaijuCoupler_RT(App, vApp) endif @@ -90,29 +68,14 @@ submodule (volttypes) raijuCplTypesSub if (.not. raiApp%State%coldStarter%doneFirstCS .or. vApp%time < raiApp%State%coldStarter%tEnd) then !! Make sure we run at least once - call setActiveDomain(raiApp%Model, raiApp%Grid, raiApp%State) ! Calc voltron dst ourselves since vApp%BSDst is only set on console output call EstDST(vApp%gApp%Model,vApp%gApp%Grid,vApp%gApp%State,BSDst0=BSDst) - call raijuGeoColdStart(raiApp%Model, raiApp%Grid, raiApp%State, vApp%time, BSDst) + raiApp%State%coldStarter%doCS_next_preAdv = .true. + raiApp%State%coldStarter%modelDst_next_preAdv = BSDst + !call setActiveDomain(raiApp%Model, raiApp%Grid, raiApp%State) + !call raijuGeoColdStart(raiApp%Model, raiApp%Grid, raiApp%State, vApp%time, BSDst) endif - !if (doFirstColdStart) then - ! ! Its happening, everybody stay calm - ! write(*,*) "RAIJU Doing first cold start..." - ! ! NOTE: By this point we have put coupling info into raiju (e.g. bVol, xyzmin, MHD moments) - ! ! But haven't calculated active domain yet because that happens in preadvancer - ! ! So we jump in and do it here so we have it for cold starting - ! call setActiveDomain(raiApp%Model, raiApp%Grid, raiApp%State) - ! ! Calc voltron dst ourselves since vApp%BSDst is only set on console output - ! call EstDST(vApp%gApp%Model,vApp%gApp%Grid,vApp%gApp%State,BSDst0=BSDst) - ! call raijuGeoColdStart(raiApp%Model, raiApp%Grid, raiApp%State, vApp%time, BSDst, doCXO=App%doColdstartCX,doPsphO=.true.) - !endif - !if (doUpdateColdStart) then - ! write(*,*)"RAIJU doing update cold start at t=",vApp%time - ! write(*,*)" (calculating model BSDst,)",vApp%time - ! call setActiveDomain(raiApp%Model, raiApp%Grid, raiApp%State) - ! call EstDST(vApp%gApp%Model,vApp%gApp%Grid,vApp%gApp%State,BSDst0=BSDst) - ! call raijuGeoColdStart(raiApp%Model, raiApp%Grid, raiApp%State, vApp%time, BSDst, doCXO=App%doColdstartCX,doPsphO=.false.) - !endif + end associate end subroutine volt2RAIJU From a46c66a774c9b01680242237a1def240aa197476 Mon Sep 17 00:00:00 2001 From: Anthony Date: Fri, 5 Sep 2025 12:45:15 -0600 Subject: [PATCH 03/15] Adding more reasonable default grid res options for voltron based on gamera grid size --- src/voltron/voltapp.F90 | 33 +++++++++++++++++++++++++++++---- 1 file changed, 29 insertions(+), 4 deletions(-) diff --git a/src/voltron/voltapp.F90 b/src/voltron/voltapp.F90 index b87e3ea0..7232ea23 100644 --- a/src/voltron/voltapp.F90 +++ b/src/voltron/voltapp.F90 @@ -209,7 +209,7 @@ module voltapp gApp%Model%t = vApp%time / gApp%Model%Units%gT0 gApp%State%time = gApp%Model%t - call genVoltShellGrid(vApp, xmlInp) + call genVoltShellGrid(vApp, xmlInp, gApp%Grid%Nkp) call initVoltState(vApp) endif @@ -754,11 +754,14 @@ module voltapp end subroutine init_volt2Chmp - subroutine genVoltShellGrid(vApp, xmlInp) + subroutine genVoltShellGrid(vApp, xmlInp, gamRes) class(voltApp_T) , intent(inout) :: vApp type(XML_Input_T), intent(in) :: xmlInp + integer, intent(in) :: gamRes character(len=strLen) :: gType + integer :: Nt_def, Np_def + !! Default number of active cells in theta and phi unless xml says otherwise integer :: Nt, Np !! Number of active cells in theta and phi integer :: Ng @@ -786,8 +789,30 @@ module voltapp ! Note: Nt is for a single hemisphere, we will manually double it in a minute ! TODO: This means we will always have even number of total cells, and a cell interfce right on the equator ! Can upgrade to allow for odd number later - call xmlInp%Set_Val(Nt, "grid/Nt", 180 ) ! 1 deg res default for uniform grid - call xmlInp%Set_Val(Np, "grid/Np", 360) ! 1 deg res default + + ! First determine defaults + if (gamRes<=64) then + ! DBL + Nt_def = 90 + Np_def = 180 + else if (gamRes<=128) then + ! QUAD + Nt_def = 180 + Np_def = 360 + else if (gamRes<=256) then + ! OCT + Nt_def = 360 + Np_def = 720 + else + ! HEX or above + ! Idk good luck + Nt_def = 540 + Np_def = 1440 + endif + + + call xmlInp%Set_Val(Nt, "grid/Nt", Nt_def) ! 1 deg res default for uniform grid + call xmlInp%Set_Val(Np, "grid/Np", Np_def) ! 1 deg res default ! Ghost cells call xmlInp%Set_Val(Ng, "grid/Ng", 4) ! # of ghosts in every direction nGhosts = 0 From 44370b4742b347e2950f04a2b154ab86c939dd99 Mon Sep 17 00:00:00 2001 From: Anthony Date: Wed, 10 Sep 2025 11:02:47 -0600 Subject: [PATCH 04/15] Debug messages and flag fix --- src/raiju/raijuBCs.F90 | 6 ++++++ src/raiju/raijuColdStartHelper.F90 | 7 ++++--- src/raiju/raijuPreAdvancer.F90 | 2 +- 3 files changed, 11 insertions(+), 4 deletions(-) diff --git a/src/raiju/raijuBCs.F90 b/src/raiju/raijuBCs.F90 index f4bd3acf..556dffe7 100644 --- a/src/raiju/raijuBCs.F90 +++ b/src/raiju/raijuBCs.F90 @@ -30,6 +30,12 @@ module raijuBCs doWholeDomain = .false. endif + if(doWholeDomain) then + write(*,*)"applyRaijuBCs: doWholeDomain=T" + else + write(*,*)"applyRaijuBCs: doWholeDomain=F" + endif + call calcMomentIngestionLocs(Model, Grid, State, doWholeDomain, doMomentIngest) call applyMomentIngestion(Model, Grid, State, doMomentIngest) diff --git a/src/raiju/raijuColdStartHelper.F90 b/src/raiju/raijuColdStartHelper.F90 index 43d28ccf..b8a85033 100644 --- a/src/raiju/raijuColdStartHelper.F90 +++ b/src/raiju/raijuColdStartHelper.F90 @@ -101,8 +101,11 @@ module raijuColdStartHelper cs%lastEval = t0 cs%lastTarget = dstTarget + cs%doneFirstCS = .true. ! Whether we do anything or not, we were at least called once - if (dstTarget > 0) then ! We got nothing to contribute + if (dstTarget > 0) then ! We've got nothing to contribute + write(*,*)"RAIJU coldstart not adding anything" + write(*,*)'doAccumulate=',doAccumulate return endif @@ -154,8 +157,6 @@ module raijuColdStartHelper State%eta = etaCS endif - State%coldStarter%doneFirstCS = .true. - end subroutine raijuGeoColdStart diff --git a/src/raiju/raijuPreAdvancer.F90 b/src/raiju/raijuPreAdvancer.F90 index 44bb54a1..ea5dede3 100644 --- a/src/raiju/raijuPreAdvancer.F90 +++ b/src/raiju/raijuPreAdvancer.F90 @@ -44,7 +44,7 @@ module raijuPreAdvancer call applyRaijuBCs(Model, Grid, State, doWholeDomainO=State%isFirstCpl) ! If fullEtaMap=True, mom2eta map is applied to the whole domain call Toc("BCs") - ! Handle plasmaasphere refilling for the full step about to happen + ! Handle plasmasphere refilling for the full step about to happen call plasmasphereRefill(Model,Grid,State) ! Handle edge cases that may effect the validity of information carried over from last coupling period From 30aab6296f83f9752ce8a9d37ce336eda2f0f7a5 Mon Sep 17 00:00:00 2001 From: "Anthony M. Sciola" Date: Wed, 10 Sep 2025 12:29:28 -0700 Subject: [PATCH 05/15] fixing (hopefully) some coldstart logic around when it runs and what it adds. removing debug output --- src/raiju/raijuBCs.F90 | 5 -- src/raiju/raijuColdStartHelper.F90 | 88 ++++++++++++++++--- .../modelInterfaces/raiCplTypesSub.F90 | 21 ++--- 3 files changed, 88 insertions(+), 26 deletions(-) diff --git a/src/raiju/raijuBCs.F90 b/src/raiju/raijuBCs.F90 index 556dffe7..cc412ad9 100644 --- a/src/raiju/raijuBCs.F90 +++ b/src/raiju/raijuBCs.F90 @@ -30,11 +30,6 @@ module raijuBCs doWholeDomain = .false. endif - if(doWholeDomain) then - write(*,*)"applyRaijuBCs: doWholeDomain=T" - else - write(*,*)"applyRaijuBCs: doWholeDomain=F" - endif call calcMomentIngestionLocs(Model, Grid, State, doWholeDomain, doMomentIngest) call applyMomentIngestion(Model, Grid, State, doMomentIngest) diff --git a/src/raiju/raijuColdStartHelper.F90 b/src/raiju/raijuColdStartHelper.F90 index b8a85033..fe270f20 100644 --- a/src/raiju/raijuColdStartHelper.F90 +++ b/src/raiju/raijuColdStartHelper.F90 @@ -56,7 +56,7 @@ module raijuColdStartHelper !! If false, replace State%eta with coldStart info !! Default: false - logical :: isFirstCS, doAccumulate + logical :: doInitRC, doAccumulate integer :: s, sIdx_p, sIdx_e real(rp) :: dstReal, dstTarget real(rp) :: dps_current, dps_preCX, dps_postCX, dps_rescale, dps_ele @@ -88,32 +88,33 @@ module raijuColdStartHelper if (.not. cs%doneFirstCS) then ! On first try, we assume there is no existing ring current, and its our job to make up the entire difference dstTarget = dstReal - dstModel - + doInitRC = .true. else if (t0 > (cs%lastEval + cs%evalCadence)) then ! If we are updating, there should already be some ring current ! If dstReal - dstModel is still < 0, we need to add ADDITIONAL pressure to get them to match dps_current = spcEta2DPS(Model, Grid, State%bvol_cc, State%eta, Grid%spc(sIdx_p), isGood) + spcEta2DPS(Model, Grid, State%bvol_cc, State%eta, Grid%spc(sIdx_e), isGood) dstTarget = dstReal - (dstModel - dps_current) - else - ! Otherwise we have nothing to do, just chill til next update time - return + doInitRC = .false. endif cs%lastEval = t0 cs%lastTarget = dstTarget cs%doneFirstCS = .true. ! Whether we do anything or not, we were at least called once + + ! Init psphere (will always run if this function is called) + call setRaijuInitPsphere(Model, Grid, State, Model%psphInitKp) - if (dstTarget > 0) then ! We've got nothing to contribute + ! Now decide if we need to add a starter ring current + if (dstTarget >= 0) then ! We've got nothing to contribute write(*,*)"RAIJU coldstart not adding anything" write(*,*)'doAccumulate=',doAccumulate return endif - if (.not. cs%doneFirstCS) then - ! Init psphere (doesn't care about accumulation, always hard resets) - call setRaijuInitPsphere(Model, Grid, State, Model%psphInitKp) + if (doInitRC) then ! Init hot protons - call raiColdStart_initHOTP(Model, Grid, State, t0, dstTarget, etaCS) + !call raiColdStart_initHOTP(Model, Grid, State, t0, dstTarget, etaCS) + call raiColdStart_initHOTP_RCOnly(Model, Grid, State, t0, dstTarget, etaCS) dps_preCX = spcEta2DPS(Model, Grid, State%bvol_cc, etaCS, Grid%spc(sIdx_p), isGood) ! Hit it with some charge exchange if (cs%doCX) then @@ -130,7 +131,7 @@ module raijuColdStartHelper etaCS(:,:,Grid%spc(sIdx_p)%kStart:Grid%spc(sIdx_p)%kEnd) = etaScale*etaCS(:,:,Grid%spc(sIdx_p)%kStart:Grid%spc(sIdx_p)%kEnd) dps_rescale = spcEta2DPS(Model, Grid, State%bvol_cc, etaCS, Grid%spc(sIdx_p), isGood) - if (isfirstCS) then + if (doInitRC) then write(*,*) "RAIJU Cold starting..." write(*,'(a,f7.2)') " Real Dst : ",dstReal write(*,'(a,f7.2)') " Model Dst : ",dstModel @@ -160,6 +161,71 @@ module raijuColdStartHelper end subroutine raijuGeoColdStart + subroutine raiColdStart_initHOTP_RCOnly(Model, Grid, State, t0, dstTarget, etaCS) + type(raijuModel_T), intent(in) :: Model + type(raijuGrid_T), intent(in) :: Grid + type(raijuState_T), intent(in) :: State + real(rp), intent(in) :: t0 + !! Target time to pull SW values from + real(rp), intent(in) :: dstTarget + real(rp), dimension(Grid%shGrid%isg:Grid%shGrid%ieg, Grid%shGrid%jsg:Grid%shGrid%jeg, Grid%Nk), intent(inout) :: etaCS + + real(rp) :: dstTarget_p + logical :: isInTM03 + integer :: i,j,sIdx + integer, dimension(2) :: ij_TM + real(rp) :: vSW, dSW, dPS_emp, pPS_emp, ktPS_emp + real(rp) :: x0_TM, y0_TM, T0_TM, Bvol0_TM, P0_ps, N0_ps + real(rp) :: L, vm, P_rc, D_rc, kt_rc + + sIdx = spcIdx(Grid, F_HOTP) + associate(sh=>Grid%shGrid, spc=>Grid%spc(sIdx)) + + ! Set everything to zero to start + etaCS(:,:,spc%kStart:spc%kEnd) = 0.0_rp + + ! Scale target Dst down to account for electrons contributing stuff later + dstTarget_p = dstTarget / (1.0 + 1.0/Model%tiote) + call SetQTRC(dstTarget_p,doVerbO=.false.) ! This sets a global QTRC_P0 inside earthhelper.F90 + + + ! Get reference TM value at -10 Re + x0_TM = -10.0-TINY + y0_TM = 0.0 + ! Empirical temperature + call EvalTM03([x0_TM,y0_TM,0.0_rp],N0_ps,P0_ps,isInTM03) + T0_TM = DP2kT(N0_ps, P0_ps) + ! Model FTV + ij_TM = minloc( sqrt( (State%xyzMincc(:,:,XDIR)-x0_TM)**2 + (State%xyzMincc(:,:,YDIR)**2) ) ) + Bvol0_TM = State%bvol_cc(ij_TM(IDIR), ij_TM(JDIR)) + + ! Now set our initial density and pressure profile + do j=sh%jsg,sh%jeg + do i=sh%isg,sh%ie ! Note: Not setting low lat ghosts, we want them to be zero + + if (State%active(i,j) .eq. RAIJUINACTIVE) cycle + + L = norm2(State%xyzMincc(i,j,XDIR:YDIR)) + vm = State%bvol_cc(i,j)**(-2./3.) + + kt_rc = T0_TM*(Bvol0_TM/State%bvol_cc(i,j))**(2./3.) + kt_rc = min(kt_rc, 4.0*T0_TM) ! Limit cap. Not a big fan, but without cap we get stuff that's too energetic and won't go away (until FLC maybe) + + + P_rc = P_QTRC(L) ! From earthhelper.F90 + D_rc = PkT2Den(P_rc, kt_rc) + + + ! Finally map it to HOTP etas + call DkT2SpcEta(Model, spc, etaCS(i,j,spc%kStart:spc%kEnd), D_rc, kt_rc, vm) + + enddo + enddo + + end associate + + end subroutine raiColdStart_initHOTP_RCOnly + subroutine raiColdStart_initHOTP(Model, Grid, State, t0, dstTarget, etaCS) type(raijuModel_T), intent(in) :: Model type(raijuGrid_T), intent(in) :: Grid diff --git a/src/voltron/modelInterfaces/raiCplTypesSub.F90 b/src/voltron/modelInterfaces/raiCplTypesSub.F90 index fbda3c80..03a4d88e 100644 --- a/src/voltron/modelInterfaces/raiCplTypesSub.F90 +++ b/src/voltron/modelInterfaces/raiCplTypesSub.F90 @@ -66,16 +66,17 @@ submodule (volttypes) raijuCplTypesSub ! Someone updated raiCpl's coupling variables by now, stuff it into RAIJU proper call raiCpl2RAIJU(App) - if (.not. raiApp%State%coldStarter%doneFirstCS .or. vApp%time < raiApp%State%coldStarter%tEnd) then - !! Make sure we run at least once - ! Calc voltron dst ourselves since vApp%BSDst is only set on console output - call EstDST(vApp%gApp%Model,vApp%gApp%Grid,vApp%gApp%State,BSDst0=BSDst) - raiApp%State%coldStarter%doCS_next_preAdv = .true. - raiApp%State%coldStarter%modelDst_next_preAdv = BSDst - !call setActiveDomain(raiApp%Model, raiApp%Grid, raiApp%State) - !call raijuGeoColdStart(raiApp%Model, raiApp%Grid, raiApp%State, vApp%time, BSDst) - endif - + associate(cs=>raiApp%State%coldStarter) + if (.not. cs%doneFirstCS .or. (cs%doUpdate .and. vApp%time < cs%tEnd) ) then + !! Make sure we run at least once + ! Calc voltron dst ourselves since vApp%BSDst is only set on console output + call EstDST(vApp%gApp%Model,vApp%gApp%Grid,vApp%gApp%State,BSDst0=BSDst) + raiApp%State%coldStarter%doCS_next_preAdv = .true. + raiApp%State%coldStarter%modelDst_next_preAdv = BSDst + !call setActiveDomain(raiApp%Model, raiApp%Grid, raiApp%State) + !call raijuGeoColdStart(raiApp%Model, raiApp%Grid, raiApp%State, vApp%time, BSDst) + endif + end associate end associate end subroutine volt2RAIJU From 59cc61037611e4ebf6a454d1ebf36f4de6e67c6d Mon Sep 17 00:00:00 2001 From: "Anthony M. Sciola" Date: Thu, 11 Sep 2025 16:21:11 -0700 Subject: [PATCH 06/15] More improvements to coldstart behavior. Fixing issue with gradVM calc that can lead to negative bvol. Couple other improvements just in case --- src/raiju/raijuBCs.F90 | 20 +++++--- src/raiju/raijuColdStartHelper.F90 | 75 ++++++++++++++++++++---------- src/raiju/raijuIO.F90 | 2 +- src/raiju/raijuPreAdvancer.F90 | 10 ++-- src/raiju/raijuStarter.F90 | 2 + 5 files changed, 71 insertions(+), 38 deletions(-) diff --git a/src/raiju/raijuBCs.F90 b/src/raiju/raijuBCs.F90 index cc412ad9..84b41258 100644 --- a/src/raiju/raijuBCs.F90 +++ b/src/raiju/raijuBCs.F90 @@ -68,9 +68,18 @@ module raijuBCs doMomentIngest = .false. ! Determine where to do BCs if(doWholeDomain) then - where (State%active .ne. RAIJUINACTIVE) - doMomentIngest = .true. - end where + !where (State%active .ne. RAIJUINACTIVE) + ! doMomentIngest = .true. + !end where + associate(sh=>Grid%shGrid) + do j=sh%jsg,sh%jeg + do i=sh%isg,sh%ieg + if (State%active(i,j) .ne. RAIJUINACTIVE) then + doMomentIngest(i,j) = .true. + endif + enddo + enddo + end associate else associate(sh=>Grid%shGrid) @@ -122,7 +131,6 @@ module raijuBCs psphIdx = spcIdx(Grid, F_PSPH) eleIdx = spcIdx(Grid, F_HOTE) !$OMP PARALLEL DO default(shared) & - !$OMP schedule(dynamic) & !$OMP private(i,j,s,fIdx,fm,vm,kT,etaBelow,tmp_kti,tmp_kte,eMin,tmp_D,tmp_P) do j=Grid%shGrid%jsg,Grid%shGrid%jeg do i=Grid%shGrid%isg,Grid%shGrid%ieg @@ -219,8 +227,8 @@ module raijuBCs press2D = State%Pavg (i-2:i+2, j-2:j+2, fIdx) wgt2D = State%domWeights(i-2:i+2, j-2:j+2) isG2D = State%active (i-2:i+2, j-2:j+2) .ne. RAIJUINACTIVE - D = sum(den2D * wgt2D, mask=isG2D)/sum(wgt2D, mask=isG2D) - P = sum(press2D* wgt2D, mask=isG2D)/sum(wgt2D, mask=isG2D) + D = sum(den2D * wgt2D, mask=isG2D)/max( sum(wgt2D, mask=isG2D), TINY) + P = sum(press2D* wgt2D, mask=isG2D)/max( sum(wgt2D, mask=isG2D), TINY) endif end subroutine getDomWeightedMoments diff --git a/src/raiju/raijuColdStartHelper.F90 b/src/raiju/raijuColdStartHelper.F90 index fe270f20..b98781d9 100644 --- a/src/raiju/raijuColdStartHelper.F90 +++ b/src/raiju/raijuColdStartHelper.F90 @@ -33,6 +33,8 @@ module raijuColdStartHelper call iXML%Set_Val(coldStarter%tEnd,'coldStarter/tEnd',coldStarter%evalCadence-TINY) ! Don't do any updates as default endif + coldStarter%doneFirstCS = .false. + end subroutine initRaijuColdStarter @@ -57,6 +59,7 @@ module raijuColdStartHelper !! Default: false logical :: doInitRC, doAccumulate + integer :: i,j,k integer :: s, sIdx_p, sIdx_e real(rp) :: dstReal, dstTarget real(rp) :: dps_current, dps_preCX, dps_postCX, dps_rescale, dps_ele @@ -74,7 +77,7 @@ module raijuColdStartHelper associate(cs=>State%coldStarter) - where (State%active .eq. RAIJUACTIVE) + where (State%active .ne. RAIJUINACTIVE) isGood = .true. elsewhere isGood = .false. @@ -105,16 +108,16 @@ module raijuColdStartHelper call setRaijuInitPsphere(Model, Grid, State, Model%psphInitKp) ! Now decide if we need to add a starter ring current - if (dstTarget >= 0) then ! We've got nothing to contribute - write(*,*)"RAIJU coldstart not adding anything" - write(*,*)'doAccumulate=',doAccumulate - return - endif + !if (dstTarget >= 0) then ! We've got nothing to contribute + ! write(*,*)"RAIJU coldstart not adding anything" + ! write(*,*)'doAccumulate=',doAccumulate + ! return + !endif if (doInitRC) then ! Init hot protons - !call raiColdStart_initHOTP(Model, Grid, State, t0, dstTarget, etaCS) - call raiColdStart_initHOTP_RCOnly(Model, Grid, State, t0, dstTarget, etaCS) + call raiColdStart_initHOTP(Model, Grid, State, t0, dstTarget, etaCS) + !call raiColdStart_initHOTP_RCOnly(Model, Grid, State, t0, dstTarget, etaCS) dps_preCX = spcEta2DPS(Model, Grid, State%bvol_cc, etaCS, Grid%spc(sIdx_p), isGood) ! Hit it with some charge exchange if (cs%doCX) then @@ -127,19 +130,23 @@ module raijuColdStartHelper dps_current = dps_postCX ! Note: if using fudge we're gonna lose electrons immediately, don't include them in current dst for now endif - etaScale = abs(dstTarget / dps_current) - etaCS(:,:,Grid%spc(sIdx_p)%kStart:Grid%spc(sIdx_p)%kEnd) = etaScale*etaCS(:,:,Grid%spc(sIdx_p)%kStart:Grid%spc(sIdx_p)%kEnd) - dps_rescale = spcEta2DPS(Model, Grid, State%bvol_cc, etaCS, Grid%spc(sIdx_p), isGood) + if (dstTarget < 0) then + etaScale = abs(dstTarget / dps_current) + etaCS(:,:,Grid%spc(sIdx_p)%kStart:Grid%spc(sIdx_p)%kEnd) = etaScale*etaCS(:,:,Grid%spc(sIdx_p)%kStart:Grid%spc(sIdx_p)%kEnd) + dps_rescale = spcEta2DPS(Model, Grid, State%bvol_cc, etaCS, Grid%spc(sIdx_p), isGood) + else + dps_rescale = dps_current + endif if (doInitRC) then write(*,*) "RAIJU Cold starting..." - write(*,'(a,f7.2)') " Real Dst : ",dstReal - write(*,'(a,f7.2)') " Model Dst : ",dstModel - write(*,'(a,f7.2)') " Target DPS-Dst : ",dstTarget - write(*,'(a,f7.2)') " Hot proton pre-loss : ",dps_preCX - write(*,'(a,f7.2)') " post-loss : ",dps_postCX - write(*,'(a,f7.2)') " post-rescale : ",dps_rescale - write(*,'(a,f7.2)') " Hot electron DPS-Dst : ",dps_ele + write(*,'(a,f10.6)') " Real Dst : ",dstReal + write(*,'(a,f10.6)') " Model Dst : ",dstModel + write(*,'(a,f10.6)') " Target DPS-Dst : ",dstTarget + write(*,'(a,f10.6)') " Hot proton pre-loss : ",dps_preCX + write(*,'(a,f10.6)') " post-loss : ",dps_postCX + write(*,'(a,f10.6)') " post-rescale : ",dps_rescale + write(*,'(a,f10.6)') " Hot electron DPS-Dst : ",dps_ele else write(*,'(a,f7.2)') " Real Dst : ",dstReal write(*,'(a,f7.2)') " Model Dst : ",dstModel @@ -153,7 +160,20 @@ module raijuColdStartHelper ! finally, put it into raiju state if(doAccumulate) then - State%eta = State%eta + etaCS + !State%eta = State%eta + etaCS + associate(sh=>Grid%shGrid) + !$OMP PARALLEL DO default(shared) & + !$OMP private(i,j,k) + do k=1,Grid%Nk + do j=sh%jsg,sh%jeg + do i=sh%isg,sh%ieg + if (etaCS(i,j,k) > State%eta(i,j,k)) then + State%eta(i,j,k) = etaCS(i,j,k) + endif + enddo + enddo + enddo + end associate else State%eta = etaCS endif @@ -253,8 +273,10 @@ module raijuColdStartHelper call InitTM03(Model%tsF,t0) ! Scale target Dst down to account for electrons contributing stuff later - dstTarget_p = dstTarget / (1.0 + 1.0/Model%tiote) - call SetQTRC(dstTarget_p,doVerbO=.false.) ! This sets a global QTRC_P0 inside earthhelper.F90 + if (dstTarget < 0) then + dstTarget_p = dstTarget / (1.0 + 1.0/Model%tiote) + call SetQTRC(dstTarget_p,doVerbO=.false.) ! This sets a global QTRC_P0 inside earthhelper.F90 + endif ! Get Borovsky statistical values vSW = abs(GetSWVal("Vx",Model%tsF,t0)) @@ -295,7 +317,11 @@ module raijuColdStartHelper call EvalTM03_SM(State%xyzMincc(i,j,:),N0_ps,P0_ps,isInTM03) - P_rc = P_QTRC(L) ! From earthhelper.F90 + if (dstTarget < 0) then + P_rc = P_QTRC(L) ! From earthhelper.F90 + else + P_rc = 0.0 + endif if (.not. isInTM03) then N0_ps = dPS_emp @@ -337,8 +363,7 @@ module raijuColdStartHelper call lossCX%doInit(Model, Grid, nullXML) !$OMP PARALLEL DO default(shared) & - !$OMP schedule(dynamic) & - !$OMP private(i,j,tau) + !$OMP private(i,j,k,tau) do j=Grid%shGrid%jsg,Grid%shGrid%jeg do i=Grid%shGrid%isg,Grid%shGrid%ieg do k = spc%kStart,spc%kEnd @@ -372,7 +397,7 @@ module raijuColdStartHelper !$OMP PARALLEL DO default(shared) & !$OMP schedule(dynamic) & - !$OMP private(i,j,vm,den,kt_p,kt_e) + !$OMP private(i,j,vm,den_p, press_p,kt_p,kt_e) do j=Grid%shGrid%jsg,Grid%shGrid%jeg do i=Grid%shGrid%isg,Grid%shGrid%ie ! Note: Not setting low lat ghosts, we want them to be zero if (State%active(i,j) .eq. RAIJUINACTIVE) cycle diff --git a/src/raiju/raijuIO.F90 b/src/raiju/raijuIO.F90 index 3fe18b36..f4d03529 100644 --- a/src/raiju/raijuIO.F90 +++ b/src/raiju/raijuIO.F90 @@ -16,7 +16,7 @@ module raijuIO implicit none - integer, parameter, private :: MAXIOVAR = 70 + integer, parameter, private :: MAXIOVAR = 100 logical, private :: doRoot = .true. !Whether root variables need to be written logical, private :: doFat = .false. !Whether to output lots of extra datalogical, private :: doRoot = .true. !Whether root variables need to be written diff --git a/src/raiju/raijuPreAdvancer.F90 b/src/raiju/raijuPreAdvancer.F90 index ea5dede3..a320a4b8 100644 --- a/src/raiju/raijuPreAdvancer.F90 +++ b/src/raiju/raijuPreAdvancer.F90 @@ -307,7 +307,6 @@ module raijuPreAdvancer associate(sh=>Grid%shGrid) !$OMP PARALLEL DO default(shared) & - !$OMP schedule(dynamic) & !$OMP private(i,j,qLow,qHigh) do j=sh%jsg,sh%jeg do i=sh%isg,sh%ieg @@ -335,7 +334,6 @@ module raijuPreAdvancer allocate(gradQtmp(sh%isg:sh%ieg,sh%jsg:sh%jeg, 2)) gradQtmp = gradQ !$OMP PARALLEL DO default(shared) & - !$OMP schedule(dynamic) & !$OMP private(i,j) do j=sh%jsg+1,sh%jeg-1 do i=sh%isg+1,sh%ieg-1 @@ -486,12 +484,11 @@ module raijuPreAdvancer gradVM(:,:,RAI_TH) = gradVM(:,:,RAI_TH) + dV0_dth !$OMP PARALLEL DO default(shared) & - !$OMP schedule(dynamic) & !$OMP private(i,j,bVolcc) do j=sh%jsg,sh%jeg do i=sh%isg,sh%ieg - bVolcc = toCenter2D(dV(i:i+1,j:j+1)) + DipFTV_colat(Grid%thcRp(i), B0) ! Will include smoothing of dV if enabled - !bVolcc = toCenter2D(V(i:i+1,j:j+1)) + !bVolcc = toCenter2D(dV(i:i+1,j:j+1)) + DipFTV_colat(Grid%thcRp(i), B0) ! Will include smoothing of dV if enabled + bVolcc = toCenter2D(V(i:i+1,j:j+1)) gradVM(i,j,:) = (-2./3.)*bVolcc**(-5./3.)*gradVM(i,j,:) enddo enddo @@ -567,7 +564,6 @@ module raijuPreAdvancer enddo ! Now everyone else !$OMP PARALLEL DO default(shared) & - !$OMP schedule(dynamic) & !$OMP private(i,j) do j=jsg+1,jeg do i=isg+1,ieg @@ -575,6 +571,8 @@ module raijuPreAdvancer enddo enddo + call wrapJcorners(sh, Vsm) + ! Write back to provided array V = Vsm end associate diff --git a/src/raiju/raijuStarter.F90 b/src/raiju/raijuStarter.F90 index 31eccc3d..7d174bdf 100644 --- a/src/raiju/raijuStarter.F90 +++ b/src/raiju/raijuStarter.F90 @@ -553,6 +553,8 @@ module raijustarter ! Similarly, set vaFrac to safe value in case stand-alone never writes to it State%vaFrac = 1.0 + State%isFirstCpl = .true. + ! Init State sub-modules if (Model%isSA) then ! If we are standalone, this is the place to get coldStarter settings From f7e4b2d062f0e31f42ffaa1c89826c42d327ae30 Mon Sep 17 00:00:00 2001 From: "Anthony M. Sciola" Date: Thu, 11 Sep 2025 17:04:51 -0700 Subject: [PATCH 07/15] Exposing option to smooth bvol and electric potential before gradient calculation in xml --- src/base/defs/raijudefs.F90 | 1 + src/base/types/raijuTypes.F90 | 2 ++ src/raiju/raijuPreAdvancer.F90 | 4 ++-- src/raiju/raijuStarter.F90 | 1 + 4 files changed, 6 insertions(+), 2 deletions(-) diff --git a/src/base/defs/raijudefs.F90 b/src/base/defs/raijudefs.F90 index 051cb799..d73d4122 100644 --- a/src/base/defs/raijudefs.F90 +++ b/src/base/defs/raijudefs.F90 @@ -70,6 +70,7 @@ module raijudefs real(rp), parameter :: def_cfl = 0.3 real(rp), parameter :: cflMax = 0.3 logical, parameter :: def_doUseVelLRs = .true. + logical, parameter :: def_doSmoothGrads = .true. ! Domain limits ! Buffer not allowed beyond min of maxTail and maxSun diff --git a/src/base/types/raijuTypes.F90 b/src/base/types/raijuTypes.F90 index 9a35b0fc..1f576138 100644 --- a/src/base/types/raijuTypes.F90 +++ b/src/base/types/raijuTypes.F90 @@ -226,6 +226,8 @@ module raijutypes !! For debug logical :: writeGhosts !! For debug + logical :: doSmoothGrads + !! Whether or not we smooth variables (bvol and electric potential) before taking gradients logical :: doClockConsoleOut !! If we are driving, output clock info logical :: doOutput_potGrads diff --git a/src/raiju/raijuPreAdvancer.F90 b/src/raiju/raijuPreAdvancer.F90 index a320a4b8..36819d69 100644 --- a/src/raiju/raijuPreAdvancer.F90 +++ b/src/raiju/raijuPreAdvancer.F90 @@ -264,7 +264,7 @@ module raijuPreAdvancer associate(sh=>Grid%shGrid) ! Gauss-Green calculation of cell-averaged gradients - call potExB(Grid%shGrid, State, pExB, doSmoothO=.true., isGCornerO=isGCorner) ! [V] + call potExB(Grid%shGrid, State, pExB, doSmoothO=Model%doSmoothGrads, isGCornerO=isGCorner) ! [V] call potCorot(Model%planet, Grid%shGrid, pCorot, Model%doGeoCorot) ! [V] call calcGradIJ_cc(Model%planet%rp_m, Grid, isGCorner, pExB , State%gradPotE_cc , doLimO=.true. ) ! [V/m] call calcGradIJ_cc(Model%planet%rp_m, Grid, isGCorner, pCorot, State%gradPotCorot_cc, doLimO=.false.) ! [V/m] @@ -273,7 +273,7 @@ module raijuPreAdvancer ! lambda is constant, so just need grad(V^(-2/3) ) call calcGradVM_cc(Model%planet%rp_m, Model%planet%ri_m, Model%planet%magMoment, & Grid, isGCorner, State%bvol, State%gradVM_cc, & - doSmoothO=.true., doLimO=.true.) + doSmoothO=Model%doSmoothGrads, doLimO=.true.) end associate end subroutine calcPotGrads_cc diff --git a/src/raiju/raijuStarter.F90 b/src/raiju/raijuStarter.F90 index 7d174bdf..54f9d5fe 100644 --- a/src/raiju/raijuStarter.F90 +++ b/src/raiju/raijuStarter.F90 @@ -175,6 +175,7 @@ module raijustarter Model%activeDomRad = abs(Model%activeDomRad) !---Solver ---! + call iXML%Set_Val(Model%doSmoothGrads,'sim/doSmoothGrads',def_doSmoothGrads) call iXML%Set_Val(Model%doUseVelLRs,'sim/useVelLRs',def_doUseVelLRs) call iXML%Set_Val(Model%maxItersPerSec,'sim/maxIter',def_maxItersPerSec) call iXML%Set_Val(Model%maxOrder,'sim/maxOrder',7) From 49425f09893d8a5c27401773189c294dac256c99 Mon Sep 17 00:00:00 2001 From: "Anthony M. Sciola" Date: Fri, 12 Sep 2025 12:13:11 -0700 Subject: [PATCH 08/15] Fixing DPS calc bug and resetting coldstart output format --- src/raiju/raijuColdStartHelper.F90 | 16 ++++++++-------- src/raiju/raijuEtautils.F90 | 6 +++--- 2 files changed, 11 insertions(+), 11 deletions(-) diff --git a/src/raiju/raijuColdStartHelper.F90 b/src/raiju/raijuColdStartHelper.F90 index b98781d9..a57aa450 100644 --- a/src/raiju/raijuColdStartHelper.F90 +++ b/src/raiju/raijuColdStartHelper.F90 @@ -140,15 +140,15 @@ module raijuColdStartHelper if (doInitRC) then write(*,*) "RAIJU Cold starting..." - write(*,'(a,f10.6)') " Real Dst : ",dstReal - write(*,'(a,f10.6)') " Model Dst : ",dstModel - write(*,'(a,f10.6)') " Target DPS-Dst : ",dstTarget - write(*,'(a,f10.6)') " Hot proton pre-loss : ",dps_preCX - write(*,'(a,f10.6)') " post-loss : ",dps_postCX - write(*,'(a,f10.6)') " post-rescale : ",dps_rescale - write(*,'(a,f10.6)') " Hot electron DPS-Dst : ",dps_ele - else write(*,'(a,f7.2)') " Real Dst : ",dstReal + write(*,'(a,f7.2)') " Model Dst : ",dstModel + write(*,'(a,f7.2)') " Target DPS-Dst : ",dstTarget + write(*,'(a,f7.2)') " Hot proton pre-loss : ",dps_preCX + write(*,'(a,f7.2)') " post-loss : ",dps_postCX + write(*,'(a,f7.2)') " post-rescale : ",dps_rescale + write(*,'(a,f7.2)') " Hot electron DPS-Dst : ",dps_ele + el7 +2 write(*,'(a,f7.2)') " Real Dst : ",dstReal write(*,'(a,f7.2)') " Model Dst : ",dstModel write(*,'(a,f7.2)') " Current DPS-Dst : ",dps_current write(*,'(a,f7.2)') " Target DPS-Dst : ",dstTarget diff --git a/src/raiju/raijuEtautils.F90 b/src/raiju/raijuEtautils.F90 index 3c6d248e..65a3fef7 100644 --- a/src/raiju/raijuEtautils.F90 +++ b/src/raiju/raijuEtautils.F90 @@ -189,8 +189,8 @@ module raijuetautils !! Eval mask, true = point is included in calculation real(rp) :: dpsdst - integer :: i,j,k - real(rp) :: press, bVol, energyDen, energy + integer :: i,j + real(rp) :: press, energyDen, energy logical :: isDead = .false. dpsdst = 0.0 @@ -199,7 +199,7 @@ module raijuetautils do i=Grid%shGrid%isg,Grid%shGrid%ieg if (.not. isGood(i,j)) cycle press = SpcEta2Press(spc, eta(i,j,spc%kStart:spc%kEnd), bvol_cc(i,j)) ! [nPa] - energyDen = (press*1.0D-9) * (bVol*Model%planet%rp_m*1.0D9) * (Grid%Brcc(i,j)*1.0D-9)/kev2J ! p[J/m^3] * bVol[m/T] * B[T] = [J/m^2] * keV/J = [keV/m^2] + energyDen = (press*1.0D-9) * (bVol_cc(i,j)*Model%planet%rp_m*1.0D9) * (Grid%Brcc(i,j)*1.0D-9)/kev2J ! p[J/m^3] * bVol[m/T] * B[T] = [J/m^2] * keV/J = [keV/m^2] energy = energyDen*(Grid%areaCC(i,j)*Model%planet%ri_m**2) ! [keV/m^2]* Re^2[m^2] = [keV] dpsdst = dpsdst - 4.2*(1.0D-30)*energy ! [nT] enddo From 0dd7c427698f437cf244e00444529619b8d7f0a5 Mon Sep 17 00:00:00 2001 From: "Anthony M. Sciola" Date: Fri, 12 Sep 2025 16:21:44 -0700 Subject: [PATCH 09/15] fixing weird typo, turning vApp verbosity back on after gam init (ty Jeff) --- src/raiju/raijuColdStartHelper.F90 | 3 +-- src/voltron/voltapp.F90 | 2 ++ 2 files changed, 3 insertions(+), 2 deletions(-) diff --git a/src/raiju/raijuColdStartHelper.F90 b/src/raiju/raijuColdStartHelper.F90 index a57aa450..ad431cd7 100644 --- a/src/raiju/raijuColdStartHelper.F90 +++ b/src/raiju/raijuColdStartHelper.F90 @@ -147,8 +147,7 @@ module raijuColdStartHelper write(*,'(a,f7.2)') " post-loss : ",dps_postCX write(*,'(a,f7.2)') " post-rescale : ",dps_rescale write(*,'(a,f7.2)') " Hot electron DPS-Dst : ",dps_ele - el7 -2 write(*,'(a,f7.2)') " Real Dst : ",dstReal + write(*,'(a,f7.2)') " Real Dst : ",dstReal write(*,'(a,f7.2)') " Model Dst : ",dstModel write(*,'(a,f7.2)') " Current DPS-Dst : ",dps_current write(*,'(a,f7.2)') " Target DPS-Dst : ",dstTarget diff --git a/src/voltron/voltapp.F90 b/src/voltron/voltapp.F90 index 7232ea23..bbb3b0ab 100644 --- a/src/voltron/voltapp.F90 +++ b/src/voltron/voltapp.F90 @@ -99,6 +99,8 @@ module voltapp ! adjust XMl reader root call xmlInp%SetRootStr('Kaiju/Voltron') + ! Make sure verbosity is still right after others do stuff with the reader + call xmlInp%SetVerbose(vApp%isLoud) !Initialize planet information call getPlanetParams(vApp%planet, xmlInp) From 6d8c916213bcec37bc6b88ccafb2b771807a681a Mon Sep 17 00:00:00 2001 From: Anthony Date: Tue, 16 Sep 2025 12:14:57 -0600 Subject: [PATCH 10/15] updating volt shellGrid unit test --- tests/voltron/testVoltGridGen.pf | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/tests/voltron/testVoltGridGen.pf b/tests/voltron/testVoltGridGen.pf index cc079811..b0fe5105 100644 --- a/tests/voltron/testVoltGridGen.pf +++ b/tests/voltron/testVoltGridGen.pf @@ -9,6 +9,8 @@ module testVoltGridGen implicit none character(len=strLen) :: xmlName = 'voltGridTests.xml' + integer, parameter :: gamRes = 64 + !! genVoltShellGrid expects a gamera resolution to determine defaults for its own contains @@ -36,7 +38,7 @@ module testVoltGridGen call xmlInp%Set_Val(Nt, "grid/Nt", -1) ! -1 so things blow up if xml isn't set properly call xmlInp%Set_Val(Np, "grid/Np", -1) ! -1 so things blow up if xml isn't set properly - call genVoltShellGrid(vApp, xmlInp) + call genVoltShellGrid(vApp, xmlInp, gamRes) @assertEqual(vApp%shGrid%Nt, 2*Nt, "Wrong amount of theta cells") @assertEqual(vApp%shGrid%Np, Np , "Wrong amount of phi cells") @@ -62,7 +64,7 @@ module testVoltGridGen call xmlInp%Set_Val(Nt, "grid/Nt", -1) ! -1 so things blow up if xml isn't set properly call xmlInp%Set_Val(Np, "grid/Np", -1) ! -1 so things blow up if xml isn't set properly - call genVoltShellGrid(vApp, xmlInp) + call genVoltShellGrid(vApp, xmlInp, gamRes) @assertEqual(2*Nt, vApp%shGrid%Nt, "Wrong amount of theta cells") @assertEqual( Np, vApp%shGrid%Np, "Wrong amount of phi cells") From 85f419350bedeb966b319bfe4c80ce20cf6a2003 Mon Sep 17 00:00:00 2001 From: "Anthony M. Sciola" Date: Fri, 19 Sep 2025 09:32:50 -0700 Subject: [PATCH 11/15] (Not tested) Changing coldstart to only add RC if needed --- src/raiju/raijuColdStartHelper.F90 | 92 ++++++++++-------------------- src/raiju/raijuPreAdvancer.F90 | 3 + 2 files changed, 34 insertions(+), 61 deletions(-) diff --git a/src/raiju/raijuColdStartHelper.F90 b/src/raiju/raijuColdStartHelper.F90 index ad431cd7..8eb0d368 100644 --- a/src/raiju/raijuColdStartHelper.F90 +++ b/src/raiju/raijuColdStartHelper.F90 @@ -88,48 +88,38 @@ module raijuColdStartHelper ! Update Dst target dstReal = GetSWVal('symh', Model%tsF, t0) - if (.not. cs%doneFirstCS) then + if (cs%doneFirstCS) then + write(*,*)"Already coldstarted once, you shouldn't be here " + return + else ! On first try, we assume there is no existing ring current, and its our job to make up the entire difference dstTarget = dstReal - dstModel - doInitRC = .true. - else if (t0 > (cs%lastEval + cs%evalCadence)) then - ! If we are updating, there should already be some ring current - ! If dstReal - dstModel is still < 0, we need to add ADDITIONAL pressure to get them to match - dps_current = spcEta2DPS(Model, Grid, State%bvol_cc, State%eta, Grid%spc(sIdx_p), isGood) + spcEta2DPS(Model, Grid, State%bvol_cc, State%eta, Grid%spc(sIdx_e), isGood) - dstTarget = dstReal - (dstModel - dps_current) - doInitRC = .false. endif cs%lastEval = t0 cs%lastTarget = dstTarget cs%doneFirstCS = .true. ! Whether we do anything or not, we were at least called once - - ! Init psphere (will always run if this function is called) - call setRaijuInitPsphere(Model, Grid, State, Model%psphInitKp) ! Now decide if we need to add a starter ring current - !if (dstTarget >= 0) then ! We've got nothing to contribute - ! write(*,*)"RAIJU coldstart not adding anything" - ! write(*,*)'doAccumulate=',doAccumulate - ! return - !endif - - if (doInitRC) then - ! Init hot protons - call raiColdStart_initHOTP(Model, Grid, State, t0, dstTarget, etaCS) - !call raiColdStart_initHOTP_RCOnly(Model, Grid, State, t0, dstTarget, etaCS) - dps_preCX = spcEta2DPS(Model, Grid, State%bvol_cc, etaCS, Grid%spc(sIdx_p), isGood) - ! Hit it with some charge exchange - if (cs%doCX) then - call raiColdStart_applyCX(Model, Grid, State, Grid%spc(sIdx_p), etaCS) - endif - dps_postCX = spcEta2DPS(Model, Grid, State%bvol_cc, etaCS, Grid%spc(sIdx_p), isGood) - ! Use HOTP moments to set electrons - call raiColdStart_initHOTE(Model, Grid, State, etaCS) - dps_ele = spcEta2DPS(Model, Grid, State%bvol_cc, etaCS, Grid%spc(sIdx_e), isGood) - dps_current = dps_postCX ! Note: if using fudge we're gonna lose electrons immediately, don't include them in current dst for now + if (dstTarget >= 0) then ! We've got nothing to contribute + write(*,*)"RAIJU coldstart not adding starter ring current" + return endif + ! Init hot protons + !call raiColdStart_initHOTP(Model, Grid, State, t0, dstTarget, etaCS) + call raiColdStart_initHOTP_RCOnly(Model, Grid, State, t0, dstTarget, etaCS) + dps_preCX = spcEta2DPS(Model, Grid, State%bvol_cc, etaCS, Grid%spc(sIdx_p), isGood) + ! Hit it with some charge exchange + if (cs%doCX) then + call raiColdStart_applyCX(Model, Grid, State, Grid%spc(sIdx_p), etaCS) + endif + dps_postCX = spcEta2DPS(Model, Grid, State%bvol_cc, etaCS, Grid%spc(sIdx_p), isGood) + ! Use HOTP moments to set electrons + call raiColdStart_initHOTE(Model, Grid, State, etaCS) + dps_ele = spcEta2DPS(Model, Grid, State%bvol_cc, etaCS, Grid%spc(sIdx_e), isGood) + dps_current = dps_postCX ! Note: if using fudge we're gonna lose electrons immediately, don't include them in current dst for now + if (dstTarget < 0) then etaScale = abs(dstTarget / dps_current) etaCS(:,:,Grid%spc(sIdx_p)%kStart:Grid%spc(sIdx_p)%kEnd) = etaScale*etaCS(:,:,Grid%spc(sIdx_p)%kStart:Grid%spc(sIdx_p)%kEnd) @@ -138,41 +128,20 @@ module raijuColdStartHelper dps_rescale = dps_current endif - if (doInitRC) then - write(*,*) "RAIJU Cold starting..." - write(*,'(a,f7.2)') " Real Dst : ",dstReal - write(*,'(a,f7.2)') " Model Dst : ",dstModel - write(*,'(a,f7.2)') " Target DPS-Dst : ",dstTarget - write(*,'(a,f7.2)') " Hot proton pre-loss : ",dps_preCX - write(*,'(a,f7.2)') " post-loss : ",dps_postCX - write(*,'(a,f7.2)') " post-rescale : ",dps_rescale - write(*,'(a,f7.2)') " Hot electron DPS-Dst : ",dps_ele - write(*,'(a,f7.2)') " Real Dst : ",dstReal - write(*,'(a,f7.2)') " Model Dst : ",dstModel - write(*,'(a,f7.2)') " Current DPS-Dst : ",dps_current - write(*,'(a,f7.2)') " Target DPS-Dst : ",dstTarget - write(*,'(a,f7.2)') " post-rescale : ",dps_rescale - write(*,'(a,f7.2)') " Hot electron DPS-Dst : ",dps_ele - endif + write(*,*) "RAIJU Cold starting..." + write(*,'(a,f7.2)') " Real Dst : ",dstReal + write(*,'(a,f7.2)') " Model Dst : ",dstModel + write(*,'(a,f7.2)') " Target DPS-Dst : ",dstTarget + write(*,'(a,f7.2)') " Hot proton pre-loss : ",dps_preCX + write(*,'(a,f7.2)') " post-loss : ",dps_postCX + write(*,'(a,f7.2)') " post-rescale : ",dps_rescale + write(*,'(a,f7.2)') " Hot electron DPS-Dst : ",dps_ele end associate ! finally, put it into raiju state if(doAccumulate) then - !State%eta = State%eta + etaCS - associate(sh=>Grid%shGrid) - !$OMP PARALLEL DO default(shared) & - !$OMP private(i,j,k) - do k=1,Grid%Nk - do j=sh%jsg,sh%jeg - do i=sh%isg,sh%ieg - if (etaCS(i,j,k) > State%eta(i,j,k)) then - State%eta(i,j,k) = etaCS(i,j,k) - endif - enddo - enddo - enddo - end associate + State%eta = State%eta + etaCS else State%eta = etaCS endif @@ -245,6 +214,7 @@ module raijuColdStartHelper end subroutine raiColdStart_initHOTP_RCOnly + subroutine raiColdStart_initHOTP(Model, Grid, State, t0, dstTarget, etaCS) type(raijuModel_T), intent(in) :: Model type(raijuGrid_T), intent(in) :: Grid diff --git a/src/raiju/raijuPreAdvancer.F90 b/src/raiju/raijuPreAdvancer.F90 index 36819d69..582450b6 100644 --- a/src/raiju/raijuPreAdvancer.F90 +++ b/src/raiju/raijuPreAdvancer.F90 @@ -41,6 +41,9 @@ module raijuPreAdvancer ! Moments to etas, initial active shell calculation call Tic("BCs") + if (State%isFirstCpl) then + call setRaijuInitPsphere(Model, Grid, State, Model%psphInitKp) + endif call applyRaijuBCs(Model, Grid, State, doWholeDomainO=State%isFirstCpl) ! If fullEtaMap=True, mom2eta map is applied to the whole domain call Toc("BCs") From d1cae5259b6a4d8ebe3b2935b36c7f3d9b91813c Mon Sep 17 00:00:00 2001 From: Anthony Date: Wed, 24 Sep 2025 15:06:04 -0600 Subject: [PATCH 12/15] Setting some default dipole grads for points on the corners --- src/raiju/raijuPreAdvancer.F90 | 12 +++++++++--- 1 file changed, 9 insertions(+), 3 deletions(-) diff --git a/src/raiju/raijuPreAdvancer.F90 b/src/raiju/raijuPreAdvancer.F90 index 36819d69..d5578852 100644 --- a/src/raiju/raijuPreAdvancer.F90 +++ b/src/raiju/raijuPreAdvancer.F90 @@ -487,9 +487,15 @@ module raijuPreAdvancer !$OMP private(i,j,bVolcc) do j=sh%jsg,sh%jeg do i=sh%isg,sh%ieg - !bVolcc = toCenter2D(dV(i:i+1,j:j+1)) + DipFTV_colat(Grid%thcRp(i), B0) ! Will include smoothing of dV if enabled - bVolcc = toCenter2D(V(i:i+1,j:j+1)) - gradVM(i,j,:) = (-2./3.)*bVolcc**(-5./3.)*gradVM(i,j,:) + if(all(isGcorner(i:i+1,j:j+1))) then + !bVolcc = toCenter2D(dV(i:i+1,j:j+1)) + DipFTV_colat(Grid%thcRp(i), B0) ! Will include smoothing of dV if enabled + bVolcc = toCenter2D(V(i:i+1,j:j+1)) + gradVM(i,j,:) = (-2./3.)*bVolcc**(-5./3.)*gradVM(i,j,:) + else + ! gradVM should be zero for this point coming out of calcGradIJ_cc, but set to dipole value just in case + gradVM(i,j,RAI_PH) = 0.0 + gradVM(i,j,RAI_TH) = (-2./3.)*DipFTV_colat(Grid%thcRp(i), B0)**(-5./3.)*dV0_dth(i,j) + endif enddo enddo From 2b16b820b3823cf0c14daf57ba689ad6ecfcf4fe Mon Sep 17 00:00:00 2001 From: Anthony Date: Thu, 25 Sep 2025 13:56:38 -0600 Subject: [PATCH 13/15] Putting the plasma back in plasmasphere --- src/raiju/raijuPreAdvancer.F90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/raiju/raijuPreAdvancer.F90 b/src/raiju/raijuPreAdvancer.F90 index 45e8898f..558dc6a1 100644 --- a/src/raiju/raijuPreAdvancer.F90 +++ b/src/raiju/raijuPreAdvancer.F90 @@ -41,10 +41,10 @@ module raijuPreAdvancer ! Moments to etas, initial active shell calculation call Tic("BCs") + call applyRaijuBCs(Model, Grid, State, doWholeDomainO=State%isFirstCpl) ! If fullEtaMap=True, mom2eta map is applied to the whole domain if (State%isFirstCpl) then call setRaijuInitPsphere(Model, Grid, State, Model%psphInitKp) endif - call applyRaijuBCs(Model, Grid, State, doWholeDomainO=State%isFirstCpl) ! If fullEtaMap=True, mom2eta map is applied to the whole domain call Toc("BCs") ! Handle plasmasphere refilling for the full step about to happen From 450c953455d52fcae83ac791f052049a4ef80456 Mon Sep 17 00:00:00 2001 From: Anthony Date: Fri, 26 Sep 2025 13:32:51 -0600 Subject: [PATCH 14/15] Reverting coldstart method to setting full domain. Zero gradVM for non-good points. NaN catch for EAVG calculation in CoupleIMaxToMix --- src/raiju/raijuColdStartHelper.F90 | 16 +++++++++++++--- src/raiju/raijuPreAdvancer.F90 | 3 ++- .../modelInterfaces/imag2mix_interface.F90 | 2 +- 3 files changed, 16 insertions(+), 5 deletions(-) diff --git a/src/raiju/raijuColdStartHelper.F90 b/src/raiju/raijuColdStartHelper.F90 index 8eb0d368..f318666c 100644 --- a/src/raiju/raijuColdStartHelper.F90 +++ b/src/raiju/raijuColdStartHelper.F90 @@ -107,8 +107,8 @@ module raijuColdStartHelper endif ! Init hot protons - !call raiColdStart_initHOTP(Model, Grid, State, t0, dstTarget, etaCS) - call raiColdStart_initHOTP_RCOnly(Model, Grid, State, t0, dstTarget, etaCS) + call raiColdStart_initHOTP(Model, Grid, State, t0, dstTarget, etaCS) + !call raiColdStart_initHOTP_RCOnly(Model, Grid, State, t0, dstTarget, etaCS) dps_preCX = spcEta2DPS(Model, Grid, State%bvol_cc, etaCS, Grid%spc(sIdx_p), isGood) ! Hit it with some charge exchange if (cs%doCX) then @@ -141,7 +141,17 @@ module raijuColdStartHelper ! finally, put it into raiju state if(doAccumulate) then - State%eta = State%eta + etaCS + !$OMP PARALLEL DO default(shared) & + !$OMP private(i,j,k) + do j=Grid%shGrid%jsg,Grid%shGrid%jeg + do i=Grid%shGrid%isg,Grid%shGrid%ieg + do k=1,Grid%Nk + if(etaCS(i,j,k) > State%eta(i,j,k)) then + State%eta(i,j,k) = etaCS(i,j,k) + endif + enddo + enddo + enddo else State%eta = etaCS endif diff --git a/src/raiju/raijuPreAdvancer.F90 b/src/raiju/raijuPreAdvancer.F90 index 558dc6a1..f87cb897 100644 --- a/src/raiju/raijuPreAdvancer.F90 +++ b/src/raiju/raijuPreAdvancer.F90 @@ -497,7 +497,8 @@ module raijuPreAdvancer else ! gradVM should be zero for this point coming out of calcGradIJ_cc, but set to dipole value just in case gradVM(i,j,RAI_PH) = 0.0 - gradVM(i,j,RAI_TH) = (-2./3.)*DipFTV_colat(Grid%thcRp(i), B0)**(-5./3.)*dV0_dth(i,j) + gradVM(i,j,RAI_TH) = 0.0 + !gradVM(i,j,RAI_TH) = (-2./3.)*DipFTV_colat(Grid%thcRp(i), B0)**(-5./3.)*dV0_dth(i,j) endif enddo enddo diff --git a/src/voltron/modelInterfaces/imag2mix_interface.F90 b/src/voltron/modelInterfaces/imag2mix_interface.F90 index d54f2363..d6976fe5 100644 --- a/src/voltron/modelInterfaces/imag2mix_interface.F90 +++ b/src/voltron/modelInterfaces/imag2mix_interface.F90 @@ -231,7 +231,7 @@ module imag2mix_interface imP_avg(RAI_EDEN) = imP_avg(RAI_EDEN )/nPnts**2 imP_avg(RAI_EPRE) = imP_avg(RAI_EPRE )/nPnts**2 imP_avg(RAI_NPSP) = imP_avg(RAI_NPSP )/nPnts**2 - imP_avg(RAI_EAVG) = imP_avg(RAI_EFLUX) / imP_avg(RAI_ENFLX) + imP_avg(RAI_EAVG) = imP_avg(RAI_EFLUX) / max(TINY, imP_avg(RAI_ENFLX)) imP_avg(RAI_GTYPE) = imP_avg(RAI_GTYPE)/nGood imP_avg(RAI_THCON) = imP_avg(RAI_THCON)/nGood imP_avg(RAI_PHCON) = imP_avg(RAI_PHCON)/nGood From 3c92508ed61b7e105d636099bb0da090a92e6c66 Mon Sep 17 00:00:00 2001 From: Kareem Sorathia Date: Tue, 30 Sep 2025 15:37:42 -0400 Subject: [PATCH 15/15] Applying band-aid to fix dblfix --- .../modelInterfaces/raiCplTypesSub.F90 | 29 ++++++++++--------- 1 file changed, 15 insertions(+), 14 deletions(-) diff --git a/src/voltron/modelInterfaces/raiCplTypesSub.F90 b/src/voltron/modelInterfaces/raiCplTypesSub.F90 index 03a4d88e..01d3bd32 100644 --- a/src/voltron/modelInterfaces/raiCplTypesSub.F90 +++ b/src/voltron/modelInterfaces/raiCplTypesSub.F90 @@ -263,21 +263,22 @@ submodule (volttypes) raijuCplTypesSub endif enddo - ! Mock up cold electrons balancing hot protons and see if it produces meaningful flux - call InterpShellVar_TSC_pnt(sh, State%Den_avg(idx_proton) , th, ph, d_ion) - call InterpShellVar_TSC_pnt(sh, State%Press_avg(idx_proton) , th, ph, p_ion) - pie_frac = 0.05 ! Fraction of ion pressure contained by these neutralizing electrons - pe_kT = DP2kT(d_ion, p_ion*pie_frac) ! [keV] - pe_nflux = imP(RAI_ENFLX)*d_ion/d_hot ! Scale number flux to same loss processes except there were d_ion density instead of d_electron density - pe_eflux = (pe_kT*kev2erg)*pe_nflux ! [erg/cm^2/s] - if (pe_eflux > imP(RAI_EFLUX)) then ! Use in place of normal flux only if energy flux for these are greater than hot electron channel fluxes - imP(RAI_EFLUX) = pe_eflux - imP(RAI_ENFLX) = pe_nflux - imP(RAI_EDEN ) = d_ion*1.0e6 ! [#/m^3] - imP(RAI_EPRE ) = p_ion*pie_frac*1.0e-9 ! [Pa] - endif - endif + !K: Removing this code for now, should be rewritten to use the MHD D,P => precip routines + ! call InterpShellVar_TSC_pnt(sh, State%Den_avg(idx_proton) , th, ph, d_ion) + ! call InterpShellVar_TSC_pnt(sh, State%Press_avg(idx_proton) , th, ph, p_ion) + ! pie_frac = 0.05 ! Fraction of ion pressure contained by these neutralizing electrons + ! pe_kT = DP2kT(d_ion, p_ion*pie_frac) ! [keV] + ! pe_nflux = imP(RAI_ENFLX)*d_ion/d_hot ! Scale number flux to same loss processes except there were d_ion density instead of d_electron density + ! pe_eflux = (pe_kT*kev2erg)*pe_nflux ! [erg/cm^2/s] + ! if (pe_eflux > imP(RAI_EFLUX)) then ! Use in place of normal flux only if energy flux for these are greater than hot electron channel fluxes + ! imP(RAI_EFLUX) = pe_eflux + ! imP(RAI_ENFLX) = pe_nflux + ! imP(RAI_EDEN ) = d_ion*1.0e6 ! [#/m^3] + ! imP(RAI_EPRE ) = p_ion*pie_frac*1.0e-9 ! [Pa] + ! endif + + endif !spcList(s)%spcType == X enddo ! derive mean energy where nflux is non-trivial. if (imP(RAI_ENFLX) > TINY) imP(RAI_EAVG) = imP(RAI_EFLUX)/imP(RAI_ENFLX) * erg2kev ! Avg E [keV]