Line data Source code
1 : ! ***********************************************************************
2 : !
3 : ! Copyright (C) 2010-2021 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 thermohaline
21 :
22 : use const_def, only: dp, one_third, pi, qe, amu, boltzm, crad, clight
23 : use num_lib
24 : use utils_lib
25 : use auto_diff
26 :
27 : implicit none
28 :
29 : private
30 : public :: get_D_thermohaline
31 :
32 : contains
33 :
34 : !> Computes the diffusivity of thermohaline mixing when the
35 : !! thermal gradient is stable and the composition gradient is unstable.
36 : !!
37 : !! @param thermohaline_option A string specifying which thermohaline prescription to use.
38 : !! @param grada Adiabatic gradient dlnT/dlnP
39 : !! @param gradr Radiative temperature gradient dlnT/dlnP, equals the actual gradient because there's no convection
40 : !! @param T Temperature
41 : !! @param opacity opacity
42 : !! @param rho Density
43 : !! @param Cp Heat capacity at constant pressure
44 : !! @param gradL_composition_term dlnMu/dlnP where Mu is the mean molecular weight.
45 : !! @param iso The index of the species that drives thermohaline mixing.
46 : !! @param XH1 Mass fraction of H1.
47 : !! @param thermohaline_coeff Free parameter multiplying the thermohaline diffusivity.
48 : !! @param D_thrm Output, diffusivity.
49 : !! @param ierr Output, error index.
50 0 : subroutine get_D_thermohaline(thermohaline_option, &
51 : grada, gradr, T, opacity, rho, Cp, gradL_composition_term, &
52 : iso, XH1, thermohaline_coeff, D_thrm, ierr)
53 : character(len=*), intent(in) :: thermohaline_option
54 : real(dp), intent(in) :: &
55 : grada, gradr, T, opacity, rho, Cp, gradL_composition_term, XH1, &
56 : thermohaline_coeff
57 : integer, intent(in) :: iso
58 : real(dp), intent(out) :: D_thrm
59 : integer, intent(out) :: ierr
60 0 : real(dp) :: dgrad, K_therm, K_T, K_mu, nu, R0, Pr, tau, r_th
61 : include 'formats'
62 0 : dgrad = max(1d-40, grada - gradr) ! positive since Schwarzschild stable
63 0 : K_therm = 4d0*crad*clight*pow3(T)/(3d0*opacity*rho) ! thermal conductivity
64 0 : if (thermohaline_option == 'Kippenhahn') then
65 : ! Kippenhahn, R., Ruschenplatt, G., & Thomas, H.-C. 1980, A&A, 91, 175
66 0 : D_thrm = -3d0*K_therm/(2*rho*cp)*gradL_composition_term/dgrad
67 0 : else if (thermohaline_option == 'Traxler_Garaud_Stellmach_11' .or. &
68 : thermohaline_option == 'Brown_Garaud_Stellmach_13') then
69 0 : call get_diff_coeffs(K_therm, Cp, rho, T, opacity, iso, XH1, K_T, K_mu, nu)
70 0 : R0 = (gradr - grada)/gradL_composition_term
71 0 : Pr = nu/K_T
72 0 : tau = K_mu/K_T
73 0 : r_th = (R0 - 1d0)/(1d0/tau - 1d0)
74 0 : if (r_th >= 1d0) then ! stable if R0 >= 1/tau
75 0 : D_thrm = 0d0
76 0 : else if (Pr < 0d0) then
77 : ! Bad results from get_diff_coeffs will just result in NaNs from thermohaline options, so skip
78 0 : D_thrm = 0d0
79 0 : else if (thermohaline_option == 'Traxler_Garaud_Stellmach_11') then
80 : ! Traxler, Garaud, & Stellmach, ApJ Letters, 728:L29 (2011).
81 : ! also see Denissenkov. ApJ 723:563–579, 2010.
82 0 : D_thrm = 101d0*sqrt(K_mu*nu)*exp(-3.6d0*r_th)*pow(1d0 - r_th, 1.1d0) ! eqn 24
83 : else ! if (s% thermohaline_option == 'Brown_Garaud_Stellmach_13') then
84 0 : D_thrm = K_mu*(Numu(R0, r_th, pr, tau) - 1d0)
85 : end if
86 : else
87 0 : D_thrm = 0
88 0 : ierr = -1
89 0 : write (*, *) 'unknown for MLT thermohaline_option'//trim(thermohaline_option)
90 : end if
91 0 : D_thrm = thermohaline_coeff*D_thrm
92 0 : end subroutine get_D_thermohaline
93 :
94 0 : subroutine get_diff_coeffs(K_therm, Cp, rho, T, opacity, iso, XH1, kt, kmu, vis)
95 : use chem_def, only: chem_isos
96 : real(dp), intent(in) :: K_therm, Cp, rho, T, opacity, XH1
97 : integer, intent(in) :: iso
98 : real(dp), intent(out) :: kt, kmu, vis
99 0 : real(dp) :: loglambdah, loglambdacx, loglambdacy, ccx, ccy, qe4
100 0 : real(dp) :: Bcoeff, chemA, chemZ, acx, acy, nu_mol, nu_rad
101 : real(dp), parameter :: sqrt5 = sqrt(5d0)
102 0 : kt = K_therm/(Cp*rho) ! thermal diffusivity (assumes radiatively dominated)
103 0 : qe4 = pow4(qe)
104 :
105 : ! Log Lambda for pure H (equation 10 from Proffitt Michaud 93)
106 0 : loglambdah = -19.26d0 - 0.5d0*log(rho) + 1.5d0*log(T) - 0.5d0*log(1d0 + 0.5d0*(1 + XH1))
107 0 : nu_rad = 4d0*crad*pow4(T)/(15d0*clight*opacity*pow2(rho)) ! radiative viscosity
108 0 : nu_mol = 0.406d0*sqrt(amu)*pow(boltzm*T, 2.5d0)/(qe4*loglambdah*rho)
109 : ! From Spitzer "Physics of Fully Ionized Gases equation 5-54
110 : ! Assumes pure H. Still trying to work out what it would be for a mixture.
111 0 : vis = nu_mol + nu_rad ! total viscosity
112 :
113 : ! The following is from Proffitt & Michaud, 1993.
114 : ! Their constant B (equation 15)
115 0 : Bcoeff = (15.d0/16.d0)*sqrt(2.d0*amu/(5*pi))*pow(boltzm, 2.5d0)/qe4
116 : ! Extract what species drives the thermohaline concvection
117 0 : chemA = chem_isos%Z_plus_N(iso)
118 0 : chemZ = chem_isos%Z(iso)
119 :
120 0 : if (chemZ > 2) then
121 : ! This is if the driving chemical is NOT He.
122 : ! Log Lambda for H-dominant chem mixture (equation 10)
123 0 : loglambdacx = loglambdah - log(chemz)
124 : ! Log Lambda for He-dominant chem mixture (equation 10)
125 0 : loglambdacy = loglambdah - log(2.d0*chemz)
126 : ! Calculation of C_ij coeffs (equation 12)
127 0 : ccx = log(exp(1.2d0*loglambdacx) + 1.0d0)/1.2d0
128 0 : ccy = log(exp(1.2d0*loglambdacy) + 1.0d0)/1.2d0
129 : ! Reduced masses (I had to guess, from Bahcall & Loeb 1990), with H and He
130 0 : acx = (1.d0*chemA)/(1.d0 + chemA)
131 0 : acy = 4*chemA/(4.d0 + chemA)
132 : ! My formula (see notes) based on Proffitt and Michaud 1993
133 : kmu = 2*Bcoeff*pow(T, 2.5d0)/(sqrt5*rho*chemZ*chemZ)/ &
134 0 : (XH1*sqrt(acx)*ccx + (1 - XH1)*sqrt(acy)*ccy)
135 :
136 : else
137 : ! Log Lambda for H-He mixture (equation 10)
138 : loglambdah = -19.26d0 - log(2d0) - 0.5d0*log(rho) + &
139 0 : 1.5d0*log(T) - 0.5d0*log(1d0 + 0.5d0*(1 + XH1))
140 : ! Calculation of C_ij coeffs (equation 12)
141 0 : ccy = log(exp(1.2d0*loglambdah) + 1d0)/1.2d0
142 : ! My formula (see notes) based on Proffitt and Michaud 1993
143 0 : kmu = (Bcoeff*pow(T, 2.5d0)/(rho*ccy))*(3 + XH1)/((1 + XH1)*(3 + 5*XH1)*(0.7d0 + 0.3d0*XH1))
144 :
145 : end if
146 : ! write(57,*) kt,kmu,vis,chemZ
147 :
148 0 : end subroutine get_diff_coeffs
149 :
150 0 : real(dp) function numu(R0, r_th, prandtl, diffratio)
151 : !Function calculates Nu_mu from input parameters, following Brown et al. 2013.
152 : !Written by P. Garaud (2013). Please email pgaraud@ucsc.edu for troubleshooting.
153 :
154 : real(dp), intent(in) :: R0, r_th, prandtl, diffratio
155 0 : real(dp) :: maxl2, maxl, lambdamax
156 0 : real(dp) :: myvars(2)
157 : integer :: ierr, iter, max_iters
158 :
159 : ! Initialize guess using estimates from Brown et al. 2013
160 0 : call analytical_estimate_th(maxl, lambdamax, r_th, prandtl, diffratio)
161 :
162 0 : myvars(1) = maxl
163 0 : myvars(2) = lambdamax
164 :
165 : !Call Newton relaxation algorithm
166 0 : call NR(myvars, prandtl, diffratio, R0, ierr)
167 :
168 : !If the growth rate is negative, then try another set of parameters as first guess.
169 : !Repeat as many times as necessary until convergence is obtained.
170 0 : iter = 1
171 0 : max_iters = 200
172 0 : do while (iter <= max_iters .and. ((myvars(2) < 0) .or. (ierr /= 0)))
173 : !write(*,*) 'Alternative', r_th,prandtl,diffratio,iter
174 : !Reset guess values
175 0 : myvars(1) = maxl
176 0 : myvars(2) = lambdamax
177 : !Call relaxation for slightly different Pr, tau, R0.
178 0 : call NR(myvars, prandtl*(1d0 + iter*1.d-2), diffratio, R0/(1d0 + iter*1.d-2), ierr)
179 : !If it converged this time, call NR for the real parameters.
180 0 : if (ierr == 0) call NR(myvars, prandtl, diffratio, R0, ierr)
181 : !write(*,*) prandtl,diffratio,R0,myvars(1),myvars(2),ierr
182 : !Otherwise, increase counter and try again.
183 0 : iter = iter + 1
184 : end do
185 :
186 0 : if ((myvars(2) < 0) .or. (ierr /= 0)) then
187 0 : write (*, *) "WARNING: thermohaline Newton relaxation failed to converge, falling back to estimate"
188 0 : maxl2 = maxl*maxl
189 : else ! NR succeeded, so use results in myvars
190 : !Plug solution into "l^2" and lambda.
191 0 : maxl2 = myvars(1)*myvars(1)
192 0 : lambdamax = myvars(2)
193 : !write(*,*) prandtl,diffratio,r_th,maxl2,lambdamax
194 : end if
195 :
196 : !Calculate Nu_mu using Formula (33) from Brown et al, with C = 7.
197 0 : numu = 1.d0 + 49.d0*lambdamax*lambdamax/(diffratio*maxl2*(lambdamax + diffratio*maxl2))
198 :
199 0 : end function numu
200 :
201 0 : subroutine thermohaline_rhs(myx, myf, myj, prandtl, diffratio, R0)
202 : ! This routine is needed for the NR solver.
203 : ! Inputs the two following equations for lambda and maxl2:
204 : ! lambda^3 + a_2 lambda^2 + a_1 lambda + a_0 = 0 (eq. 19 of Brown et al.)
205 : ! b_2 lambda^2 + b_1 lambda + b_0 = 0 (eq. 20 of Brown et al.)
206 : ! Inputs f, the equations, and j, their jacobian.
207 : ! Written by P. Garaud (2013). Please email pgaraud@ucsc.edu for troubleshooting.
208 :
209 : real(dp), intent(in) :: myx(2), prandtl, diffratio, R0
210 : real(dp), intent(out) :: myf(2), myj(2, 2)
211 0 : real(dp) :: a_2, a_1, a_0, b_2, b_1, b_0, myterm, myx1_2, myx1_3, myx1_4
212 :
213 : ! This inputs the coefficients.
214 0 : b_2 = 1d0 + prandtl + diffratio
215 0 : myx1_2 = myx(1)*myx(1)
216 0 : myx1_3 = myx1_2*myx(1)
217 0 : myx1_4 = myx1_3*myx(1)
218 0 : a_2 = myx1_2*b_2
219 0 : myterm = diffratio*prandtl + prandtl + diffratio
220 0 : b_1 = 2*myx1_2*myterm
221 0 : a_1 = myx1_4*myterm + prandtl*(1.0d0 - (1d0/R0))
222 0 : b_0 = 3.d0*myx1_4*diffratio*prandtl + prandtl*(diffratio - (1d0/R0))
223 0 : a_0 = myx1_4*myx1_2*diffratio*prandtl + myx1_2*prandtl*(diffratio - (1d0/R0))
224 :
225 : ! write(*,*) a_2,a_1,a_0,b_2,b_1,b_0
226 :
227 : ! These are equations 19 and 20
228 0 : myf(1) = ((myx(2) + a_2)*myx(2) + a_1)*myx(2) + a_0
229 0 : myf(2) = b_2*myx(2)*myx(2) + b_1*myx(2) + b_0
230 :
231 : ! These are their Jacobians for the NR relaxation.
232 : myj(1, 1) = 2*myx(1)*b_2*myx(2)*myx(2) + &
233 : 4*myx1_3*myterm*myx(2) + 6*myx1_4*myx(1)*diffratio*prandtl &
234 0 : + 2*myx(1)*prandtl*(diffratio - (1d0/R0))
235 0 : myj(1, 2) = 3*myx(2)*myx(2) + 2*a_2*myx(2) + a_1
236 0 : myj(2, 1) = 4*myx(1)*myterm*myx(2) + 12.d0*myx1_3*diffratio*prandtl
237 0 : myj(2, 2) = 2*b_2*myx(2) + b_1
238 :
239 0 : end subroutine thermohaline_rhs
240 :
241 0 : subroutine analytical_estimate_th(maxl, lambdamax, r_th, prandtl, diffratio)
242 : ! Inputs analytical estimates for l and lambda from Brown et al. 2013.
243 :
244 : real(dp), intent(in) :: r_th, prandtl, diffratio
245 : real(dp), intent(out) :: maxl, lambdamax
246 0 : real(dp) :: phi, maxl4, maxl6
247 :
248 0 : phi = diffratio/prandtl
249 :
250 0 : if (r_th < 0.5d0) then
251 0 : if (r_th > prandtl) then
252 0 : maxl = pow((1.d0/(1.d0 + phi)) - 2.d0*dsqrt(r_th*phi)/pow(1d0 + phi, 2.5d0), 0.25d0)
253 : ! Equation (B14)
254 0 : maxl4 = maxl*maxl*maxl*maxl
255 0 : maxl6 = maxl4*maxl*maxl
256 0 : lambdamax = 2*prandtl*phi*maxl6/(1d0 - (1d0 + phi)*maxl4) ! Equation (B11)
257 : else
258 0 : maxl = dsqrt(dsqrt(1d0/(1d0 + phi)) - dsqrt(prandtl)*(1d0 + phi/((1d0 + phi)*(1d0 + phi))))
259 : ! Equation (B5)
260 0 : lambdamax = dsqrt(prandtl) - prandtl*dsqrt(1d0 + phi) ! Equation (B5)
261 : end if
262 : else
263 0 : maxl = pow((one_third)*(1d0 - r_th) + (1d0 - r_th)*(1d0 - r_th)*(5d0 - 4d0*phi)/27d0, 0.25d0)
264 : ! Equation (B19) carried to next order (doesn't work well otherwise)
265 0 : maxl4 = maxl*maxl*maxl*maxl
266 0 : maxl6 = maxl4*maxl*maxl
267 0 : lambdamax = 2d0*prandtl*phi*maxl6/(1d0 - (1d0 + phi)*maxl4) ! Equation (B11)
268 : end if
269 0 : if (lambdamax < 0) then ! shouldn't be needed, but just as precaution
270 0 : maxl = 0.5d0
271 0 : lambdamax = 0.5d0
272 : end if
273 :
274 0 : end subroutine analytical_estimate_th
275 :
276 0 : subroutine NR(xrk, prandtl, diffratio, R0, ierr)
277 : ! Newton Relaxation routine used to solve cubic & quadratic in thermohaline case.
278 : ! Written by P. Garaud (2013). Please email pgaraud@ucsc.edu for troubleshooting.
279 :
280 : real(dp), parameter :: acy = 1.d-13 ! accuracy of NR solution.
281 : integer, parameter :: niter = 20 ! max number of iterations allowed before giving up.
282 : integer, parameter :: & ! array dimension input parameters for dgesvx
283 : n = 2, &
284 : nrhs = 1, &
285 : lda = n, &
286 : ldaf = n, &
287 : ldb = n, &
288 : ldx = n
289 :
290 : real(dp), intent(inout) :: xrk(2)
291 : real(dp), intent(in) :: prandtl, diffratio, R0
292 : integer, intent(out) :: ierr
293 : integer :: iter
294 0 : real(dp) :: f(2) ! Functions f
295 0 : real(dp) :: j(2, 2) ! Jacobian
296 0 : real(dp) :: err, errold ! Error at each iteration
297 0 : real(dp) :: x1_sav, x2_sav
298 0 : real(dp) :: A(lda, n), AF(ldaf, n), R(n), C(n), B(ldb, nrhs), X(ldx, nrhs), &
299 0 : rcond, ferr(nrhs), berr(nrhs), work(4*n)
300 : character :: fact, trans, equed
301 : integer :: ipiv(n), iwork(n)
302 :
303 : include 'formats'
304 :
305 : ! Initialize flags and other counters.
306 0 : ierr = 0
307 0 : iter = 0
308 0 : err = 1d99
309 0 : errold = 1d99
310 : ! Save input guess
311 0 : x1_sav = xrk(1)
312 0 : x2_sav = xrk(2)
313 :
314 : ! While error is too large .and. decreasing, iterate.
315 0 : do while ((err > acy) .and. (ierr == 0) .and. (iter < niter))
316 0 : call thermohaline_rhs(xrk, f, j, prandtl, diffratio, R0)
317 :
318 0 : fact = 'E'
319 0 : trans = 'N'
320 0 : equed = ''
321 :
322 0 : A = j
323 0 : B(1, 1) = f(1)
324 0 : B(2, 1) = f(2)
325 :
326 : call dgesvx(fact, trans, n, nrhs, A, lda, AF, ldaf, ipiv, &
327 : equed, r, c, B, ldb, x, ldx, rcond, ferr, berr, &
328 0 : work, iwork, ierr)
329 :
330 0 : if (ierr /= 0) then
331 : !write(*,*) 'dgesvx failed in thermohaline routine', iter
332 : !write(*,2) j(1,1),j(1,2)
333 : !write(*,2) j(2,1),j(2,2)
334 : else
335 0 : iter = iter + 1
336 0 : f(1) = X(1, 1)
337 0 : f(2) = X(2, 1)
338 0 : err = dsqrt(f(1)*f(1) + f(2)*f(2)) ! Calculate the new error
339 : ! If, after a while, the error is still not decreasing, give up and exit NR.
340 : ! Otherwise, continue.
341 0 : if ((iter > 5) .and. (err > errold)) then
342 : ! Write(*,2) 'Error not decreasing at iter', iter, err, errold
343 0 : ierr = 1
344 : ! Reset xs and exit loop.
345 0 : xrk(1) = x1_sav
346 0 : xrk(2) = x2_sav
347 : else
348 0 : xrk = xrk - f ! The solution is now in f, so update x
349 : errold = err
350 : end if
351 : end if
352 : end do
353 :
354 0 : if (err <= acy) then
355 0 : ierr = 0
356 : else
357 : ! write(*,2) 'NR failed to converge', err, iter
358 0 : ierr = 1
359 : end if
360 :
361 0 : end subroutine NR
362 :
363 : end module thermohaline
|
