Line data Source code
1 : ! ***********************************************************************
2 : !
3 : ! Copyright (C) 2015-2022 Evan Bauer & The MESA Team
4 : !
5 : ! This program is free software: you can redistribute it and/or modify
6 : ! it under the terms of the GNU Lesser General Public License
7 : ! as published by the Free Software Foundation,
8 : ! either version 3 of the License, or (at your option) any later version.
9 : !
10 : ! This program is distributed in the hope that it will be useful,
11 : ! but WITHOUT ANY WARRANTY; without even the implied warranty of
12 : ! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
13 : ! See the GNU Lesser General Public License for more details.
14 : !
15 : ! You should have received a copy of the GNU Lesser General Public License
16 : ! along with this program. If not, see <https://www.gnu.org/licenses/>.
17 : !
18 : ! ***********************************************************************
19 :
20 : module diffusion_support
21 :
22 : use const_def, only: dp, pi, pi4, me, mp, qe, amu, fine, boltzm, hbar, msun, rsun, secyer, one_third, four_thirds_pi
23 : use chem_def
24 : use utils_lib, only: is_bad_num
25 : use star_private_def
26 :
27 : implicit none
28 :
29 : private
30 : public :: tiny_mass
31 : public :: smallx
32 : public :: get_matrix_coeffs
33 :
34 : real(dp), parameter :: Xlim = 1d-14
35 : real(dp), parameter :: tiny_mass = 1d3 ! a kilogram
36 : real(dp), parameter :: tinyX = 1d-50
37 : real(dp), parameter :: smallX = 1d-20
38 :
39 : contains
40 :
41 0 : subroutine get_matrix_coeffs( &
42 : s, nz, nc, m, nzlo, nzhi, ih1, ihe4, pure_Coulomb, &
43 0 : dt, v_advection_max, tiny_C, diffusion_factor, &
44 0 : A, X, Z, rho_face, T_face, four_pi_r2_rho_face, &
45 0 : xm_face, cell_dm, dm_bar, dlnP_dr_face, dlnT_dr_face, dlnRho_dr_face, &
46 0 : r_face, r_mid, gamma_T_limit_coeffs_face, alfa_face, &
47 0 : rad_accel_face, log10_g_rad, g_rad, &
48 : min_T_for_radaccel, max_T_for_radaccel, &
49 0 : X_init, X_face, C, C_div_X, C_div_X_face, &
50 0 : e_ap, e_at, e_ar, e_ax, E_field_face, &
51 0 : g_ap, g_at, g_ar, g_ax, g_field_face, &
52 0 : v_advection_face, v_total_face, vlnP_face, vlnT_face, v_rad_face, &
53 0 : GT_face, D_self_face, AD_face, SIG_face, sigma_lnC, ierr)
54 :
55 : type (star_info), pointer :: s
56 : integer, intent(in) :: &
57 : nz, nc, m, nzlo, nzhi, ih1, ihe4
58 : logical, intent(in) :: pure_Coulomb
59 : real(dp), intent(in) :: &
60 : dt, v_advection_max, tiny_C, min_T_for_radaccel, max_T_for_radaccel
61 : real(dp), dimension(:), intent(in) :: &
62 : diffusion_factor, A, rho_face, T_face, four_pi_r2_rho_face, &
63 : xm_face, cell_dm, dm_bar, dlnP_dr_face, dlnT_dr_face, dlnRho_dr_face, &
64 : r_face, r_mid, gamma_T_limit_coeffs_face, alfa_face
65 : real(dp), dimension(:,:), intent(in) :: &
66 : Z, X_init, rad_accel_face, log10_g_rad, g_rad
67 : real(dp), dimension(:,:), intent(inout) :: X
68 :
69 : real(dp), dimension(:), intent(out) :: AD_face, &
70 : e_ap, e_at, e_ar, E_field_face, &
71 : g_ap, g_at, g_ar, g_field_face
72 : real(dp), dimension(:,:), intent(out) :: &
73 : X_face, C, C_div_X, C_div_X_face, v_advection_face, v_total_face, &
74 : vlnP_face, vlnT_face, v_rad_face, GT_face, e_ax, g_ax, D_self_face
75 : real(dp), dimension(:,:,:), intent(out) :: SIG_face, sigma_lnC
76 : integer, intent(out) :: ierr
77 :
78 : integer :: i, j, k, op_err
79 0 : real(dp) :: tmp, tinyX, dlamch, sfmin, &
80 0 : AD_dm_full_on, AD_dm_full_off, AD_boost_factor, sum_dm, &
81 0 : Vlimit_dm_full_on, Vlimit_dm_full_off, Vlimit, sigmax, &
82 0 : SIG_factor, GT_factor
83 0 : real(dp), dimension(m) :: C_face, Z_face, dC_dr_face
84 0 : real(dp), dimension(nc) :: total_diffusion_factor
85 0 : real(dp) :: dlnne_dr_face
86 :
87 : include 'formats'
88 :
89 0 : ierr = 0
90 0 : sfmin = dlamch('S')
91 :
92 0 : tinyX = 1d-50
93 0 : do k=nzlo,nzhi
94 0 : do j=1,nc
95 0 : X(j,k) = max(X(j,k),tinyX)
96 : end do
97 0 : tmp = sum(Z(1:nc,k)*X(1:nc,k)/A(1:nc))
98 0 : do j=1,nc
99 0 : C_div_X(j,k) = 1d0/(A(j)*tmp)
100 0 : C(j,k) = X(j,k)*C_div_X(j,k)
101 : end do
102 0 : C(m,k) = 1d0
103 0 : X(m,k) = A(m)/dot_product(A(1:nc),C(1:nc,k))
104 : end do
105 :
106 0 : Vlimit_dm_full_on = s% diffusion_Vlimit_dm_full_on*Msun
107 0 : Vlimit_dm_full_off = s% diffusion_Vlimit_dm_full_off*Msun
108 0 : Vlimit = s% diffusion_Vlimit
109 :
110 :
111 0 : !$OMP PARALLEL DO PRIVATE(k, j, i, total_diffusion_factor, op_err, C_face, Z_face, dC_dr_face, dlnne_dr_face, tmp) SCHEDULE(dynamic,2)
112 : do k = nzlo+1, nzhi
113 :
114 : ! Total diffusion scaling factor is product of
115 : ! diffusion_class_factor (const throughout the star), and
116 : ! extra_diffusion_factor (profile set in other_diffusion_factor)
117 : if(s% use_other_diffusion_factor) then
118 : do j=1,nc
119 : total_diffusion_factor(j) = diffusion_factor(j)*s% extra_diffusion_factor(j,k)
120 : end do
121 : else
122 : total_diffusion_factor(1:nc) = diffusion_factor(1:nc)
123 : end if
124 :
125 : call get1_CXZn_face( &
126 : k, nz, nc, m, nzlo, nzhi, C, X, Z, A, alfa_face, tiny_C, &
127 : four_pi_r2_rho_face(k)/dm_bar(k), dlnRho_dr_face(k), &
128 : C_face, X_face, Z_face, C_div_X_face, dC_dr_face, dlnne_dr_face)
129 :
130 : op_err = 0
131 : call get1_coeffs_face( &
132 : s, k, nz, nc, m, nzlo, nzhi, ih1, ihe4, pure_Coulomb, &
133 : rho_face(k), T_face(k), gamma_T_limit_coeffs_face(k), &
134 : four_pi_r2_rho_face(k), dm_bar(k), v_advection_max, tiny_C, sfmin, &
135 : dlnP_dr_face(k), dlnT_dr_face(k), dlnRho_dr_face(k), &
136 : s% grav(k), dlnne_dr_face, &
137 : Vlimit_dm_full_on, Vlimit_dm_full_off, Vlimit, xm_face(k), r_mid, dt, &
138 : A, X_face(:,k), Z_face, C_face, C_div_X_face(:,k), &
139 : total_diffusion_factor, rad_accel_face(:,k), &
140 : s% diffusion_use_cgs_solver, s% eta(k), &
141 : s% cgs_thermal_diffusion_eta_full_on, s% cgs_thermal_diffusion_eta_full_off, &
142 : (T_face(k) <= max_T_for_radaccel .and. T_face(k) >= min_T_for_radaccel), &
143 : v_advection_face(:,k), vlnP_face(:,k), vlnT_face(:,k), v_rad_face(:,k), &
144 : e_ap(k), e_at(k), e_ar(k), e_ax(:,k), &
145 : g_ap(k), g_at(k), g_ar(k), g_ax(:,k), &
146 : sigma_lnC(:,:,k), op_err)
147 : if (op_err /= 0) ierr = op_err
148 :
149 : if(s% diffusion_use_cgs_solver) then
150 : ! Electric Field from Iben & MacDonald solve.
151 : E_field_face(k) = &
152 : boltzm*T_face(k)*dlnT_dr_face(k)*e_at(k) + &
153 : amu*(s% grav(k))*e_ap(k) ! + e_ar(k) should be the radiation term, skipping for now
154 : do j=1,nc
155 : E_field_face(k) = E_field_face(k) + &
156 : boltzm*T_face(k)*e_ax(j,k)*(dC_dr_face(j)&
157 : &/C_face(j) + dlnne_dr_face)
158 : end do
159 : ! Because we solved for eE as the unknown rather than E in the vector containing diffusion velocities.
160 : E_field_face(k) = E_field_face(k)/qe
161 : g_field_face(k) = s% grav(k)
162 : else
163 : ! Electric and gravitational field from Thoul solve.
164 : E_field_face(k) = &
165 : e_ap(k)*dlnP_dr_face(k) + &
166 : e_at(k)*dlnT_dr_face(k) ! + e_ar(k)*...... skipping the radiation term for now
167 : g_field_face(k) = &
168 : g_ap(k)*dlnP_dr_face(k) + &
169 : g_at(k)*dlnT_dr_face(k) ! + g_ar(k)*...... skipping the radiation term for now
170 : do j=1,nc
171 : if (C_face(j) < 1d-20) cycle
172 : E_field_face(k) = E_field_face(k) + &
173 : e_ax(j,k)*dC_dr_face(j)/C_face(j)
174 : g_field_face(k) = g_field_face(k) + &
175 : g_ax(j,k)*dC_dr_face(j)/C_face(j)
176 : end do
177 : ! Convert to cgs units
178 : tmp = sum(Z_face(1:nc)*X_face(1:nc,k)/A(1:nc))
179 : g_field_face(k) = g_field_face(k)*(1.144d-40)*T_face(k)*tmp/(amu*amu)
180 : e_field_face(k) = e_field_face(k)*(1.144d-40)*T_face(k)*tmp/(qe*amu)
181 : end if
182 :
183 : do i = 1,nc
184 : v_total_face(i,k) = v_advection_face(i,k)
185 : do j = 1,nc
186 : v_total_face(i,k) = v_total_face(i,k) - sigma_lnC(i,j,k)*dC_dr_face(j)/C_face(j)
187 : end do
188 : end do
189 :
190 : end do
191 : !$OMP END PARALLEL DO
192 :
193 0 : if (ierr /= 0) return
194 0 : sum_dm = cell_dm(nzlo)
195 :
196 0 : AD_dm_full_on = s% diffusion_AD_dm_full_on*Msun
197 0 : AD_dm_full_off = s% diffusion_AD_dm_full_off*Msun
198 0 : AD_boost_factor = s% diffusion_AD_boost_factor
199 :
200 0 : SIG_factor = s% diffusion_SIG_factor
201 0 : GT_factor = s% diffusion_GT_factor
202 :
203 : !write(*,1) 'GT_factor SIG_factor', GT_factor, SIG_factor
204 :
205 0 : do k = nzlo+1, nzhi
206 : call get1_flow_coeffs( &
207 : k, nc, m, v_advection_face(:,k), v_advection_max, &
208 : SIG_factor, GT_factor, sigma_lnC(:,:,k), &
209 : four_pi_r2_rho_face(k), dm_bar(k), &
210 0 : C_div_X_face(:,k), GT_face(:,k), D_self_face(:,k), SIG_face(:,:,k))
211 0 : if (sum_dm >= AD_dm_full_off) then
212 0 : AD_face(k) = 0d0
213 : else
214 : sigmax = 0d0
215 0 : do j=1,nc
216 0 : if (SIG_face(j,j,k) > sigmax) sigmax = SIG_face(j,j,k)
217 : end do
218 0 : AD_face(k) = AD_boost_factor*sigmax
219 0 : if (sum_dm > AD_dm_full_on) &
220 : AD_face(k) = AD_face(k) * &
221 : (sum_dm - AD_dm_full_off)/&
222 0 : (AD_dm_full_on - AD_dm_full_off)
223 : !write(*,2) 'boost factor AD_face', k, AD_face(k)/sigmax, AD_face(k)
224 : end if
225 0 : sum_dm = sum_dm + cell_dm(k)
226 : end do
227 :
228 0 : do j=1,nc ! not used, but copy just for sake of plotting
229 0 : D_self_face(j,nzlo) = D_self_face(j,nzlo+1)
230 0 : v_advection_face(j,nzlo) = v_advection_face(j,nzlo+1)
231 0 : v_total_face(j,nzlo) = v_total_face(j,nzlo+1)
232 0 : vlnP_face(j,nzlo) = vlnP_face(j,nzlo+1)
233 0 : vlnT_face(j,nzlo) = vlnT_face(j,nzlo+1)
234 0 : v_rad_face(j,nzlo) = v_rad_face(j,nzlo+1)
235 0 : GT_face(j,nzlo) = GT_face(j,nzlo+1)
236 0 : do i=1,nc
237 0 : SIG_face(i,j,nzlo) = SIG_face(i,j,nzlo+1)
238 : end do
239 : end do
240 :
241 : end subroutine get_matrix_coeffs
242 :
243 :
244 0 : subroutine get1_coeffs_face( &
245 : s, k, nz, nc, m, nzlo, nzhi, ih1, ihe4, pure_Coulomb, &
246 : rho_face, T_face, gamma_T_limit_coeff_face, &
247 : four_pi_r2_rho_face, dm_bar, &
248 : v_advection_max, tiny_C, sfmin, &
249 : dlnP_dr_face, dlnT_dr_face, dlnRho_dr_face, grav, dlnne_dr_face, &
250 0 : Vlimit_dm_full_on, Vlimit_dm_full_off, Vlimit, xm_face, r_mid, dt, &
251 0 : A, X_face, Z_face, C_face, C_div_X_face, &
252 0 : diffusion_factor, rad_accel_face, &
253 : use_cgs_solver, eta, eta_on, eta_off, rad, &
254 0 : v_advection_face, vlnP_face, vlnT_face, v_rad_face, &
255 0 : e_ap, e_at, e_ar, e_ax, &
256 0 : g_ap, g_at, g_ar, g_ax, &
257 0 : sigma_lnC, ierr)
258 :
259 : type (star_info), pointer :: s
260 : integer, intent(in) :: k, nz, nc, m, nzlo, nzhi, ih1, ihe4
261 : logical, intent(in) :: pure_Coulomb
262 : real(dp), intent(in) :: rho_face, T_face, gamma_T_limit_coeff_face, &
263 : xm_face, r_mid(:), Vlimit_dm_full_on, Vlimit_dm_full_off, Vlimit, dt
264 : real(dp), intent(in) :: four_pi_r2_rho_face, dm_bar
265 : real(dp), intent(in) :: v_advection_max, tiny_C, sfmin
266 : real(dp), intent(in) :: dlnP_dr_face, dlnT_dr_face, &
267 : dlnRho_dr_face, grav, dlnne_dr_face
268 : real(dp), intent(in), dimension(:) :: &
269 : A, X_face, Z_face, C_face, C_div_X_face, &
270 : rad_accel_face, diffusion_factor
271 : logical, intent(in) :: use_cgs_solver, rad
272 : real(dp), intent(in) :: eta, eta_on, eta_off
273 : real(dp), intent(inout), dimension(:) :: &
274 : v_advection_face, vlnP_face, vlnT_face, v_rad_face ! (nc)
275 : real(dp), intent(out) :: e_ap, e_at, e_ar, g_ap, g_at, g_ar
276 : real(dp), intent(inout) :: e_ax(:), g_ax(:) ! (m)
277 : real(dp), intent(inout) :: sigma_lnC(:,:) ! (nc,nc)
278 : integer, intent(out) :: ierr
279 :
280 0 : real(dp), dimension(m) :: AP, AT, AR
281 0 : real(dp), dimension(m,m) :: kappa_st, Zdiff, Zdiff1, Zdiff2, AX
282 :
283 : include 'formats'
284 :
285 0 : ierr = 0
286 :
287 : call get1_burgers_coeffs( &
288 : s, k, nc, m, A, Z_face, X_face, C_face, &
289 : rho_face, T_face, pure_Coulomb, &
290 0 : kappa_st, Zdiff, Zdiff1, Zdiff2)
291 :
292 : call get1_gradient_coeffs( &
293 : k, m, sfmin, A, Z_face, X_face, C_face, rho_face, T_face, &
294 : use_cgs_solver, eta, eta_on, eta_off, &
295 : rad, rad_accel_face, kappa_st, Zdiff, Zdiff1, Zdiff2, &
296 : AP, AT, AR, AX, &
297 : e_ap, e_at, e_ar, e_ax, &
298 : g_ap, g_at, g_ar, g_ax, &
299 0 : ierr)
300 0 : if (ierr /= 0) return
301 :
302 : call get1_diffusion_velocities( &
303 : k, nc, m, nzlo, nzhi, AP, AT, AR, AX, rho_face, T_face, &
304 : dlnP_dr_face, dlnT_dr_face, dlnRho_dr_face, &
305 : grav, dlnne_dr_face, X_face, &
306 : Vlimit_dm_full_on, Vlimit_dm_full_off, Vlimit, xm_face, r_mid, s% dt, &
307 : gamma_T_limit_coeff_face, v_advection_max, diffusion_factor, &
308 : use_cgs_solver, &
309 0 : v_advection_face, vlnP_face, vlnT_face, v_rad_face, sigma_lnC)
310 :
311 0 : end subroutine get1_coeffs_face
312 :
313 :
314 0 : subroutine get1_CXZn_face( &
315 0 : k, nz, nc, m, nzlo, nzhi, C, X, Z, A, alfa_face, tiny_C, &
316 0 : d_dr_factor, dlnRho_dr_face, C_face, X_face, Z_face, C_div_X_face, &
317 0 : dC_dr_face, dlnne_dr_face)
318 : integer, intent(in) :: k, nc, m, nz, nzlo, nzhi
319 : real(dp), dimension(:,:), intent(in) :: C, X, Z ! (m,nz)
320 : real(dp), intent(in) :: A(:) ! (m) atomic number
321 : real(dp), intent(in) :: alfa_face(:), d_dr_factor, dlnRho_dr_face
322 : real(dp), intent(in) :: tiny_C
323 : real(dp), dimension(:), intent(out) :: C_face, Z_face, dC_dr_face ! (m)
324 : real(dp), intent(out) :: dlnne_dr_face
325 : real(dp), dimension(:,:), intent(out) :: X_face, C_div_X_face ! (m,nz)
326 : integer :: j
327 0 : real(dp) :: tmp, tmp1, tmp2, dlntmp_dr_face, alfa, beta
328 :
329 0 : alfa = alfa_face(k)
330 0 : beta = 1d0 - alfa
331 0 : do j = 1, m
332 0 : X_face(j,k) = alfa*X(j,k) + beta*X(j,k-1)
333 0 : Z_face(j) = alfa*Z(j,k) + beta*Z(j,k-1)
334 : end do
335 :
336 : ! "tmp" at the face of the zone
337 0 : tmp = sum(Z_face(1:nc)*X_face(1:nc,k)/A(1:nc))
338 0 : do j = 1, m
339 0 : C_div_X_face(j,k) = 1/(A(j)*tmp)
340 0 : C_face(j) = X_face(j,k)*C_div_X_face(j,k)
341 :
342 : ! Old way of calculating the derivative. Fixed below.
343 : ! dC_dr_face(j) = (X(j,k-1) - X(j,k))*C_div_X_face(j,k)*d_dr_factor
344 :
345 0 : dC_dr_face(j) = (C(j,k-1) - C(j,k))*d_dr_factor
346 : end do
347 :
348 : ! Calculate the electron number density ln gradient used in
349 : ! converting between the Thoul and Iben/MacDonald notations.
350 : ! This is accomplished by adding two logarithmic derivatives
351 : ! evaluated at the face of the zone, one of which we already
352 : ! know, and one of which must be calculated.
353 :
354 : ! "tmp" averaged over zone k
355 0 : tmp1 = sum(Z(1:nc,k)*X(1:nc,k)/A(1:nc))
356 : ! "tmp" averaged over zone k-1
357 0 : tmp2 = sum(Z(1:nc,k-1)*X(1:nc,k-1)/A(1:nc))
358 0 : dlntmp_dr_face = ((tmp2-tmp1)/tmp)*d_dr_factor
359 : ! ne = rho*tmp/amu, so...
360 0 : dlnne_dr_face = dlnRho_dr_face + dlntmp_dr_face
361 :
362 0 : end subroutine get1_CXZn_face
363 :
364 :
365 0 : subroutine get1_burgers_coeffs( &
366 0 : s, k, nc, m, A, Z, X, C, rho, T, pure_Coulomb, &
367 0 : kappa_st, Zdiff, Zdiff1, Zdiff2)
368 :
369 : use paquette_coeffs, only: paquette_coefficients
370 :
371 : type (star_info), pointer :: s
372 : integer, intent(in) :: k, nc, m
373 : real(dp), intent(in) :: rho, T
374 : logical, intent(in) :: pure_Coulomb
375 : real(dp), intent(in), dimension(:) :: A, X, Z, C ! (m)
376 : real(dp), intent(inout), dimension(:,:) :: &
377 : kappa_st, Zdiff, Zdiff1, Zdiff2 ! (m,m)
378 :
379 : integer :: i, j
380 0 : real(dp) :: ac, ni, cz, xij, ne, ao, lambdad, lambda, alfa, Gamlo, Gamhi
381 0 : real(dp), dimension(m) :: charge, na
382 0 : real(dp), dimension(m,m) :: cl, Ath, Ddiff, Kdiff, Kdiff2
383 0 : real(dp) :: Gamma, kappa_SM
384 0 : real(dp) :: Ddiff_Caplan(nc)
385 :
386 0 : do i = 1, nc
387 0 : charge(i) = max(1d0, Z(i)) ! assume some ionization
388 : end do
389 0 : charge(m) = Z(m)
390 :
391 0 : if (.not. pure_Coulomb) then ! use Paquette coeffs
392 : ! Get number densities (per cm^3)
393 0 : do i = 1, nc
394 0 : na(i) = rho*X(i)/(A(i)*amu)
395 : end do
396 0 : na(m) = 0.d0
397 0 : do i = 1, nc
398 0 : na(m) = na(m) + charge(i)*na(i)
399 : end do
400 : ! Compute resistance coefficients from Paquette&al (1986)
401 : call paquette_coefficients( &
402 0 : rho, T, m, A, charge, na, Ddiff, Kdiff, Zdiff, Zdiff1, Zdiff2, Ath)
403 :
404 0 : if(.not. s% diffusion_use_paquette .and. .not. s% use_other_diffusion_coefficients) then
405 0 : Gamma = s% gam(k)
406 0 : Gamlo = 3d0
407 0 : Gamhi = 10d0 ! for blending over from Stanton & Murillo coeffs to CBF coeffs at high Gamma
408 0 : if(Gamma < Gamlo .or. .not. s% diffusion_use_caplan) then
409 0 : call get_SM_coeffs(nc,m,rho,T,A,charge,na,Kdiff,Zdiff,Zdiff1,Zdiff2,kappa_SM)
410 : ! This must get called after Paquette because it doesn't calculate
411 : ! the electron entries (m). It leaves them untouched while calculating
412 : ! and changing all the ion-ion terms (1:nc), so after calling this
413 : ! routine we have all ion-ion coefficients from Stanton&Murillo and
414 : ! all ion-electron coefficients from Paquette&al.
415 0 : else if(Gamma >= Gamlo .and. Gamma <= Gamhi) then ! blend
416 0 : alfa = (Gamhi - Gamma)/(Gamhi - Gamlo)
417 0 : Kdiff2(:,:) = Kdiff(:,:) ! need to initialize so that (m,m) entries are set
418 0 : call get_SM_coeffs(nc,m,rho,T,A,charge,na,Kdiff,Zdiff,Zdiff1,Zdiff2,kappa_SM) ! alfa = 1
419 0 : call get_CBF_coeffs(nc,m,rho,T,A,charge,na,Gamma,Ddiff_Caplan,Kdiff2) ! alfa = 0
420 0 : Kdiff(:,:) = alfa*Kdiff(:,:) + (1d0-alfa)*Kdiff2(:,:)
421 : else
422 : ! use Caplan, Bauer, & Freeman coefficients at moderate to strong coupling
423 0 : call get_CBF_coeffs(nc,m,rho,T,A,charge,na,Gamma,Ddiff_Caplan,Kdiff)
424 : end if
425 : end if
426 :
427 0 : if(s% use_other_diffusion_coefficients) then
428 : call s% other_diffusion_coefficients( &
429 : s% id, k, nc, m, rho, T, A, X, Z, C, charge, na, &
430 0 : Ddiff, Kdiff, Zdiff, Zdiff1, Zdiff2, Ath)
431 : end if
432 :
433 : ! Unit conversion conveniently applies to both Paquette and Stanton&Murillo
434 0 : kappa_st(:,:) = Kdiff(:,:)/(1.41D-25*pow(T,-1.5D0)*na(m)*na(m))
435 : ! = kappa_st of eq 37, Thoul&al 1994
436 : return
437 : end if
438 :
439 : ! calculate density of electrons (ne) from mass density (rho):
440 : ac=0.d0
441 0 : do i=1, m
442 0 : ac=ac+a(i)*c(i)
443 : end do
444 0 : ne=rho/(mp*ac)
445 : ! calculate interionic distance (ao):
446 0 : ni=0.d0
447 0 : do i=1, nc
448 0 : ni=ni+c(i)*ne
449 : end do
450 0 : ao=pow(0.23873d0/ni,one_third)
451 : ! calculate debye length (lambdad):
452 0 : cz=0.d0
453 0 : do i=1, m
454 0 : cz=cz+c(i)*charge(i)*charge(i)
455 : end do
456 0 : lambdad=6.9010d0*sqrt(t/(ne*cz))
457 : ! calculate lambda to use in coulomb logarithm:
458 0 : lambda=max(lambdad, ao)
459 : ! calculate coulomb logarithms:
460 0 : do i=1, m
461 0 : do j=1, m
462 0 : xij=2.3939d3*t*lambda/abs(charge(i)*charge(j))
463 0 : cl(i,j)=0.81245d0*log1p(0.18769d0*pow(xij,1.2d0))
464 : end do
465 : end do
466 :
467 : ! set coeffs for pure Coulomb potential
468 0 : do i=1, m
469 0 : do j=1, m
470 0 : Zdiff(i,j) = 0.6d0
471 0 : Zdiff1(i,j) = 1.3d0
472 0 : Zdiff2(i,j) = 2d0
473 : kappa_st(i,j) = &
474 : cl(i,j)*sqrt(a(i)*a(j)/(a(i)+a(j)))* &
475 0 : c(i)*c(j)*charge(i)*charge(i)*charge(j)*charge(j)
476 : end do
477 : end do
478 :
479 0 : end subroutine get1_burgers_coeffs
480 :
481 :
482 0 : subroutine get1_gradient_coeffs( &
483 0 : k, m, sfmin, A, Z, X, C, rho, T, use_cgs_solver, eta, eta_on, eta_off, &
484 0 : rad, rad_accel, kappa_st, Zdiff, Zdiff1, Zdiff2, &
485 0 : AP, AT, AR, AX, &
486 0 : e_ap, e_at, e_ar, e_ax, &
487 0 : g_ap, g_at, g_ar, g_ax, &
488 : ierr)
489 : integer, intent(in) :: k, m
490 : real(dp), intent(in) :: sfmin, rho, T
491 : real(dp), intent(in) :: eta, eta_on, eta_off
492 : real(dp), intent(in), dimension(:) :: A, X, Z, C, rad_accel ! (m)
493 : logical, intent(in) :: use_cgs_solver, rad
494 : real(dp), dimension(:,:), intent(in) :: &
495 : kappa_st, Zdiff, Zdiff1, Zdiff2 ! (m,m)
496 : real(dp), dimension(:), intent(out) :: AP, AT, AR ! (m)
497 : real(dp), intent(inout) :: AX(:,:) ! (m,m)
498 : real(dp), intent(inout) :: e_ap, e_at, e_ar, e_ax(:) ! (m)
499 : real(dp), intent(inout) :: g_ap, g_at, g_ar, g_ax(:) ! (m)
500 : integer, intent(out) :: ierr
501 :
502 : integer :: i, j
503 0 : real(dp) :: charge(m), nd(m), Kdiff(m,m), alfa, beta
504 :
505 : ! For blending solves with and without thermal diffusion.
506 0 : real(dp) :: AP1(m), AT1(m), AR1(m), AX1(m,m)
507 0 : real(dp) :: e_ap1, e_at1, e_ar1, e_ax1(m)
508 0 : real(dp) :: AP2(m), AT2(m), AR2(m), AX2(m,m)
509 0 : real(dp) :: e_ap2, e_at2, e_ar2, e_ax2(m)
510 :
511 : include 'formats'
512 :
513 0 : ierr = 0
514 :
515 0 : do i=1,m-1
516 0 : charge(i) = max(1d0, Z(i))
517 : end do
518 0 : charge(m) = Z(m)
519 :
520 0 : if(use_cgs_solver) then ! Use the cgs solver
521 : ! Get number densities using info contained in X,A,Z,rho
522 0 : nd(1:m) = 0d0
523 0 : do i=1,m-1
524 0 : nd(i) = rho*X(i)/(A(i)*amu)
525 0 : nd(m) = nd(m) + nd(i)*charge(i) ! Electron Number Density satisfies charge neutrality
526 : end do
527 :
528 0 : Kdiff(:,:) = kappa_st(:,:)*(1.41D-25*pow(T,-1.5D0)*nd(m)*nd(m))
529 :
530 0 : if(eta < eta_on) then
531 : call solve_burgers_cgs_with_thermal(2*m+1,m,A,charge,nd,rad_accel,rad, &
532 : Kdiff, Zdiff, Zdiff1, Zdiff2, &
533 : AP,AT,AR,AX, &
534 0 : e_ap,e_at,e_ar,e_ax,ierr)
535 0 : else if(eta > eta_off) then
536 : call solve_burgers_cgs_no_thermal(m+1,m,A,charge,nd,rad_accel,rad, &
537 : Kdiff,AP,AT,AR,AX, &
538 0 : e_ap,e_at,e_ar,e_ax,ierr)
539 : else
540 : ! Call both and do a linear blend of all coefficients.
541 0 : alfa = (eta - eta_on)/(eta_off - eta_on) ! alfa = 1 means no thermal diffusion.
542 0 : beta = 1d0 - alfa ! beta = 1 means full thermal diffusion.
543 :
544 : call solve_burgers_cgs_no_thermal(m+1,m,A,charge,nd,rad_accel,rad, &
545 : Kdiff,AP1,AT1,AR1,AX1, &
546 0 : e_ap1,e_at1,e_ar1,e_ax1,ierr)
547 :
548 0 : if (ierr /= 0) then
549 : !return
550 0 : write(*,2) 'solve_burgers_cgs_no_thermal failed', k
551 0 : do i=1,m-1
552 0 : write(*,2) 'A X Z C', i, A(i), X(i), Z(i), C(i)
553 : end do
554 0 : call mesa_error(__FILE__,__LINE__,'get1_gradient_coeffs')
555 : end if
556 :
557 : call solve_burgers_cgs_with_thermal(2*m+1,m,A,charge,nd,rad_accel,rad, &
558 : Kdiff, Zdiff, Zdiff1, Zdiff2, &
559 : AP2,AT2,AR2,AX2, &
560 0 : e_ap2,e_at2,e_ar2,e_ax2,ierr)
561 :
562 : ! Check how much different the two solutions are.
563 : ! if( abs((AT1(3) - AT2(3))/AT1(3)) > 2d0 ) then ! 3 is index for Helium in basic.net
564 : ! print *, "Thermal diffusion changing temperature coefficient by more than factor of two."
565 : ! print *, "Relative difference: ", abs((AT1(3) - AT2(3))/AT1(3))
566 : ! end if
567 :
568 : ! Blending between the two solutions.
569 0 : do i = 1,m
570 0 : AP(i) = alfa*AP1(i) + beta*AP2(i)
571 0 : AT(i) = alfa*AT1(i) + beta*AT2(i)
572 0 : AR(i) = alfa*AR1(i) + beta*AR2(i)
573 0 : do j = 1,m
574 0 : AX(i,j) = alfa*AX1(i,j) + beta*AX2(i,j)
575 : end do
576 0 : e_ax(i) = alfa*e_ax1(i) + beta*e_ax2(i)
577 : end do
578 0 : e_ap = alfa*e_ap1 + beta*e_ap2
579 0 : e_at = alfa*e_at1 + beta*e_at2
580 0 : e_ar = alfa*e_ar1 + beta*e_ar2
581 :
582 : end if
583 : ! Gravity isn't being calculated by this version of the
584 : ! diffusion routine, so set the coefficients to zero.
585 0 : g_ap = 0d0
586 0 : g_at = 0d0
587 0 : g_ar = 0d0
588 0 : g_ax(1:m) = 0d0
589 :
590 0 : if (ierr /= 0) then
591 : !return
592 0 : write(*,2) 'solve_burgers_cgs failed', k
593 0 : do i=1,m-1
594 0 : write(*,2) 'A X Z C', i, A(i), X(i), Z(i), C(i)
595 : end do
596 0 : call mesa_error(__FILE__,__LINE__,'get1_gradient_coeffs')
597 : end if
598 : else ! Use the Thoul solver
599 : call do1_solve_thoul_hu( &
600 : 2*m+2, m, sfmin, A, charge, X, C, rad_accel, rad, &
601 : kappa_st, Zdiff, Zdiff1, Zdiff2, &
602 : AP, AT, AR, AX, &
603 : e_ap, e_at, e_ar, e_ax, &
604 : g_ap, g_at, g_ar, g_ax, &
605 0 : ierr)
606 :
607 0 : if (ierr /= 0) then
608 : !return
609 0 : write(*,2) 'do1_solve_thoul_hu failed', k
610 0 : do i=1,m-1
611 0 : write(*,2) 'A X Z C', i, A(i), X(i), Z(i), C(i)
612 : end do
613 0 : call mesa_error(__FILE__,__LINE__,'get1_gradient_coeffs')
614 : end if
615 : end if
616 :
617 :
618 0 : end subroutine get1_gradient_coeffs
619 :
620 :
621 0 : subroutine get1_diffusion_velocities( &
622 0 : k, nc, m, nzlo, nzhi, AP, AT, AR, AX, rho, T, &
623 0 : dlnP_dr, dlnT_dr, dlnRho_dr, grav, dlnne_dr, X_face, &
624 0 : Vlimit_dm_full_on, Vlimit_dm_full_off, Vlimit, xm_face, r_mid, dt, &
625 0 : limit_coeff, v_advection_max, diffusion_factor, use_cgs_solver, &
626 0 : vgt, vlnP, vlnT, vrad, sigma_lnC)
627 : integer, intent(in) :: k, nc, m, nzlo, nzhi
628 : real(dp), intent(in), dimension(:) :: AP, AT, AR, r_mid
629 : real(dp), intent(in) :: AX(:,:) ! (m,m)
630 : real(dp), intent(in) :: rho, T, limit_coeff, v_advection_max, &
631 : Vlimit_dm_full_on, Vlimit_dm_full_off, Vlimit, xm_face, dt
632 : real(dp), intent(in) :: dlnP_dr, dlnT_dr, dlnRho_dr, grav, dlnne_dr
633 : real(dp), intent(in) :: X_face(:)
634 : real(dp), intent(in) :: diffusion_factor(:)
635 : logical, intent(in) :: use_cgs_solver
636 : real(dp), intent(inout), dimension(:) :: vgt, vlnP, vlnT, vrad
637 : real(dp), intent(inout) :: sigma_lnC(:,:) ! (nc,nc)
638 :
639 : integer :: i, j, im
640 0 : real(dp) :: coef, coef_vrad, dv_im, dr, T2pt5, &
641 0 : vcross, vmax, vmax_limit, frac, alfa, beta
642 0 : real(dp) :: tau0 ! = 6d13*secyer, solar diffusion time (seconds)
643 : real(dp), parameter :: rho_unit = 1d2
644 : real(dp), parameter :: T_unit = 1d7
645 :
646 : include 'formats'
647 :
648 0 : if (limit_coeff <= 0) then
649 0 : vgt(:) = 0
650 0 : sigma_lnC(:,:) = 0
651 : return
652 : end if
653 :
654 0 : dr = r_mid(k-1) - r_mid(k)
655 0 : vcross = dr/dt
656 0 : if (xm_face >= Vlimit_dm_full_off .or. Vlimit <= 0d0) then
657 : vmax_limit = 1d99
658 : alfa = 0d0
659 0 : beta = 1d0
660 0 : else if (xm_face <= Vlimit_dm_full_on) then
661 0 : vmax_limit = vcross*Vlimit
662 0 : alfa = 1d0
663 0 : beta = 0d0
664 : else ! combine
665 : alfa = (xm_face - Vlimit_dm_full_off)/&
666 0 : (Vlimit_dm_full_on - Vlimit_dm_full_off)
667 0 : beta = 1d0 - alfa ! fraction of normal v when it is > vmax
668 0 : vmax_limit = vcross*Vlimit/alfa ! Want to scale to no limit at alfa = 0
669 : end if
670 :
671 0 : if(use_cgs_solver) then ! Converts coefficients to velocities
672 : ! assuming cgs routine.
673 0 : do i=1,nc
674 0 : do j=1,nc
675 0 : sigma_lnC(j,i) = -limit_coeff*diffusion_factor(i)*boltzm*T*AX(j,i)
676 : end do
677 0 : vlnP(i) = AP(i)*amu*grav*diffusion_factor(i)*limit_coeff
678 0 : vlnT(i) = AT(i)*boltzm*T*dlnT_dr*diffusion_factor(i)*limit_coeff
679 0 : vrad(i) = AR(i)*diffusion_factor(i)*limit_coeff ! AR already contains all constants.
680 0 : vgt(i) = vlnP(i) + vlnT(i) + vrad(i)
681 : end do
682 :
683 0 : do i = 1,nc
684 : ! Converting from Iben/MacDonald notation to Thoul
685 : ! notation using electron number density gradient.
686 0 : vgt(i) = vgt(i) - dlnne_dr*sum(sigma_lnC(i,1:nc))
687 0 : if (X_face(i) < 1d-15) vgt(i) = 0d0
688 : end do
689 :
690 : else ! converts coefficients to velocities assuming Thoul.
691 0 : tau0 = 6d13*secyer
692 0 : T2pt5 = pow(T/T_unit,2.5d0)
693 0 : coef = limit_coeff*Rsun*T2pt5/(rho/rho_unit)*(Rsun/tau0)
694 0 : coef_vrad = (limit_coeff/T)*Rsun*Rsun*T2pt5/(rho/rho_unit)/tau0
695 : ! converts to cgs units
696 : ! T^(5/2)/rho takes care of the last part from Thoul equations
697 : ! (21) and (28). Rsun/tau0 converts the velocity from units of
698 : ! Rsun/tau0 to cm/s. The extra Rsun accounts for the fact that
699 : ! MESA's gradients (e.g. dlnP_dr) have units of cm^-1, whereas
700 : ! Thoul's gradients in eqn (21) are Rsun^-1.
701 :
702 0 : do i=1,nc
703 0 : do j=1,nc
704 0 : sigma_lnC(j,i) = -diffusion_factor(i)*coef*AX(j,i)
705 : end do
706 0 : vlnP(i) = AP(i)*dlnP_dr*diffusion_factor(i)*coef
707 0 : vlnT(i) = AT(i)*dlnT_dr*diffusion_factor(i)*coef
708 0 : vrad(i) = AR(i)*diffusion_factor(i)*coef_vrad
709 0 : vgt(i) = vlnP(i) + vlnT(i) + vrad(i)
710 0 : if (X_face(i) < 1d-15) vgt(i) = 0d0
711 : end do
712 : end if
713 :
714 : ! final fixup for vgt of most abundant so it gives baryon conservation.
715 0 : im = maxloc(X_face(1:nc),dim=1)
716 0 : dv_im = -dot_product(X_face(1:nc), vgt(1:nc))/X_face(im)
717 0 : vgt(im) = vgt(im) + dv_im
718 :
719 0 : vmax = maxval(abs(vgt(1:nc)))
720 0 : if (vmax > v_advection_max) then
721 0 : frac = v_advection_max/vmax
722 0 : do i=1,nc
723 0 : vgt(i) = vgt(i)*frac
724 : ! Need to rescale the sigma terms by the same factor
725 : ! or diffusive equilibrium can break in WDs.
726 0 : do j=1,nc
727 0 : sigma_lnC(j,i) = sigma_lnC(j,i)*frac
728 : end do
729 : end do
730 : vmax = v_advection_max
731 : !write(*,3) 'vmax > v_advection_max', im, k, vmax, v_advection_max
732 : !call mesa_error(__FILE__,__LINE__,'get1_diffusion_velocities')
733 : end if
734 :
735 0 : if (alfa > 0d0 .and. vmax > vmax_limit) then
736 0 : frac = vmax_limit/vmax
737 0 : do i=1,nc
738 0 : vgt(i) = vgt(i)*frac
739 : ! Need to rescale the sigma terms by the same factor
740 : ! or diffusive equilibrium can break in WDs.
741 0 : do j=1,nc
742 0 : sigma_lnC(j,i) = sigma_lnC(j,i)*frac
743 : end do
744 : end do
745 : end if
746 :
747 : end subroutine get1_diffusion_velocities
748 :
749 :
750 0 : subroutine get1_flow_coeffs( &
751 : k, nc, m, &
752 0 : v_advection_face, v_advection_max, SIG_factor, GT_factor, &
753 0 : sigma_lnC_face, four_pi_r2_rho_face, &
754 0 : dm_bar, C_div_X_face, GT_face, D_self_face, SIG_face)
755 : integer, intent(in) :: k, nc, m
756 : real(dp), intent(in) :: &
757 : v_advection_max, SIG_factor, GT_factor, v_advection_face(:) ! (nc)
758 : real(dp), intent(in) :: sigma_lnC_face(:,:) ! (nc,nc)
759 : real(dp), intent(in) :: four_pi_r2_rho_face, dm_bar
760 : real(dp), intent(in), dimension(:) :: C_div_X_face ! (m)
761 : real(dp), intent(inout) :: GT_face(:) ! (nc)
762 : real(dp), intent(inout) :: D_self_face(:) ! (nc)
763 : real(dp), intent(inout) :: SIG_face(:,:) ! (nc,nc)
764 :
765 : integer :: i, j
766 0 : real(dp) :: c
767 :
768 : include 'formats'
769 :
770 0 : c = SIG_factor*four_pi_r2_rho_face*four_pi_r2_rho_face/dm_bar
771 0 : do i = 1, nc
772 0 : GT_face(i) = GT_factor*four_pi_r2_rho_face*v_advection_face(i)
773 0 : D_self_face(i) = sigma_lnC_face(i,i)
774 0 : do j = 1, nc
775 0 : SIG_face(i,j) = c*sigma_lnC_face(i,j)/C_div_X_face(j)
776 : end do
777 : end do
778 :
779 0 : end subroutine get1_flow_coeffs
780 :
781 :
782 : ! *************************************************************
783 : ! Original of this routine was written by Anne A. Thoul, at the Institute
784 : ! for Advanced Study, Princeton, NJ 08540.
785 : ! See Thoul et al., Ap.J. 421, p. 828 (1994)
786 :
787 : ! With modifications by Hali Hu for non Coulomb and rad levitation.
788 : ! *************************************************************
789 : ! This routine inverses the burgers equations.
790 : !
791 : ! The system contains N equations with N unknowns.
792 : ! The equations are: the M momentum equations,
793 : ! the M energy equations,
794 : ! two constraints: the current neutrality
795 : ! the zero fluid velocity.
796 : ! The unknowns are: the M diffusion velocities,
797 : ! the M heat fluxes,
798 : ! the electric field E
799 : ! the gravitational force g.
800 : !
801 : ! **************************************************
802 :
803 :
804 : ! comments from Anne's original version
805 : ! Inverse the system for each possible right-hand-side, i.e.,
806 : ! if alpha is the r.h.s., we obtain the coefficient A_p
807 : ! if nu ---------------------------------------- A_T
808 : ! if gamma(i,j) ----------------------------------- A_Cj
809 : !
810 : ! If I=1, we obtain the hydrogen diffusion velocity
811 : ! If I=2, ------------- helium ------------------
812 : ! If I=3,M-1, --------- heavy element -------------
813 : ! If I=M, ------------- electrons -----------------
814 : ! For I=M,2M, we get the heat fluxes
815 : ! For I=N-1, we get the electric field
816 : ! For I=N, we get the gravitational force g
817 :
818 :
819 0 : subroutine do1_solve_thoul_hu( &
820 0 : n, m, sfmin, a, z, x, c, rad_accel, rad, &
821 0 : kappa_st, Zdiff, Zdiff1, Zdiff2, &
822 0 : ap, at, ar, ax, &
823 0 : e_ap, e_at, e_ar, e_ax, &
824 0 : g_ap, g_at, g_ar, g_ax, &
825 : ierr)
826 :
827 : ! the parameter m is the number of fluids considered (ions+electrons)
828 : ! the parameter n is the number of equations (2*m+2).
829 : !
830 : ! the vectors a,z and x contain the atomic mass numbers,
831 : ! the charges (ionization), and the mass fractions, of the elements.
832 : ! note: since m is the electron fluid, its mass and charge must be
833 : ! a(m)=m_e/m_u
834 : ! z(m)=-1.
835 : !
836 : ! the array cl contains the values of the coulomb logarithms.
837 : ! the vector ap, at, and array ax contains the results for the diffusion
838 : ! coefficients.
839 :
840 : integer, intent(in) :: m,n
841 : real(dp), intent(in) :: sfmin
842 : real(dp), intent(in), dimension(:) :: A, Z, X, C, rad_accel ! (m)
843 : logical, intent(in) :: rad
844 : real(dp), intent(in), dimension(:,:) :: &
845 : kappa_st, Zdiff, Zdiff1, Zdiff2 ! (m,m)
846 : ! kappa_st from the resistance coefficient Kdiff with eq (37) Thoul&al.
847 : ! Zdiff, Zdiff1, Zdiff2 = arrays of resistance coefficients,
848 : real(dp), intent(inout), dimension(:) :: ap, at, ar ! (m)
849 : real(dp), intent(inout) :: ax(:,:) ! (m,m)
850 : real(dp), intent(inout) :: e_ap, e_at, e_ar, e_ax(:) ! (m)
851 : real(dp), intent(inout) :: g_ap, g_at, g_ar, g_ax(:) ! (m)
852 : integer, intent(out) :: ierr
853 :
854 0 : integer :: i, j, l, indx(n)
855 0 : real(dp) :: cc, ac, ko, f
856 0 : real(dp), dimension(m,m) :: xx, y, yy, k
857 0 : real(dp), dimension(n) :: alpha, nu, ga, beta
858 0 : real(dp), dimension(n,n) :: delta, gamma
859 :
860 : ! the vector c contains the concentrations
861 : ! cc is the total concentration: cc=sum(c_s)
862 : ! ac is proportional to the mass density: ac=sum(a_s c_s)
863 : ! the arrays xx,y,yy and k are various parameters which appear in
864 : ! burgers equations.
865 : ! the vectors and arrays alpha, nu, gamma, delta, and ga represent
866 : ! the "right- and left-hand-sides" of burgers equations, and later
867 : ! the diffusion coefficients.
868 :
869 : ! initialize:
870 :
871 0 : ierr = 0
872 0 : ko = 2d0
873 0 : indx(1:n) = 0
874 :
875 : ! calculate cc and ac:
876 :
877 0 : cc=sum(c(1:m))
878 0 : ac=dot_product(a(1:m),c(1:m))
879 :
880 : ! calculate the coefficients of the burgers equations
881 :
882 0 : do i=1,m
883 0 : do j=1,m
884 0 : xx(i,j)=a(j)/(a(i)+a(j))
885 0 : y(i,j)=a(i)/(a(i)+a(j))
886 0 : yy(i,j) = 3D0*y(i,j) + Zdiff1(i,j)*xx(i,j)*a(j)/a(i)
887 0 : k(i,j) = kappa_st(i,j)
888 : end do
889 : end do
890 :
891 : ! write the burgers equations and the two constraints as
892 : ! alpha_s dp + nu_s dt + sum_t(not ihe or m) gamma_st dc_t
893 : ! = sum_t delta_st w_t
894 :
895 0 : do i=1,m
896 0 : alpha(i)=c(i)/cc
897 0 : nu(i)=0d0
898 0 : gamma(i,1:n)=0d0
899 0 : if (rad) then
900 0 : beta(i) = -(amu/boltzm)*alpha(i)*a(i)*rad_accel(i)
901 : else
902 0 : beta(i) = 0d0
903 : end if
904 0 : do j=1,m
905 0 : if (j /= m) then ! HH: Include He gradient
906 0 : gamma(i,j) = -c(j)/cc
907 0 : if (j == i) gamma(i,j) = gamma(i,j) + 1d0
908 0 : gamma(i,j) = gamma(i,j)*c(i)/cc
909 : end if
910 : end do
911 : end do
912 :
913 0 : do i=m+1,n-2
914 0 : alpha(i)=0d0
915 0 : nu(i)=2.5d0*c(i-m)/cc
916 0 : beta(i) = 0d0
917 0 : gamma(i,1:n)=0d0
918 : end do
919 :
920 0 : alpha(n-1)=0d0
921 0 : nu(n-1)=0d0
922 0 : beta(n-1)=0d0
923 0 : gamma(n-1,1:n)=0d0
924 :
925 0 : alpha(n)=0d0
926 0 : nu(n)=0d0
927 0 : beta(n)=0d0
928 0 : gamma(n,1:n)=0d0
929 :
930 0 : delta(1:n,1:n) = 0d0
931 :
932 0 : do i=1,m
933 :
934 0 : do j=1,m
935 0 : if (j == i) then
936 0 : do l=1,m
937 0 : if (l /= i) then
938 0 : delta(i,j)=delta(i,j)-k(i,l)
939 : end if
940 : end do
941 : else
942 0 : delta(i,j)=k(i,j)
943 : end if
944 : end do
945 :
946 0 : do j=m+1,n-2
947 0 : if (j-m == i) then
948 0 : do l=1,m
949 0 : if (l /= i) &
950 0 : delta(i,j) = delta(i,j) + Zdiff(i,l)*xx(i,l)*k(i,l)
951 : end do
952 : else
953 0 : delta(i,j) = -Zdiff(i,j-m)*y(i,j-m)*k(i,j-m)
954 : end if
955 : end do
956 :
957 0 : delta(i,n-1)=c(i)*z(i)
958 :
959 0 : delta(i,n)=-c(i)*a(i)
960 :
961 : end do
962 :
963 0 : do i=m+1,n-2
964 :
965 0 : do j=1,m
966 0 : if (j == i-m) then
967 0 : do l=1,m
968 0 : if (l /= i-m) delta(i,j) = &
969 0 : delta(i,j) + 2.5D0*Zdiff(i-m,l)*xx(i-m,l)*k(i-m,l)
970 : end do
971 : else
972 0 : delta(i,j) = -(2.5d0*Zdiff(i-m,j))*xx(i-m,j)*k(i-m,j)
973 : end if
974 : end do
975 :
976 0 : do j=m+1,n-2
977 0 : if (j-m == i-m) then
978 0 : do l=1,m
979 0 : if (l /= i-m) delta(i,j) = delta(i,j) - &
980 0 : y(i-m,l)*k(i-m,l)*(0.8D0*Zdiff2(i-m,l)*xx(i-m,l)+yy(i-m,l))
981 : end do
982 0 : delta(i,j) = delta(i,j) - 0.4D0*Zdiff2(i-m,i-m)*k(i-m,i-m)
983 : else
984 : delta(i,j) = k(i-m,j-m)*xx(i-m,j-m)*y(i-m,j-m) * &
985 0 : (3D0 + Zdiff1(i-m,j-m) - 0.8D0*Zdiff2(i-m,j-m))
986 : end if
987 : end do
988 :
989 0 : delta(i,n-1:n)=0d0
990 :
991 : end do
992 :
993 0 : do j=1,m
994 0 : delta(n-1,j) = c(j)*z(j)
995 : end do
996 0 : delta(n-1,m+1:n) = 0d0
997 :
998 0 : do j=1,m
999 0 : delta(n,j) = c(j)*a(j)
1000 : end do
1001 0 : delta(n,m+1:n) = 0d0
1002 :
1003 0 : call dgetrf(n, n, delta, n, indx, ierr)
1004 0 : if (ierr /= 0) return
1005 :
1006 0 : call dgetrs( 'n', n, 1, delta, n, indx, alpha, n, ierr )
1007 0 : if (ierr /= 0) return
1008 :
1009 0 : call dgetrs( 'n', n, 1, delta, n, indx, nu, n, ierr )
1010 0 : if (ierr /= 0) return
1011 :
1012 0 : if (rad) then
1013 0 : call dgetrs( 'n', n, 1, delta, n, indx, beta, n, ierr )
1014 0 : if (ierr /= 0) return
1015 : end if
1016 :
1017 0 : do j=1,n
1018 0 : do i=1,n
1019 0 : ga(i)=gamma(i,j)
1020 : end do
1021 0 : call dgetrs( 'n', n, 1, delta, n, indx, ga, n, ierr )
1022 0 : if (ierr /= 0) return
1023 0 : do i=1,n
1024 0 : gamma(i,j)=ga(i)
1025 : end do
1026 : end do
1027 :
1028 0 : f = ko*ac*cc
1029 0 : do j=1,m
1030 0 : ap(j)=alpha(j)*f
1031 0 : at(j)=nu(j)*f
1032 0 : ar(j)=beta(j)*f
1033 0 : do i=1,m
1034 0 : ax(i,j)=gamma(i,j)*f
1035 : end do
1036 : end do
1037 :
1038 0 : e_ap=alpha(n-1)*f
1039 0 : g_ap=alpha(n)*f
1040 :
1041 0 : e_at=nu(n-1)*f
1042 0 : g_at=nu(n)*f
1043 :
1044 0 : e_ar=beta(n-1)*f
1045 0 : g_ar=beta(n)*f
1046 :
1047 0 : do i=1,m
1048 0 : e_ax(i)=gamma(n-1,i)*f
1049 0 : g_ax(i)=gamma(n,i)*f
1050 : end do
1051 :
1052 0 : end subroutine do1_solve_thoul_hu
1053 :
1054 :
1055 0 : subroutine solve_burgers_cgs_no_thermal( &
1056 0 : n, m, A, Z, nd, rad_accel, rad, &
1057 0 : Kdiff, ap, at, ar, ax, &
1058 0 : e_ap, e_at, e_ar, e_ax, ierr)
1059 :
1060 : ! nd = array of number densities
1061 : ! m = # of species including electrons
1062 : ! n = m+1 without thermal diffusion
1063 : ! is the number of equations and unknowns
1064 : ! Thermal diffusion off: (m-1) diffusion equations,
1065 : ! 2 conservation equations.
1066 : ! Thermal diffusion on: 2m diffusion equations (maybe 2*m-1)
1067 : ! 2 conservation equations.
1068 :
1069 : integer, intent(in) :: m,n
1070 : real(dp), intent(in), dimension(:) :: A, Z, nd, rad_accel ! (m)
1071 : logical, intent(in) :: rad
1072 : real(dp), intent(in), dimension(:,:) :: Kdiff ! (m,m)
1073 : real(dp), intent(inout), dimension(:) :: ap, at, ar ! (m)
1074 : real(dp), intent(inout) :: ax(:,:) ! (m,m)
1075 : real(dp), intent(inout) :: e_ap, e_at, e_ar, e_ax(:)
1076 : integer, intent(out) :: ierr
1077 :
1078 0 : integer :: i,j,l,indx(n)
1079 0 : real(dp), dimension(n) :: alpha, beta, nu, ga
1080 : ! Add in beta later for rad lev. ga is the temp holder for the
1081 : ! columns of gamm when doing matrix solve one column at a time.
1082 0 : real(dp), dimension(n,n) :: delta, gamm
1083 :
1084 : ! initialize
1085 0 : ierr = 0
1086 0 : indx(1:n) = 0
1087 0 : delta(1:n,1:n) = 0d0
1088 0 : alpha(1:n) = 0d0
1089 0 : beta(1:n) = 0d0
1090 0 : nu(1:n) = 0d0
1091 0 : ga(1:n) = 0d0
1092 0 : gamm(1:n,1:n) = 0d0
1093 :
1094 0 : e_ap = 0d0
1095 0 : e_at = 0d0
1096 0 : e_ar = 0d0
1097 0 : e_ax(1:m) = 0d0
1098 :
1099 : ! Assign the RHS Matrix multiplying the unknown quantities.
1100 : ! Right now this is for thermal diffusion off, assuming gravity
1101 : ! is a known, so there are m unknown diffusion velocities and
1102 : ! one unknown electric field, corresponding to (m-1) ion equaitions
1103 : ! (no electron equation) and 2 conservation equations.
1104 0 : do i = 1,(m-1)
1105 :
1106 0 : do j = 1,m
1107 0 : if (j == i) then
1108 0 : do l = 1,m
1109 0 : if (l /= j) then
1110 0 : delta(i,j) = delta(i,j) - Kdiff(i,l)
1111 : end if
1112 : end do
1113 : else
1114 0 : delta(i,j) = Kdiff(i,j)
1115 : end if
1116 : end do
1117 :
1118 0 : delta(i,n) = nd(i)*Z(i)
1119 :
1120 : end do
1121 :
1122 0 : do j = 1,m
1123 0 : delta(m,j) = A(j)*nd(j)
1124 0 : delta(n,j) = Z(j)*nd(j)
1125 : end do
1126 : ! Delta assignment finished.
1127 :
1128 : ! Assign LHS quantities. (Everything not assigned here stays zero thanks to earlier initialization)
1129 0 : do i = 1,(m-1)
1130 0 : alpha(i) = nd(i)*A(i)
1131 0 : nu(i) = nd(i)
1132 0 : gamm(i,i) = nd(i)
1133 0 : if (rad) then
1134 0 : beta(i) = -1d0*alpha(i)*amu*rad_accel(i)
1135 : end if
1136 : end do
1137 :
1138 : ! alpha,nu,gamm assignment finished
1139 :
1140 :
1141 : ! Invert the matrix and solve to get the contribution of each
1142 : ! vector/column on the LHS. This takes alpha, nu, ga from being
1143 : ! the knowns on the LHS to vectors containing information on
1144 : ! contributions for the unknowns. These then get assigned to the
1145 : ! output arrays ap, at, ax, which will then be multiplied by
1146 : ! gravity, k_b*T, and number density gradients to provide the
1147 : ! full solution for diffusion velocities.
1148 :
1149 0 : call dgetrf(n,n,delta,n,indx,ierr)
1150 :
1151 0 : if( ierr /= 0 ) then
1152 : ! print *, "Factoring failed!"
1153 : return
1154 : end if
1155 :
1156 0 : call dgetrs('N',n,1,delta,n,indx,alpha,n,ierr)
1157 0 : if( ierr /= 0 ) then
1158 : ! print *, "solve failed on alpha"
1159 : return
1160 : end if
1161 :
1162 0 : call dgetrs('N',n,1,delta,n,indx,nu,n,ierr)
1163 0 : if( ierr /= 0 ) then
1164 : ! print *, "solve failed on nu"
1165 : return
1166 : end if
1167 :
1168 0 : if (rad) then
1169 0 : call dgetrs('N',n,1,delta,n,indx,beta,n,ierr)
1170 0 : if( ierr /= 0 ) then
1171 : ! print *, "solve failed on beta"
1172 : return
1173 : end if
1174 : else
1175 0 : beta(1:n) = 0d0
1176 : end if
1177 :
1178 0 : do j=1,n
1179 0 : do i=1,n
1180 0 : ga(i) = gamm(i,j)
1181 : end do
1182 0 : call dgetrs('N',n,1,delta,n,indx,ga,n,ierr)
1183 0 : if( ierr /= 0 ) then
1184 : ! print *, "solve failed on gamma"
1185 : return
1186 : end if
1187 0 : do i=1,n
1188 0 : gamm(i,j) = ga(i)
1189 : end do
1190 : end do
1191 :
1192 : ! Assign the results of the matrix solve to the output
1193 : ! arrays/matrix.
1194 :
1195 0 : do j = 1,m
1196 0 : ap(j) = alpha(j)
1197 0 : at(j) = nu(j)
1198 0 : ar(j) = beta(j)
1199 0 : do i=1,m
1200 0 : ax(i,j) = gamm(i,j)
1201 : end do
1202 : end do
1203 :
1204 0 : e_ap = alpha(n)
1205 0 : e_at = nu(n)
1206 0 : e_ar = beta(n)
1207 0 : do i=1,m
1208 0 : e_ax(i) = gamm(n,i)
1209 : end do
1210 :
1211 0 : end subroutine solve_burgers_cgs_no_thermal
1212 :
1213 :
1214 0 : subroutine solve_burgers_cgs_with_thermal( &
1215 0 : n, m, A, Z, nd, rad_accel, rad, &
1216 0 : Kdiff, zdiff, zdiff1, zdiff2, &
1217 0 : ap, at, ar, ax, &
1218 0 : e_ap, e_at, e_ar, e_ax, ierr)
1219 :
1220 : ! nd = array of number densities
1221 : ! m = # of species including electrons
1222 : ! n = 2*m+1 with thermal diffusion
1223 : ! is the number of equations and unknowns
1224 : ! Thermal diffusion on: 2m-1 diffusion equations (m-1 momentum + m energy)
1225 : ! 2 conservation equations.
1226 : ! There is one less momentum equation because we are neglecting the electron
1227 : ! equation (which is wrong when electrons are degenerate),
1228 : ! and making up for losing that equation by treating grav as known.
1229 : ! We can't do the same with the energy equation because there are no more
1230 : ! unknowns that can be dropped, so this solver shouldn't really be used when
1231 : ! there is degeneracy. Having it in this form makes it easier to transition
1232 : ! between ideal gas (where this solver is valid) and solve_burgers_cgs
1233 : ! (which is much better when things are degenerate).
1234 :
1235 : integer, intent(in) :: m,n
1236 : real(dp), intent(in), dimension(:) :: A, Z, nd, rad_accel ! (m)
1237 : logical, intent(in) :: rad
1238 : real(dp), intent(in), dimension(:,:) :: Kdiff ! (m,m)
1239 : real(dp), intent(in), dimension(:,:) :: zdiff, zdiff1, zdiff2 ! (m,m)
1240 : real(dp), intent(inout), dimension(:) :: ap, at, ar ! (m)
1241 : real(dp), intent(inout) :: ax(:,:) ! (m,m)
1242 : real(dp), intent(inout) :: e_ap, e_at, e_ar, e_ax(:)
1243 : integer, intent(out) :: ierr
1244 :
1245 0 : integer :: i,j,l,indx(n), rightshift, downshift
1246 0 : real(dp), dimension(n) :: alpha, beta, nu, ga
1247 : ! Add in beta later for rad lev. ga is the temp holder for the
1248 : ! columns of gamm when doing matrix solve one column at a time.
1249 0 : real(dp), dimension(n,n) :: delta, gamm
1250 :
1251 : ! initialize
1252 0 : ierr = 0
1253 0 : indx(1:n) = 0
1254 0 : delta(1:n,1:n) = 0d0
1255 0 : alpha(1:n) = 0d0
1256 0 : beta(1:n) = 0d0
1257 0 : nu(1:n) = 0d0
1258 0 : ga(1:n) = 0d0
1259 0 : gamm(1:n,1:n) = 0d0
1260 :
1261 0 : e_ap = 0d0
1262 0 : e_at = 0d0
1263 0 : e_ar = 0d0
1264 0 : e_ax(1:m) = 0d0
1265 : ! Just a way to make it clear what's going on with the thermal terms in delta
1266 : ! by making block matrices and then shifting them into the proper position.
1267 : ! This makes the comparisons for checking the subdiagonals easier, as well
1268 : ! as indexing of the coefficients.
1269 0 : rightshift = m
1270 0 : downshift = m-1 ! Because electron momentum equation dropped.
1271 :
1272 : ! Assign the RHS Matrix multiplying the unknown quantities.
1273 : ! Assuming gravity is a known, so there are m unknown
1274 : ! diffusion velocities, m unknown heat flow vectors, and
1275 : ! one unknown electric field, corresponding to (m-1) momentum equations
1276 : ! (no electron equation), m energy equations, and 2 conservation equations.
1277 :
1278 : ! Terms corresponding to the momentum equations.
1279 0 : do i = 1,(m-1)
1280 :
1281 : ! Terms that multiply the diffusion velocities.
1282 0 : do j = 1,m
1283 0 : if (j == i) then
1284 0 : do l = 1,m
1285 0 : if (l /= j) then
1286 0 : delta(i,j) = delta(i,j) - Kdiff(i,l)
1287 : end if
1288 : end do
1289 : else
1290 0 : delta(i,j) = Kdiff(i,j)
1291 : end if
1292 : end do
1293 :
1294 : ! Terms that multiply the heat flow vectors.
1295 0 : do j = 1,m
1296 0 : if (j == i) then
1297 0 : do l = 1,m
1298 0 : if (l /= j) then
1299 0 : delta(i,j+rightshift) = delta(i,j+rightshift) + Kdiff(i,l)*zdiff(i,l)*A(l)/(A(i)+A(l))
1300 : end if
1301 : end do
1302 : else
1303 0 : delta(i,j+rightshift) = -1d0*Kdiff(i,j)*zdiff(i,j)*A(i)/(A(i) + A(j))
1304 : end if
1305 : end do
1306 :
1307 : ! Term multiplying the electric field.
1308 0 : delta(i,n) = nd(i)*Z(i)
1309 :
1310 : end do
1311 :
1312 : ! Terms corresponding to the energy equations.
1313 0 : do i = 1,m ! All these entries get shifted lower into the i+downshift position.
1314 :
1315 : ! Terms that multiply the diffusion velocities.
1316 0 : do j = 1,m
1317 0 : if (j == i) then
1318 0 : do l = 1,m
1319 0 : if (l /= j) then
1320 0 : delta(i+downshift,j) = delta(i+downshift,j) + 2.5d0*Kdiff(i,l)*zdiff(i,l)*A(l)/(A(i) + A(l))
1321 : end if
1322 : end do
1323 : else
1324 0 : delta(i+downshift,j) = -2.5d0*Kdiff(i,j)*zdiff(i,j)*A(j)/(A(i) + A(j))
1325 : end if
1326 : end do
1327 :
1328 : ! Terms that multiply the heat flow vectors.
1329 0 : do j = 1,m
1330 : ! Shift these entries to be in the lower-right half of the matrix.
1331 0 : if (j == i) then
1332 0 : delta(i+downshift,j+rightshift) = -0.4d0*Kdiff(i,j)*zdiff2(i,j) ! 1st term on RHS of eqn (90) in MESA3
1333 0 : do l=1,m ! second line of eqn (90) in MESA3
1334 0 : if (l /= j) then
1335 : delta(i+downshift,j+rightshift) = delta(i+downshift,j+rightshift) - &
1336 : Kdiff(i,l)* &
1337 : ( &
1338 : (3d0*pow2(A(i)) + pow2(A(l))*zdiff1(i,l))/pow2(A(i)+A(l)) + &
1339 : 0.8d0*zdiff2(i,l)*A(i)*A(l)/pow2(A(i) + A(l)) &
1340 0 : )
1341 : end if
1342 : end do
1343 : else ! third line of eqn (90) in MESA3
1344 : delta(i+downshift,j+rightshift) = Kdiff(i,j)* &
1345 : ( 3d0 + zdiff1(i,j) - 0.8d0*zdiff2(i,j) )* &
1346 0 : A(i)*A(j)/pow2(A(i)+A(j))
1347 : end if
1348 : end do
1349 :
1350 : ! Term multiplying the electric field. (doesn't appear in energy equations)
1351 0 : delta(i+downshift,n) = 0d0
1352 : end do
1353 :
1354 : ! Terms from the mass and charge conservation laws, the last two rows of the matrix.
1355 : ! These only involve the diffusion velocities, not heat flow or e-field.
1356 0 : do j = 1,m
1357 0 : delta(n-1,j) = A(j)*nd(j)
1358 0 : delta(n,j) = Z(j)*nd(j)
1359 : end do
1360 : ! Delta assignment finished.
1361 :
1362 : ! Assign LHS quantities. (Everything not assigned here stays zero thanks to earlier initialization)
1363 : ! LHS of momentum equations.
1364 0 : do i = 1,(m-1)
1365 0 : alpha(i) = nd(i)*A(i)
1366 0 : nu(i) = nd(i)
1367 0 : gamm(i,i) = nd(i)
1368 0 : if (rad) then
1369 0 : beta(i) = -1d0*alpha(i)*amu*rad_accel(i)
1370 : end if
1371 : end do
1372 : ! LHS of energy equations.
1373 0 : do i = 1,m
1374 0 : nu(i+downshift) = 2.5d0*nd(i)
1375 : end do
1376 :
1377 : ! alpha,nu,gamm assignment finished
1378 :
1379 :
1380 : ! Invert the matrix and solve to get the contribution of each
1381 : ! vector/column on the LHS. This takes alpha, nu, ga from being
1382 : ! the knowns on the LHS to vectors containing information on
1383 : ! contributions for the unknowns. These then get assigned to the
1384 : ! output arrays ap, at, ax, which will then be multiplied by
1385 : ! gravity, k_b*T, and number density gradients to provide the
1386 : ! full solution for diffusion velocities.
1387 :
1388 0 : call dgetrf(n,n,delta,n,indx,ierr)
1389 :
1390 0 : if( ierr /= 0 ) then
1391 : ! print *, "Factoring failed!"
1392 : return
1393 : end if
1394 :
1395 0 : call dgetrs('N',n,1,delta,n,indx,alpha,n,ierr)
1396 0 : if( ierr /= 0 ) then
1397 : ! print *, "solve failed on alpha"
1398 : return
1399 : end if
1400 :
1401 0 : call dgetrs('N',n,1,delta,n,indx,nu,n,ierr)
1402 0 : if( ierr /= 0 ) then
1403 : ! print *, "solve failed on nu"
1404 : return
1405 : end if
1406 :
1407 0 : if (rad) then
1408 0 : call dgetrs('N',n,1,delta,n,indx,beta,n,ierr)
1409 0 : if( ierr /= 0 ) then
1410 : ! print *, "solve failed on beta"
1411 : return
1412 : end if
1413 : else
1414 0 : beta(1:n) = 0d0
1415 : end if
1416 :
1417 0 : do j=1,n
1418 0 : do i=1,n
1419 0 : ga(i) = gamm(i,j)
1420 : end do
1421 0 : call dgetrs('N',n,1,delta,n,indx,ga,n,ierr)
1422 0 : if( ierr /= 0 ) then
1423 : ! print *, "solve failed on gamma"
1424 : return
1425 : end if
1426 0 : do i=1,n
1427 0 : gamm(i,j) = ga(i)
1428 : end do
1429 : end do
1430 :
1431 :
1432 : ! Assign the results of the matrix solve from the diffusion
1433 : ! velocity part of the arrays to the solution vectors.
1434 :
1435 0 : do j = 1,m
1436 0 : ap(j) = alpha(j)
1437 0 : at(j) = nu(j)
1438 0 : ar(j) = beta(j)
1439 0 : do i=1,m
1440 0 : ax(i,j) = gamm(i,j)
1441 : end do
1442 : end do
1443 :
1444 : ! Entries m+1-2m are heat flow vectors, so those don't get output.
1445 : ! Skip to final entry (n) for electric field.
1446 0 : e_ap = alpha(n)
1447 0 : e_at = nu(n)
1448 0 : e_ar = beta(n)
1449 0 : do i=1,m
1450 0 : e_ax(i) = gamm(n,i)
1451 : end do
1452 :
1453 0 : end subroutine solve_burgers_cgs_with_thermal
1454 :
1455 :
1456 : ! Calculate coefficients given in Appendix C.3 of Stanton & Murillo, PR E 93, 043203 (2016)
1457 0 : subroutine get_SM_coeffs(nc,m,rho,T,A,Z,nd,Kdiff,zdiff,zdiff1,zdiff2,kappa)
1458 : integer, intent(in) :: nc, m
1459 : real(dp), intent(in) :: rho, T
1460 : real(dp), dimension(:), intent(in) :: A, Z, nd ! m
1461 : real(dp), intent(out) :: kappa ! ion separation relative to electron screening length.
1462 : real(dp), dimension(:,:), intent(inout) :: Kdiff,zdiff,zdiff1,zdiff2 ! (m,m) but electron entries shouldn't be used.
1463 : ! Only the ion terms are modified by this routine right now.
1464 :
1465 0 : real(dp), dimension(nc,nc) :: g_plasma, mu
1466 0 : real(dp), dimension(nc,nc,2,3) :: Keff, Omega ! Indexed by all species and different orders
1467 : ! Fitting coefficients from Stanton & Murillo
1468 0 : real(dp), dimension(2,3) :: a1,a2,a3,a4,a5,b0,b1,b2,b3,b4
1469 0 : real(dp) :: lambda ! The screening length
1470 : integer :: i,j,facmo,no,mo ! Last two are used for different orders of collisions.
1471 :
1472 0 : real(dp) :: lgp, gp1, gp2, gp3, gp4, gp5, kbT32, tmp
1473 :
1474 : ! Initialize all to 0, since some entries never get set or used.
1475 0 : a1(:,:) = 0d0
1476 0 : a2(:,:) = 0d0
1477 0 : a3(:,:) = 0d0
1478 0 : a4(:,:) = 0d0
1479 0 : a5(:,:) = 0d0
1480 0 : b0(:,:) = 0d0
1481 0 : b1(:,:) = 0d0
1482 0 : b2(:,:) = 0d0
1483 0 : b3(:,:) = 0d0
1484 0 : b4(:,:) = 0d0
1485 :
1486 : ! Table IV, coefficients for fits (C23)-(C24)
1487 0 : a1(1,1) = 1.4660d0
1488 0 : a1(1,2) = 0.52094d0
1489 0 : a1(1,3) = 0.30346d0
1490 0 : a1(2,2) = 0.85401d0
1491 :
1492 0 : a2(1,1) = -1.7836d0
1493 0 : a2(1,2) = 0.25153d0
1494 0 : a2(1,3) = 0.23739d0
1495 0 : a2(2,2) = -0.22898d0
1496 :
1497 0 : a3(1,1) = 1.4313d0
1498 0 : a3(1,2) = -1.1337d0
1499 0 : a3(1,3) = -0.62167d0
1500 0 : a3(2,2) = -0.60059d0
1501 :
1502 0 : a4(1,1) = -0.55833d0
1503 0 : a4(1,2) = 1.2155d0
1504 0 : a4(1,3) = 0.56110d0
1505 0 : a4(2,2) = 0.80591d0
1506 :
1507 0 : a5(1,1) = 0.061162d0
1508 0 : a5(1,2) = -0.43784d0
1509 0 : a5(1,3) = -0.18046d0
1510 0 : a5(2,2) = -0.30555d0
1511 :
1512 0 : b0(1,1) = 0.081033d0
1513 0 : b0(1,2) = 0.20572d0
1514 0 : b0(1,3) = 0.68375d0
1515 0 : b0(2,2) = 0.43475d0
1516 :
1517 0 : b1(1,1) = -0.091336d0
1518 0 : b1(1,2) = -0.16536d0
1519 0 : b1(1,3) = -0.38459d0
1520 0 : b1(2,2) = -0.21147d0
1521 :
1522 0 : b2(1,1) = 0.051760d0
1523 0 : b2(1,2) = 0.061572d0
1524 0 : b2(1,3) = 0.10711d0
1525 0 : b2(2,2) = 0.11116d0
1526 :
1527 0 : b3(1,1) = -0.50026d0
1528 0 : b3(1,2) = -0.12770d0
1529 0 : b3(1,3) = 0.10649d0
1530 0 : b3(2,2) = 0.19665d0
1531 :
1532 0 : b4(1,1) = 0.17044d0
1533 0 : b4(1,2) = 0.066993d0
1534 0 : b4(1,3) = 0.028760d0
1535 0 : b4(2,2) = 0.15195d0
1536 :
1537 : ! Get the screening length
1538 0 : call lam_SM(nc,m,rho,T,Z,nd,lambda,kappa)
1539 :
1540 : ! Calculate g_plasma
1541 0 : do i=1,nc
1542 0 : do j=1,nc
1543 0 : g_plasma(i,j) = Z(i)*Z(j)*qe*qe/(lambda*boltzm*T)
1544 : end do
1545 : end do
1546 :
1547 : ! Calculate Keff using Eqns (C22-C24)
1548 0 : do i=1,nc
1549 0 : do j=1,nc
1550 : ! cache g powers
1551 0 : lgp = safe_log(g_plasma(i,j))
1552 0 : gp1 = g_plasma(i,j)
1553 0 : gp2 = pow2(gp1)
1554 0 : gp3 = pow3(gp1)
1555 0 : gp4 = pow4(gp1)
1556 0 : gp5 = pow5(gp1)
1557 0 : if( g_plasma(i,j) < 0d0) then ! Don't calculate for attractive potentials, set to 0
1558 0 : Keff(i,j,:,:) = 0d0
1559 0 : else if( g_plasma(i,j) < 1d0) then ! Use eqn C23 for weakly coupled
1560 0 : do no=1,2
1561 0 : do mo=1,3 ! Implementing the (m-1)! term with a simple if statement.
1562 0 : if(mo == 3) then
1563 : facmo = 2 ! (3-1)!
1564 : else
1565 0 : facmo = 1 ! (1-1)! and (2-1)!
1566 : end if
1567 : Keff(i,j,no,mo) = (-1d0*no/4d0)*facmo*safe_log( &
1568 : a1(no,mo)*gp1 &
1569 : + a2(no,mo)*gp2 &
1570 : + a3(no,mo)*gp3 &
1571 : + a4(no,mo)*gp4 &
1572 0 : + a5(no,mo)*gp5)
1573 : end do
1574 : end do
1575 : else ! Use eqn C24 for strongly coupled
1576 0 : do no=1,2
1577 0 : do mo=1,3
1578 : Keff(i,j,no,mo) = &
1579 : (b0(no,mo) + b1(no,mo)*lgp + b2(no,mo)*pow2(lgp) ) / &
1580 0 : (1d0 + b3(no,mo)*gp1 + b4(no,mo)*gp2 )
1581 : end do
1582 : end do
1583 : end if
1584 : end do
1585 : end do
1586 :
1587 : ! Calculate the collision integrals using (C19)
1588 0 : kBT32 = pow(boltzm*T,1.5d0)
1589 0 : do i=1,nc
1590 0 : do j=1,nc
1591 0 : mu(i,j) = amu*A(i)*A(j)/(A(i) + A(j)) ! Reduced mass for collision
1592 0 : tmp = sqrt(2d0*pi/(mu(i,j))) * (pow2(Z(i)*Z(j)*qe*qe)/kBT32)
1593 0 : do no=1,2
1594 0 : do mo=1,3
1595 0 : Omega(i,j,no,mo) = tmp * Keff(i,j,no,mo)
1596 : end do
1597 : end do
1598 : end do
1599 : end do
1600 :
1601 : ! Resistance coefficient is (16/3) ni nj mu Omega11
1602 : ! (compare Paquette eqn 22 with SM eq B8)
1603 : ! The zdiffs are given in terms of collision integrals in Paquette eqns (14-16), (23-25)
1604 0 : do i=1,nc
1605 0 : do j=1,nc
1606 0 : Kdiff(i,j) = (16d0/3d0)*nd(i)*nd(j)*mu(i,j)*Omega(i,j,1,1)
1607 0 : zdiff(i,j) = 1d0 - 2d0*Omega(i,j,1,2)/(5d0*Omega(i,j,1,1))
1608 0 : zdiff1(i,j) = 2.5d0 + (2d0*Omega(i,j,1,3) - 10d0*Omega(i,j,1,2))/(5d0*Omega(i,j,1,1))
1609 0 : zdiff2(i,j) = Omega(i,j,2,2)/Omega(i,j,1,1)
1610 : end do
1611 : end do
1612 :
1613 : ! Note that SM don't give fits for attractive potentials, so we haven't
1614 : ! touched the electron entries. They exit this routine unchanged, so they
1615 : ! either need to be initialized before this routine is called or somehow
1616 : ! calculated later.
1617 :
1618 0 : end subroutine get_SM_coeffs
1619 :
1620 : ! Screening Length according to Stanton & Murillo
1621 0 : subroutine lam_SM(nc,m,rho,T,Z,nd,lam_eff,kappa)
1622 : integer, intent(in) :: nc, m
1623 : real(dp), intent(in) :: rho, T
1624 : real(dp), dimension(:), intent(in) :: Z, nd ! m. charges, number densities of all species
1625 : real(dp), intent(out) :: lam_eff
1626 : real(dp), intent(out) :: kappa ! ion separation relative to electron screening length.
1627 :
1628 0 : real(dp) :: tiny_n, ne, EF, lam_e, lam_sum, rhotot, ni_tot
1629 0 : real(dp), dimension(nc) :: ai, gam_is, lam_i ! ion sphere radius, coupling, screening for each type of ion
1630 : integer :: i
1631 :
1632 0 : tiny_n = 1d-20 ! g/cc
1633 0 : ne = nd(m)
1634 : ! Electron Fermi energy
1635 0 : EF = (hbar*hbar*pow(3d0*pi*pi*ne,2d0/3d0))/(2d0*me)
1636 : ! Electron screening length accounting for degeneracy correction
1637 0 : lam_e = pow(pi4*qe*qe*ne/sqrt(pow2(boltzm*T) + pow2((2d0/3d0)*EF)),-0.5d0)
1638 :
1639 : ! Compute kappa
1640 0 : ni_tot = 0d0
1641 0 : do i = 1,nc
1642 0 : ni_tot = ni_tot + nd(i)
1643 : end do
1644 0 : kappa = pow(3d0/(pi4*ni_tot),one_third)/lam_e
1645 :
1646 0 : rhotot = 0d0
1647 0 : do i = 1,nc
1648 : ! rhotot is a CHARGE density. Just trying to follow the notation of SM eq (34)
1649 0 : rhotot = rhotot + Z(i)*qe*nd(i) ! = qe*ne?
1650 : end do
1651 0 : do i = 1,nc
1652 0 : ai(i) = pow(Z(i)*qe/(four_thirds_pi*rhotot), one_third)
1653 0 : gam_is(i) = Z(i)*Z(i)*qe*qe/(ai(i)*boltzm*T)
1654 : ! Number densities that are 0 or tiny cause screening length to diverge.
1655 : ! This is physical; nothing is there to screen anything.
1656 : ! But numerically, I don't want to rely on fortran to handle division by zero and infinity,
1657 : ! so just set these screening lengths to be huge and they won't
1658 : ! contribute anything to overall screening length.
1659 0 : if(nd(i) < tiny_n) then
1660 0 : lam_i(i) = 1d99 ! cm. This won't contribute to any screening.
1661 : else
1662 0 : lam_i(i) = sqrt(boltzm*T/( pi4*pow2(Z(i)*qe)*nd(i) ))
1663 : end if
1664 : end do
1665 :
1666 0 : lam_sum = 1d0/pow2(lam_e) ! The electron part of the screening length.
1667 0 : do i = 1,nc ! Sum over all the ions.
1668 0 : lam_sum = lam_sum + pow( pow2(lam_i(i))*(1d0+3d0*gam_is(i)), -1d0)
1669 : end do
1670 0 : lam_eff = pow(lam_sum,-0.5d0)
1671 0 : end subroutine lam_SM
1672 :
1673 : ! Calculate coefficients given in Section 4 of Caplan, Bauer, & Freeman 2022
1674 0 : subroutine get_CBF_coeffs(nc,m,rho,T,A,Z,nd,Gamma,Ddiff,Kdiff)
1675 : integer, intent(in) :: nc, m
1676 : real(dp), intent(in) :: rho, T, Gamma
1677 : real(dp), dimension(:), intent(in) :: A, Z, nd ! m
1678 : real(dp), dimension(:), intent(out) :: Ddiff ! nc
1679 : real(dp), dimension(:,:), intent(inout) :: Kdiff ! (m,m), but only touch (nc,nc)
1680 :
1681 0 : real(dp) :: kappa, DstarOCP, omegap, ai, lam_e, Docp
1682 0 : real(dp) :: Abar, Zbar, ni_sum, Zp6bar, Kfac
1683 : integer :: i, j
1684 :
1685 0 : ni_sum = sum(nd(1:nc)) ! total ion number density (excludes electron entry m)
1686 0 : Abar = sum(A(1:nc)*nd(1:nc))/ni_sum
1687 0 : Zbar = sum(Z(1:nc)*nd(1:nc))/ni_sum
1688 0 : Zp6bar = 0d0 ! average of Z^0.6, needed for constructing Kij
1689 0 : do j = 1,nc
1690 0 : Zp6bar = Zp6bar + pow(Z(j),0.6d0)*nd(j)/ni_sum
1691 : end do
1692 :
1693 : ! Set kappa as input for Dstar calculation,
1694 : ! also sets omegap, ai, lam_e for converting to cgs units later.
1695 0 : call kappa_CBF(nc,m,rho,T,Abar,Zbar,nd,omegap,ai,lam_e,kappa)
1696 :
1697 : ! Calculate Dstar using fit of Eqn (5)
1698 0 : call get_Dstar_OCP(Gamma,kappa,DstarOCP)
1699 0 : Docp = DstarOCP*ai*ai*omegap ! convert from dimensionless to cgs units
1700 :
1701 : ! Eqn (11)
1702 0 : do j = 1,nc
1703 0 : Ddiff(j) = Docp*pow(Z(j)/Zbar,-0.6d0)
1704 : end do
1705 :
1706 : ! Translate into Kij coefficients for Burgers equations
1707 0 : Kfac = boltzm*T/(ni_sum*Zp6bar*pow(Zbar,0.6d0)*Docp)
1708 0 : do i = 1,nc
1709 0 : do j = 1,nc
1710 0 : Kdiff(i,j) = nd(i)*nd(j)*pow(Z(i)*Z(j),0.6d0)*Kfac
1711 : end do
1712 : end do
1713 :
1714 0 : end subroutine get_CBF_coeffs
1715 :
1716 0 : subroutine get_Dstar_OCP(Gamma,kappa,DstarOCP)
1717 : real(dp), intent(in) :: Gamma, kappa
1718 : real(dp), intent(out) :: DstarOCP
1719 :
1720 0 : real(dp) :: a0, a1, a2, b0, b1, b2, c0, c1, c2, c3
1721 0 : real(dp) :: Ak, Bk, Ck
1722 :
1723 0 : a0 = 1.55973d0
1724 0 : a1 = 1.10941d0
1725 0 : a2 = 1.36909d0
1726 :
1727 0 : b0 = 0.0070782d0
1728 0 : b1 = 0.80499d0
1729 0 : b2 = 4.53523d0
1730 :
1731 0 : c0 = 2.20689d0
1732 0 : c1 = 1.351594d0
1733 0 : c2 = 1.57138d0
1734 0 : c3 = 3.34187d0
1735 :
1736 0 : Ak = sqrt(pi/3d0)*(a0 + a1*pow(kappa,a2))
1737 0 : Bk = b0*exp(-b1*pow(kappa,b2))
1738 0 : Ck = c0 + c1*bitsafe_erf_fit(c2*pow(kappa,c3))
1739 :
1740 : ! Eqn (5)
1741 0 : DstarOCP = sqrt(pi/3d0)*Ak*pow(Gamma,-2.5d0)*exp(-Bk*Gamma)/log(1d0 + Ck*pow(Gamma,-1.5d0)/sqrt(3d0))
1742 :
1743 0 : end subroutine get_Dstar_OCP
1744 :
1745 : ! Screening Length for CBF diffusion coefficients.
1746 : ! Only electron screening, accounts for potentially relativistic electron degeneracy.
1747 0 : subroutine kappa_CBF(nc,m,rho,T,Abar,Zbar,nd,omegap,ai,lam_e,kappa)
1748 : integer, intent(in) :: nc, m
1749 : real(dp), intent(in) :: rho, T, Abar, Zbar
1750 : real(dp), dimension(:), intent(in) :: nd ! m. number densities of all species
1751 : real(dp), intent(out) :: omegap, ai, lam_e
1752 : real(dp), intent(out) :: kappa ! ion separation relative to electron screening length.
1753 :
1754 0 : real(dp) :: tiny_n, ne, EF, kF, lam_e_SM, lam_e_rel, ni_tot
1755 : integer :: i
1756 :
1757 0 : tiny_n = 1d-20 ! g/cc
1758 0 : ne = nd(m) ! electron number density
1759 : ! Electron Fermi energy
1760 0 : EF = (hbar*hbar*pow(3d0*pi*pi*ne,2d0/3d0))/(2d0*me)
1761 : ! Electron screening length accounting for degeneracy correction
1762 0 : lam_e_SM = pow(pi4*qe*qe*ne/sqrt(pow2(boltzm*T) + pow2((2d0/3d0)*EF)),-0.5d0)
1763 :
1764 0 : if (rho < 1d6) then
1765 0 : lam_e = lam_e_SM
1766 : else ! calculate and use relativistic screening length if rho > 1d6
1767 0 : kF = pow(3d0*pi*pi*ne,one_third) ! electron Fermi momentum
1768 0 : lam_e_rel = 1d0/(2d0*kF*sqrt(fine/pi)) ! fine = fine structure constant ~1/137
1769 0 : lam_e = min(lam_e_SM,lam_e_rel) ! min for continuity
1770 : end if
1771 :
1772 : ! Compute kappa
1773 0 : ni_tot = 0d0
1774 0 : do i = 1,nc
1775 0 : ni_tot = ni_tot + nd(i)
1776 : end do
1777 0 : omegap = sqrt(pi4*ni_tot*Zbar*Zbar*qe*qe/(Abar*amu)) ! plasma frequency
1778 0 : ai = pow(3d0/(pi4*ni_tot),one_third)
1779 0 : kappa = ai/lam_e
1780 :
1781 0 : end subroutine kappa_CBF
1782 :
1783 : ! crlibm does not include error function, so implement an
1784 : ! erf fit good to 2.5e-5 (Abramowitz and Stegun 1964)
1785 : ! in terms of other bitsafe functions.
1786 0 : real(dp) function bitsafe_erf_fit(x)
1787 : real(dp) :: x
1788 0 : real(dp) :: t, p, a1, a2, a3
1789 :
1790 0 : p = 0.47047d0
1791 0 : a1 = 0.3480242d0
1792 0 : a2 = -0.0958798d0
1793 0 : a3 = 0.7478556d0
1794 :
1795 0 : t = 1d0/(1d0 + p*x)
1796 :
1797 0 : bitsafe_erf_fit = 1d0 - (a1*t + a2*t*t + a3*t*t*t)*exp(-x*x)
1798 0 : end function bitsafe_erf_fit
1799 :
1800 : end module diffusion_support
|
