Blog: Update 4/18: charge_exchange_updates

File charge_exchange_updates, 8.8 KB (added by adebrech, 8 years ago)
Line 
1diff --git a/src/physics/EOS.f90 b/src/physics/EOS.f90
2index d55d914..aaad039 100644
3--- a/src/physics/EOS.f90
4+++ b/src/physics/EOS.f90
5@@ -180,11 +180,11 @@ CONTAINS
6 gamma15 = gamma/gamma1
7
8
9- IF (iEOS /= EOS_MULTISPECIES_GAS .AND. (lTrackHydrogen .OR. lTrackHelium)) THEN
10+ IF (iEOS /= EOS_MULTISPECIES_GAS .AND. (lTrackHydrogen .OR. lTrackHelium .or. lChargeExchange)) THEN
11 IF (MPI_ID == 0) THEN
12 PRINT*, "iEOS should = 1 (MultiSpecies gas) for tracking hydrogen and helium."
13 PRINT*, "Setting iEOS to 1. If you wish to run with iEOS /= 1, then stop this"
14- PRINT*, "run, set lTrackHydrogen and lTrackHelium to false, and run again."
15+ PRINT*, "run, set lTrackHydrogen, lTrackHelium, and lChargeExchange to false, and run again."
16 END IF
17 iEOS = EOS_MULTISPECIES_GAS
18 END IF
19@@ -226,6 +226,13 @@ CONTAINS
20 CALL AddSpecies(iHII, 'HII', 1.0001d0, muHII, 1)
21
22 END IF
23+
24+ if (lChargeExchange) then
25+ call AddSpecies(iHhot, 'hot H', 1.0001d0, muH, 0)
26+ call AddSpecies(iHcold, 'cold H', 1.0001d0, muH, 0)
27+ call AddSpecies(iHIIhot, 'hot HII', 1.0001d0, muHII, 1)
28+ call AddSpecies(iHIIcold, 'cold HII', 1.0001d0, muHII, 1)
29+ end if
30
31 IF (lTrackHelium) THEN
32 CALL AddSpecies(iHe, 'He', 5d0/3d0, muHe, 0)
33diff --git a/src/physics/abundances.f90 b/src/physics/abundances.f90
34index 3903cf0..f85afbb 100644
35--- a/src/physics/abundances.f90
36+++ b/src/physics/abundances.f90
37@@ -89,7 +89,7 @@ CONTAINS
38 untracked_n=Z_n !ZFrac*(1d0+xZ)
39 gamma7_ut=Z_n/(z_gamma-1d0)
40
41- IF (lTrackHydrogen) THEN
42+ IF (lTrackHydrogen .or. lChargeExchange) THEN
43 tracked_mass=tracked_mass+H_m !muH
44 tracked_n=tracked_n+H_n
45 ELSE
46diff --git a/src/physics/physics_control.f90 b/src/physics/physics_control.f90
47index 0c1ae36..1e5587f 100644
48--- a/src/physics/physics_control.f90
49+++ b/src/physics/physics_control.f90
50@@ -59,7 +59,7 @@ CONTAINS
51
52 NAMELIST /PhysicsData/ &
53 InterpOpts, refineVariableFactor, &
54- iSelfGravity, DefaultAccretionRoutine, UniformGravity, OmegaRot, iCylindrical, iCooling, FloorTemp,SrcPrecision, lTrackHydrogen, lTrackHelium, lTrackSulfur, SolarIonizingFlux, &
55+ iSelfGravity, DefaultAccretionRoutine, UniformGravity, OmegaRot, iCylindrical, iCooling, FloorTemp,SrcPrecision, lTrackHydrogen, lChargeExchange, lTrackHelium, lTrackSulfur, SolarIonizingFlux, &
56 lRestartOnDensityProtections, lTrackDensityProtections, iDensityProtect, iMomentumProtect, MinDensity,&
57 lRestartOnPressureProtections, lTrackPressureProtections, iPressureProtect, MinTemp, lTrackMaxSpeeds, &
58 lRestartOnHypreMaxIters, &
59@@ -113,6 +113,11 @@ CONTAINS
60 END IF!icyl
61 R2DEff=(sqrt(product(GxBounds(1:2,2)-GxBounds(1:2,1)))/sqrt(pi)*exp(-.5d0))
62 ! Source Terms are not implemented in this iteration
63+
64+ if (lChargeExchange .and. lTrackHydrogen) then
65+ if (MPI_ID .eq. 0) print *, 'Hydrogen totals can be computed from hot and cold effects. Setting lTrackHydrogen to false.'
66+ lTrackHydrogen = .false.
67+ end if
68
69 CALL SetConservedVariableIndices()
70 CALL AbundanceInit()
71diff --git a/src/physics/physics_declarations.f90 b/src/physics/physics_declarations.f90
72index ebb7ff4..92c1d35 100644
73--- a/src/physics/physics_declarations.f90
74+++ b/src/physics/physics_declarations.f90
75@@ -111,7 +111,7 @@ MODULE PhysicsDeclarations
76 INTEGER :: irho = 0, iE = 0,ivx = 0, ivy = 0, ivz = 0, iBx = 0, iBy = 0, iBz = 0, iPhiDot=0, iPhiGas=0, iEntropy=0, iAngMom=0,imom(3),ivOld(3),irhoOld, iB(3), iE_Rad=0, iE_RadDot=0, iKappa, iKappaDot
77 INTEGER :: iTemp=0, iTempDot=0
78
79- INTEGER, PUBLIC :: iH2=0,iH=0,iHII=0,iHe=0,iHeII=0,iHeIII=0,iSII=0,iSIII=0,iSIV=0
80+ INTEGER, PUBLIC :: iH2=0,iH=0,iHII=0,iHe=0,iHeII=0,iHeIII=0,iSII=0,iSIII=0,iSIV=0,iHhot=0,iHcold=0,iHIIhot=0,iHIIcold=0
81 INTEGER :: iRiemannMaxSpeed=0,iEigenMaxSpeed=0,iFluidMaxSpeed=0
82 INTEGER :: iDensityProtections=0,iPressureProtections=0
83
84@@ -152,6 +152,7 @@ MODULE PhysicsDeclarations
85
86 INTEGER :: iCylindrical, iCooling !, iSymAxis ! Source switches
87 LOGICAL :: lTrackHydrogen=.false.
88+ logical :: lChargeExchange=.false.
89 LOGICAL :: lTrackHelium=.false.
90 LOGICAL :: lTrackSulfur=.false.
91 REAL(KIND=qPREC) :: ne_min=0.01 !minimum ratio of electrons to total particles
92diff --git a/src/source/chemistry.f90 b/src/source/chemistry.f90
93index af6c0c2..85d127e 100644
94--- a/src/source/chemistry.f90
95+++ b/src/source/chemistry.f90
96@@ -20,10 +20,10 @@ CONTAINS
97
98
99 SUBROUTINE ChemistryInit
100- lChemistry = ANY((/lTrackHydrogen, lTrackHelium, lTrackSulfur/))
101+ lChemistry = ANY((/lTrackHydrogen, lChargeExchange, lTrackHelium, lTrackSulfur/))
102 IF (lChemistry) THEN
103 iErr=0
104- IF (lTrackHydrogen) THEN
105+ IF (lTrackHydrogen .or. lChargeExchange) THEN
106 CALL ReadTable('TABLES/H_Ioniz.tab', H_IonizTab, iErr)
107 CALL ReadTable('TABLES/HII_Recomb.tab', HII_RecombTab, iErr)
108 END IF
109@@ -52,7 +52,7 @@ CONTAINS
110 REAL(KIND=qPrec), DIMENSION(1:NrHydroVars) :: q
111 REAL(KIND=qPrec), DIMENSION(1:NrHydroVars) :: dqdt
112 REAL(KIND=qPrec), OPTIONAL :: Temp
113- REAL(KIND=qPREC) :: Hion, HIIrec, Heion, HeIIrec, HeIIion, HeIIIrec, ne, dn, logT, &
114+ REAL(KIND=qPREC) :: Hion, HIIrec, Heion, HeIIrec, HeIIion, HeIIIrec, ne, dn, dnCold, dnHot, logT, &
115 SIIion, SIIIion, SIIIrec, SIVrec
116
117 logT=log10(Temp)
118@@ -65,9 +65,27 @@ CONTAINS
119 dqdt(iH) = dqdt(iH) - dn
120 dqdt(iHII) = dqdt(iHII) + dn
121 dqdt(iE) = dqdt(iE) - (IonH*Hion + (Boltzmann*Temp)*HIIrec)*eDotScale
122-
123-
124 END IF
125+
126+ if (lChargeExchange) then
127+ Hion = 10d0**GetTableValue(H_IonizTab, logT)*ne*(q(iHhot) + q(iHcold))*nScale**2d0
128+ HIIrec = 10d0**GetTableValue(HII_RecombTab, logT)*ne*(q(iHIIhot) + q(iHIIcold))*nScale**2d0
129+ dnCold = (Hion*q(iHcold)/(q(iHhot) + q(iHcold))-HIIrec*q(iHIIcold)/(q(iHIIhot) + q(iHIIcold)))*ChemScale
130+ dnHot = (Hion*q(iHhot)/(q(iHhot) + q(iHcold))-HIIrec*q(iHIIhot)/(q(iHIIhot) + q(iHIIcold)))*ChemScale
131+
132+ dqdt(iHhot) = dqdt(iHhot) - dnHot
133+ dqdt(iHcold) = dqdt(iHcold) - dnCold
134+ dqdt(iHIIhot) = dqdt(iHIIhot) + dnHot
135+ dqdt(iHIIcold) = dqdt(iHIIcold) + dnCold
136+
137+ dqdt(iE) = dqdt(iE) - (IonH*Hion + (Boltzmann*Temp)*HIIrec)*eDotScale
138+
139+ netExchange = chargeExchangeRate*(q(iHIIhot)*q(iHcold) - q(iHhot)*q(iHIIcold))
140+ dqdt(iHhot) = dqdt(iHhot) + netExchange
141+ dqdt(iHIIhot) = dqdt(iHIIhot) - netExchange
142+ dqdt(iHcold) = dqdt(iHcold) - netExchange
143+ dqdt(iHIIcold) = dqdt(iHIIcold) + netExchange
144+ end if
145
146 IF (lTrackHelium) THEN
147 Heion = 10d0**GetTableValue(He_IonizTab, logT)*ne*q(iHe)*nScale**2d0
148diff --git a/src/source/source_control.f90 b/src/source/source_control.f90
149index 719dbfe..86ca0c2 100644
150--- a/src/source/source_control.f90
151+++ b/src/source/source_control.f90
152@@ -431,24 +431,34 @@ CONTAINS
153 TYPE(InfoDef) :: Info
154 INTEGER :: ip(3)
155 REAL(KIND=qPrec) :: q(:), ne, Temp
156- REAL(KIND=qPrec) :: pos(3), SrcDerivs(NrHydroVars), t, r, dn
157+ REAL(KIND=qPrec) :: pos(3), SrcDerivs(NrHydroVars), t, r, dn, dnCold, dnHot
158 ! Internal declarations
159 REAL(KIND=qPrec) :: dqdt(NrHydroVars)
160 dqdt=0d0
161
162 IF (iCooling /= 0 .OR. lChemistry) THEN
163 Temp=SourceTemperature(q)
164- IF (lTrackHydrogen .OR. lTrackHelium) ne=get_ne(q)
165+ IF (lTrackHydrogen .or. lChargeExchange .OR. lTrackHelium) ne=get_ne(q)
166 IF (lChemistry) CALL ChemReactions(q, dqdt, ne, Temp)
167+ !!! Would be good to generalize this eventually
168 IF (SolarIonizingFlux > 0) THEN
169 r=sqrt(sum(pos**2))
170- dn=q(iH)*nScale*SolarIonizingFlux/(r*lscale)**2
171- dqdt(iH)=dqdt(iH)-dn*ChemScale
172- dqdt(iHII)=dqdt(iHII)+dn*ChemScale
173- if (r < .4) THEN
174- dqdt(iH)=0d0
175- dqdt(iHII)=0d0
176- END if
177+ if (lChargeExchange) then
178+ dnHot=q(iHhot)*nScale*SolarIonizingFlux/(r*lscale)**2
179+ dnCold=q(iHcold)*nScale*SolarIonizingFlux/(r*lscale)**2
180+ dqdt(iHhot)=dqdt(iHhot)-dnHot*ChemScale
181+ dqdt(iHIIhot)=dqdt(iHIIhot)+dnHot*ChemScale
182+ dqdt(iHcold)=dqdt(iHcold)-dnCold*ChemScale
183+ dqdt(iHIIcold)=dqdt(iHIIcold)+dnCold*ChemScale
184+ else
185+ dn=q(iH)*nScale*SolarIonizingFlux/(r*lscale)**2
186+ dqdt(iH)=dqdt(iH)-dn*ChemScale
187+ dqdt(iHII)=dqdt(iHII)+dn*ChemScale
188+ if (r < .325) THEN
189+ dqdt(iH)=0d0
190+ dqdt(iHII)=0d0
191+ END if
192+ end if
193 END IF
194 IF (Temp > FloorTemp) THEN
195 IF (iCooling /= 0) CALL Cooling(q, dqdt, ne, Temp)