Skip to content

Commit 779d3a8

Browse files
authored
Merge pull request #1 from chongchonghe/chong/photo-fix
fix: correct energy normalization and electron mass in photoionizatio…
2 parents e177e7a + 31925d8 commit 779d3a8

1 file changed

Lines changed: 13 additions & 5 deletions

File tree

EOS/photoionization/actual_eos.H

Lines changed: 13 additions & 5 deletions
Original file line numberDiff line numberDiff line change
@@ -41,6 +41,12 @@ void actual_eos_init ()
4141
GET_SPECIES_PARAMS(3);
4242

4343
#undef GET_SPECIES_PARAMS
44+
45+
// Electrons have actual electron mass, not proton mass
46+
int const e_idx = network_spec_index("Electron");
47+
if (e_idx >= 0) {
48+
spmasses[e_idx] = C::m_e;
49+
}
4450
}
4551

4652

@@ -77,7 +83,9 @@ void actual_eos (I input, T& state)
7783
sum_ni_fi_over_2 += state.xn[n] / (gammas[n] - 1.0);
7884
sum_Abarinv += state.xn[n];
7985
}
80-
sum_ni_fi_over_2 /= sum_Abarinv;
86+
// fi_over_2_avg is the per-particle average heat capacity (used for cs, G, pressure)
87+
// sum_ni_fi_over_2 remains the total (unnormalized) heat capacity (used for e, T, dedT)
88+
const amrex::Real fi_over_2_avg = sum_ni_fi_over_2 / sum_Abarinv;
8189
sum_Abarinv *= protonmass/rhotot;
8290
state.mu = 1.0 / sum_Abarinv;
8391

@@ -139,7 +147,7 @@ void actual_eos (I input, T& state)
139147

140148
// stop the integration if the pressure < 0
141149
AMREX_ASSERT(state.p > 0.);
142-
eint = state.p * sum_ni_fi_over_2 / dens;
150+
eint = state.p * fi_over_2_avg / dens;
143151
temp = eint * rhotot / (sum_ni_fi_over_2 * C::k_B);
144152
}
145153
break;
@@ -191,14 +199,14 @@ void actual_eos (I input, T& state)
191199
}
192200

193201
if constexpr (has_pressure<T>::value) {
194-
amrex::Real pressure = state.rho * eint / sum_ni_fi_over_2;
202+
amrex::Real pressure = state.rho * eint / fi_over_2_avg;
195203
state.p = pressure;
196204

197205
state.dpdT = pressure / temp;
198206
state.dpdr = pressure / dens;
199-
state.cs = std::sqrt((1.0 + 1.0/sum_ni_fi_over_2) * state.p /state.rho);
207+
state.cs = std::sqrt((1.0 + 1.0/fi_over_2_avg) * state.p /state.rho);
200208
if constexpr (has_G<T>::value) {
201-
state.G = 0.5 * (1.0 + (1.0 + 1.0/sum_ni_fi_over_2));
209+
state.G = 0.5 * (1.0 + (1.0 + 1.0/fi_over_2_avg));
202210
}
203211
}
204212

0 commit comments

Comments
 (0)