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 tdc_support
21 :
22 : use const_def, only: dp, pi, sqrt_2_div_3, boltz_sigma
23 : use num_lib
24 : use utils_lib
25 : use auto_diff
26 :
27 : implicit none
28 :
29 : private
30 : public :: set_Y
31 : public :: Q_bisection_search
32 : public :: dQdZ_bisection_search
33 : public :: Af_bisection_search
34 : public :: convert
35 : public :: unconvert
36 : public :: safe_tanh
37 : public :: tdc_info
38 : public :: eval_Af
39 : public :: eval_xis
40 : public :: compute_Q
41 :
42 : !> Stores the information which is required to evaluate TDC-related quantities and which
43 : !! do not depend on Y.
44 : !!
45 : !! @param report Write debug output if true, not if false.
46 : !! @param mixing_length_alpha Mixing length parameter
47 : !! @param alpha_TDC_DAMP TDC turbulent damping parameter
48 : !! @param alpha_TDC_DAMPR TDC radiative damping parameter
49 : !! @param alpha_TDC_PtdVdt TDC coefficient on P_turb*dV/dt. Physically should probably be 1.
50 : !! @param dt Time-step
51 : !! @param c0 A proportionality factor for the convective luminosity
52 : !! @param L luminosity
53 : !! @param L0 L0 = (Lrad / grad_rad) is the luminosity radiation would carry if dlnT/dlnP = 1.
54 : !! @param A0 Initial convection speed
55 : !! @param T Temperature
56 : !! @param rho Density (g/cm^3)
57 : !! @param dV 1/rho_face - 1/rho_start_face (change in specific volume at the face)
58 : !! @param Cp Heat capacity
59 : !! @param kap Opacity
60 : !! @param Hp Pressure scale height
61 : !! @param gradL gradL is the neutrally buoyant dlnT/dlnP (= grad_ad + grad_mu),
62 : !! @param grada grada is the adiabatic dlnT/dlnP,
63 : !! @param Gamma Gamma is the MLT Gamma efficiency parameter, which we evaluate in steady state from MLT.
64 : type tdc_info
65 : logical :: report
66 : real(dp) :: mixing_length_alpha, alpha_TDC_DAMP, alpha_TDC_DAMPR, alpha_TDC_PtdVdt, dt
67 : type(auto_diff_real_tdc) :: A0, c0, L, L0, gradL, grada
68 : type(auto_diff_real_star_order1) :: T, rho, dV, Cp, kap, Hp, Gamma
69 : end type tdc_info
70 :
71 : contains
72 :
73 : !> Y = +- exp(Z)
74 : !! If Y > 0, Y = exp(Z)
75 : !! If Y < 0, Y = -exp(Z)
76 : !!
77 : !! @param Y_is_positive True if Y > 0, False otherwise.
78 : !! @param Z log|Y|
79 : !! @param Y (output)
80 197007 : type(auto_diff_real_tdc) function set_Y(Y_is_positive, Z) result(Y)
81 : logical, intent(in) :: Y_is_positive
82 : type(auto_diff_real_tdc), intent(in) :: Z
83 197007 : if (Y_is_positive) then
84 195031 : Y = exp(Z)
85 : else
86 1976 : Y = -exp(Z)
87 : end if
88 197007 : end function set_Y
89 :
90 : !> This routine performs a bisection search for Q=0 over a domain in Z for which Q is monotone.
91 : !! Monotonicity is assumed, not verified! This means failures can occur if monotonicity fails.
92 : !! The search continues until the domain is narrowed to less than a width of bracket_tolerance,
93 : !! or until more than max_iter iterations have been taken. Because this is just used to get us in
94 : !! the right ballpark, bracket_tolerance is set quite wide, to 1.
95 : !!
96 : !! There is a check at the start to verify that Q takes on opposite signs on either end of the
97 : !! domain. This is allows us to bail early if there is no root in the domain.
98 : !!
99 : !! @param info tdc_info type storing various quantities that are independent of Y.
100 : !! @param Y_is_positive True if Y > 0, False otherwise.
101 : !! @param lower_bound_Z Lower bound on Z. Note that this is an input *and* an output.
102 : !! @param upper_bound_Z Upper bound on Z. Note that this is an input *and* an output.
103 : !! @param Z Output: the midpoint of the final bounds.
104 : !! @param ierr 0 if everything worked, 1 if the domain does not contain a root.
105 26991 : subroutine Q_bisection_search(info, Y_is_positive, lower_bound_Z, upper_bound_Z, Z, ierr)
106 : ! Inputs
107 : type(tdc_info), intent(in) :: info
108 : logical, intent(in) :: Y_is_positive
109 :
110 : ! Outputs
111 : type(auto_diff_real_tdc), intent(inout) :: lower_bound_Z, upper_bound_Z
112 : type(auto_diff_real_tdc), intent(out) :: Z
113 : integer, intent(out) :: ierr
114 :
115 : ! Parameters
116 : real(dp), parameter :: bracket_tolerance = 1d0
117 : integer, parameter :: max_iter = 30
118 :
119 : ! Intermediates
120 : type(auto_diff_real_tdc) :: Y, Af, Q, Q_ub, Q_lb
121 : integer :: iter
122 :
123 8997 : ierr = 0
124 :
125 8997 : Y = set_Y(Y_is_positive, lower_bound_Z)
126 8997 : call compute_Q(info, Y, Q_lb, Af)
127 :
128 8997 : Y = set_Y(Y_is_positive, upper_bound_Z)
129 8997 : call compute_Q(info, Y, Q_ub, Af)
130 :
131 : ! Check to make sure that the lower and upper bounds on Z actually bracket
132 : ! a solution to Q(Y(Z)) = 0.
133 8997 : if (Q_lb * Q_ub > 0d0) then
134 0 : if (info%report) then
135 0 : write(*,*) 'Q bisection error. Initial Z window does not bracket a solution.'
136 0 : write(*,*) 'Q(Lower Z)',Q_lb%val
137 0 : write(*,*) 'Q(Upper Z)',Q_ub%val
138 0 : write(*,*) 'tolerance', bracket_tolerance
139 0 : write(*,*) 'Y', Y%val
140 0 : write(*,*) 'dYdZ', Y%d1val1
141 0 : write(*,*) 'exp(Z)', exp(Z%val)
142 0 : write(*,*) 'Z', Z%val
143 0 : write(*,*) 'A0', info%A0%val
144 0 : write(*,*) 'c0', info%c0%val
145 0 : write(*,*) 'L', info%L%val
146 0 : write(*,*) 'L0', info%L0%val
147 0 : write(*,*) 'grada', info%grada%val
148 0 : write(*,*) 'gradL', info%gradL%val
149 0 : write(*,'(A)')
150 : end if
151 0 : ierr = 1
152 0 : return
153 : end if
154 :
155 71960 : do iter=1,max_iter
156 71960 : Z = (upper_bound_Z + lower_bound_Z) / 2d0
157 71960 : Y = set_Y(Y_is_positive, Z)
158 :
159 71960 : call compute_Q(info, Y, Q, Af)
160 :
161 71960 : if (Q > 0d0 .and. Q_ub > 0d0) then
162 10 : upper_bound_Z = Z
163 10 : Q_ub = Q
164 71950 : else if (Q > 0d0 .and. Q_lb > 0d0) then
165 38058 : lower_bound_Z = Z
166 38058 : Q_lb = Q
167 33892 : else if (Q < 0d0 .and. Q_ub < 0d0) then
168 33790 : upper_bound_Z = Z
169 33790 : Q_ub = Q
170 102 : else if (Q < 0d0 .and. Q_lb < 0d0) then
171 102 : lower_bound_Z = Z
172 102 : Q_lb = Q
173 : end if
174 :
175 143920 : if (upper_bound_Z - lower_bound_Z < bracket_tolerance) then
176 8997 : Z = (upper_bound_Z + lower_bound_Z) / 2d0
177 8997 : call compute_Q(info, Y, Q, Af)
178 8997 : return
179 : end if
180 : end do
181 :
182 : end subroutine Q_bisection_search
183 :
184 : !> This routine performs a bisection search for dQ/dZ=0 with Y < 0.
185 : !! The domain is assumed to be restricted to have Af > 0, so that dQ/dZ is
186 : !! continuous and monotone. This is checked, and if it fails a MESA ERROR is called
187 : !! rather than throwing an ierr because if the rest of the TDC logic is correct this should
188 : !! not happen.
189 : !!
190 : !! The search continues until the domain is narrowed to less than a width of bracket_tolerance (1d-4),
191 : !! or until more than max_iter iterations have been taken.
192 : !!
193 : !! There is a check at the start to verify that dQ/dZ takes on opposite signs on either end of the
194 : !! domain. This is allows us to bail early if there is no root in the domain.
195 : !!
196 : !! @param info tdc_info type storing various quantities that are independent of Y.
197 : !! @param lower_bound_Z Lower bound on Z. Note that this is an input *and* an output.
198 : !! @param upper_bound_Z Upper bound on Z. Note that this is an input *and* an output.
199 : !! @param Z Output: the midpoint of the final bounds.
200 : !! @param has_root True if a root was found, False otherwise.
201 136 : subroutine dQdZ_bisection_search(info, lower_bound_Z_in, upper_bound_Z_in, Z, has_root)
202 : ! Inputs
203 : type(tdc_info), intent(in) :: info
204 : type(auto_diff_real_tdc), intent(in) :: lower_bound_Z_in, upper_bound_Z_in
205 :
206 : ! Outputs
207 : type(auto_diff_real_tdc), intent(out) :: Z
208 : logical, intent(out) :: has_root
209 :
210 : ! Parameters
211 : real(dp), parameter :: bracket_tolerance = 1d-4
212 : integer, parameter :: max_iter = 50
213 :
214 : ! Intermediates
215 : type(auto_diff_real_tdc) :: lower_bound_Z, upper_bound_Z
216 : type(auto_diff_real_tdc) :: Y, Af, Q, Q_lb, Q_ub, dQdZ, dQdZ_lb, dQdZ_ub
217 : integer :: iter
218 :
219 : ! Set up
220 34 : lower_bound_Z = lower_bound_Z_in !lower_bound_Z_in
221 34 : lower_bound_Z%d1val1 = 1d0
222 34 : upper_bound_Z = upper_bound_Z_in
223 34 : upper_bound_Z%d1val1 = 1d0
224 :
225 : ! Check bounds
226 34 : Y = set_Y(.false., lower_bound_Z)
227 34 : call compute_Q(info, Y, Q_lb, Af)
228 34 : if (Af == 0) then
229 0 : write(*,*) 'Z_lb, A0, Af', lower_bound_Z%val, info%A0%val, Af%val
230 0 : call mesa_error(__FILE__,__LINE__,'bad call to tdc_support dQdZ_bisection_search: Af == 0')
231 : end if
232 34 : dQdZ_lb = differentiate_1(Q_lb)
233 :
234 34 : Y = set_Y(.false., upper_bound_Z)
235 34 : call compute_Q(info, Y, Q_ub, Af)
236 34 : if (Af == 0) then
237 0 : write(*,*) 'Z_ub, A0, Af', lower_bound_Z%val, info%A0%val, Af%val
238 0 : call mesa_error(__FILE__,__LINE__,'bad call to tdc_support dQdZ_bisection_search: Af == 0')
239 : end if
240 34 : dQdZ_ub = differentiate_1(Q_ub)
241 :
242 : ! Check to make sure that the lower and upper bounds on Z actually bracket
243 : ! a solution to dQ/dZ = 0.
244 34 : has_root = .true.
245 34 : if (dQdZ_lb * dQdZ_ub > 0d0) then
246 0 : if (info%report) then
247 0 : write(*,*) 'dQdZ bisection error. Initial Z window does not bracket a solution.'
248 0 : write(*,*) 'Q(Lower Z)',Q_lb%val
249 0 : write(*,*) 'Q(Upper Z)',Q_ub%val
250 0 : write(*,*) 'dQdZ(Lower Z)',dQdZ_lb%val
251 0 : write(*,*) 'dQdZ(Upper Z)',dQdZ_ub%val
252 0 : write(*,*) 'tolerance', bracket_tolerance
253 0 : write(*,*) 'Y', Y%val
254 0 : write(*,*) 'dYdZ', Y%d1val1
255 0 : write(*,*) 'exp(Z)', exp(Z%val)
256 0 : write(*,*) 'Z', Z%val
257 0 : write(*,*) 'A0', info%A0%val
258 0 : write(*,*) 'c0', info%c0%val
259 0 : write(*,*) 'L', info%L%val
260 0 : write(*,*) 'L0', info%L0%val
261 0 : write(*,*) 'grada', info%grada%val
262 0 : write(*,*) 'gradL', info%gradL%val
263 0 : write(*,'(A)')
264 : end if
265 0 : has_root = .false.
266 0 : return
267 : end if
268 :
269 : ! Bisection search
270 680 : do iter=1,max_iter
271 680 : Z = (upper_bound_Z + lower_bound_Z) / 2d0
272 680 : Z%d1val1 = 1d0
273 680 : Y = set_Y(.false., Z)
274 :
275 680 : call compute_Q(info, Y, Q, Af)
276 680 : dQdZ = differentiate_1(Q)
277 :
278 : ! We only ever call this when Y < 0.
279 : ! In this regime, dQ/dZ can take on either sign, and has at most one stationary point.
280 :
281 680 : if (info%report) write(*,*) 'Bisecting dQdZ. Z, dQdZ, Z_lb, dQdZ_lb, Z_ub, dQdZ_ub', &
282 0 : Z%val, dQdZ%val, lower_bound_Z%val, dQdZ_lb%val, upper_bound_Z%val, dQdZ_ub%val
283 :
284 680 : if (dQdZ > 0d0 .and. dQdZ_ub > 0d0) then
285 0 : upper_bound_Z = Z
286 0 : dQdZ_ub = dQdZ
287 680 : else if (dQdZ > 0d0 .and. dQdZ_lb > 0d0) then
288 417 : lower_bound_Z = Z
289 417 : dQdZ_lb = dQdZ
290 263 : else if (dQdZ < 0d0 .and. dQdZ_ub < 0d0) then
291 263 : upper_bound_Z = Z
292 263 : dQdZ_ub = dQdZ
293 0 : else if (dQdZ < 0d0 .and. dQdZ_lb < 0d0) then
294 0 : lower_bound_Z = Z
295 0 : dQdZ_lb = dQdZ
296 : end if
297 :
298 1360 : if (upper_bound_Z - lower_bound_Z < bracket_tolerance) then
299 34 : Z = (upper_bound_Z + lower_bound_Z) / 2d0
300 34 : call compute_Q(info, Y, Q, Af)
301 34 : return
302 : end if
303 : end do
304 :
305 : end subroutine dQdZ_bisection_search
306 :
307 : !> This routine performs a bisection search for the least-negative Y such that Af(Y) = 0.
308 : !! Once we find that, we return a Y that is just slightly less negative so that Af > 0.
309 : !! This is important for our later bisection of dQ/dZ, because we want the upper-Z end of
310 : !! the domain we bisect to have dQ/dZ < 0, which means it has to capture the fact that d(Af)/dZ < 0,
311 : !! and so the Z we return has to be on the positive-Af side of the discontinuity in dQdZ.
312 : !!
313 : !! Note that monotonicity is assumed, not verified!
314 : !! Af(Y) is monotonic in Y for negative Y by construction, so this shouldn't be a problem.
315 : !!
316 : !! The search continues until the domain is narrowed to less than a width of bracket_tolerance (1d-4),
317 : !! or until more than max_iter iterations have been taken.
318 : !!
319 : !! There is a check at the start to verify that Af == 0 at the most-negative end of the domain.
320 : !! This is allows us to bail early if there is no root in the domain.
321 : !!
322 : !! @param info tdc_info type storing various quantities that are independent of Y.
323 : !! @param lower_bound_Z Lower bound on Z. Note that this is an input *and* an output.
324 : !! @param upper_bound_Z Upper bound on Z. Note that this is an input *and* an output.
325 : !! @param Z Output: The bound on the positive-Af side of the root (which always is the lower bound).
326 : !! @param Af Output: Af(Z).
327 : !! @param ierr 0 if everything worked, 1 if the domain does not contain a root.
328 136 : subroutine Af_bisection_search(info, lower_bound_Z_in, upper_bound_Z_in, Z, Af, ierr)
329 : ! Inputs
330 : type(tdc_info), intent(in) :: info
331 : type(auto_diff_real_tdc), intent(in) :: lower_bound_Z_in, upper_bound_Z_in
332 :
333 : ! Outputs
334 : integer, intent(out) :: ierr
335 : type(auto_diff_real_tdc), intent(out) :: Z, Af
336 :
337 : ! Parameters
338 : real(dp), parameter :: bracket_tolerance = 1d-4
339 : integer, parameter :: max_iter = 50
340 :
341 : ! Intermediates
342 : type(auto_diff_real_tdc) :: lower_bound_Z, upper_bound_Z
343 : type(auto_diff_real_tdc) :: Y, Q
344 : integer :: iter
345 :
346 : ! Set up
347 34 : ierr = 0
348 34 : lower_bound_Z = lower_bound_Z_in
349 34 : upper_bound_Z = upper_bound_Z_in
350 :
351 34 : Y = set_Y(.false., upper_bound_Z)
352 34 : call compute_Q(info, Y, Q, Af)
353 34 : if (Af > 0) then ! d(Af)/dZ < 0, so if Af(upper_bound_Z) > 0 there's no solution in this interval.
354 0 : ierr = 1
355 34 : return
356 : end if
357 :
358 34 : Y = set_Y(.false., lower_bound_Z)
359 34 : call compute_Q(info, Y, Q, Af)
360 34 : if (Af == 0) then
361 : ! We want to find Z such that Af(Z) is just barely above zero.
362 : ! We do this by finding Z such that Af(Z) == 0, then backing off to slightly smaller Z.
363 : ! Because d(Af)/dZ < 0, this gives a Z such that Af(Z) > 0.
364 : ! Hence, if Af(lower_bound_Z) == 0 then Af = 0 uniformly in this interval and we cannot
365 : ! return Z such that Af(Z) is just barely non-zero.
366 0 : ierr = 2
367 0 : return
368 : end if
369 :
370 714 : do iter=1,max_iter
371 714 : Z = (upper_bound_Z + lower_bound_Z) / 2d0
372 714 : Y = set_Y(.false., Z)
373 :
374 714 : call compute_Q(info, Y, Q, Af)
375 :
376 : ! Y < 0 so increasing Y means decreasing Z.
377 : ! d(Af)/dY > 0 so d(Af)/dZ < 0.
378 :
379 714 : if (Af > 0d0) then ! Means we are at too-low Z.
380 348 : lower_bound_Z = Z
381 : else
382 366 : upper_bound_Z = Z
383 : end if
384 :
385 1428 : if (upper_bound_Z - lower_bound_Z < bracket_tolerance) then
386 : ! We return the lower bound because this is guaranteed to have Af > 0 (just barely).
387 : ! This is important for our later bisection of dQ/dZ, because we want the upper-Z end of
388 : ! the domain we bisect to have dQ/dZ < 0, which means it has to capture the fact that d(Af)/dZ < 0.
389 34 : Z = lower_bound_Z
390 34 : call compute_Q(info, Y, Q, Af)
391 34 : return
392 : end if
393 : end do
394 :
395 : end subroutine Af_bisection_search
396 :
397 : !> Computes the hyperbolic tangent of x in a way that is numerically safe.
398 : !!
399 : !! @param x Input
400 : !! @param z Output
401 361826 : type(auto_diff_real_tdc) function safe_tanh(x) result(z)
402 : type(auto_diff_real_tdc), intent(in) :: x
403 :
404 361826 : if (x > 50d0) then
405 324742 : z = 1d0
406 37084 : else if (x < -50d0) then
407 0 : z = -1d0
408 : else
409 37084 : z = tanh(x)
410 : end if
411 361826 : end function safe_tanh
412 :
413 : !> The TDC newton solver needs higher-order partial derivatives than
414 : !! the star newton solver, because the TDC one needs to pass back a result
415 : !! which itself contains the derivatives that the star solver needs.
416 : !! These additional derivatives are provided by the auto_diff_real_tdc type.
417 : !!
418 : !! This method converts a auto_diff_real_star_order1 variable into a auto_diff_real_tdc,
419 : !! setting the additional partial derivatives to zero. This 'upgrades' variables storing
420 : !! stellar structure to a form the TDC solver can use.
421 : !!
422 : !! @param K_in, input, an auto_diff_real_star_order1 variable
423 : !! @param K, output, an auto_diff_real_tdc variable.
424 1733824 : type(auto_diff_real_tdc) function convert(K_in) result(K)
425 : type(auto_diff_real_star_order1), intent(in) :: K_in
426 1733824 : K%val = K_in%val
427 58950016 : K%d1Array(1:SIZE(K_in%d1Array)) = K_in%d1Array
428 : K%d1val1 = 0d0
429 58950016 : K%d1val1_d1Array(1:SIZE(K_in%d1Array)) = 0d0
430 1733824 : end function convert
431 :
432 : !> The TDC newton solver needs higher-order partial derivatives than
433 : !! the star newton solver, because the TDC one needs to pass back a result
434 : !! which itself contains the derivatives that the star solver needs.
435 : !! These additional derivatives are provided by the auto_diff_real_tdc type.
436 : !!
437 : !! This method converts a auto_diff_real_tdc variable into a auto_diff_real_star_order1,
438 : !! dropping the additional partial derivatives which (after the TDC solver is done) are
439 : !! no longer needed. This allows the output of the TDC solver to be passed back to the star solver.
440 : !!
441 : !! @param K_in, input, an auto_diff_real_tdc variable
442 : !! @param K, output, an auto_diff_real_star_order1 variable.
443 139481 : type(auto_diff_real_star_order1) function unconvert(K_in) result(K)
444 : type(auto_diff_real_tdc), intent(in) :: K_in
445 139481 : K%val = K_in%val
446 4742354 : K%d1Array = K_in%d1Array(1:SIZE(K%d1Array))
447 139481 : end function unconvert
448 :
449 : !> Q is the residual in the TDC equation, namely:
450 : !!
451 : !! Q = (L - L0 * gradL) - (L0 + c0 * Af) * Y
452 : !!
453 : !! @param info tdc_info type storing various quantities that are independent of Y.
454 : !! @param Y superadiabaticity
455 : !! @param Q The residual of the above equation (an output).
456 : !! @param Af The final convection speed (an output).
457 608852 : subroutine compute_Q(info, Y, Q, Af)
458 : type(tdc_info), intent(in) :: info
459 : type(auto_diff_real_tdc), intent(in) :: Y
460 : type(auto_diff_real_tdc), intent(out) :: Q, Af
461 : type(auto_diff_real_tdc) :: xi0, xi1, xi2, Y_env
462 :
463 : ! Y = grad-gradL
464 : ! Gamma=(grad-gradE)/(gradE-gradL)
465 : ! So
466 : ! Y_env = grad-gradE = (grad-gradL)*Gamma/(1+Gamma) = Y*Gamma/(1+Gamma)
467 : ! So overall we just multiply the Y by Gamma/(1+Gamma) to get Y_env.
468 : !
469 : ! We only use Y_env /= Y when Y > 0 (i.e. the system is convectively unstable)
470 : ! because we only have a Gamma from MLT in that case.
471 : ! so when Y < 0 we just use Y_env = Y.
472 304426 : if (Y > 0) then
473 180913 : Y_env = Y * convert(info%Gamma/(1+info%Gamma))
474 : else
475 123513 : Y_env = Y
476 : end if
477 :
478 : ! Y_env sets the acceleration of blobs.
479 304426 : call eval_xis(info, Y_env, xi0, xi1, xi2)
480 304426 : Af = eval_Af(info%dt, info%A0, xi0, xi1, xi2)
481 :
482 : ! Y_env sets the convective flux but not the radiative flux.
483 304426 : Q = (info%L - info%L0*info%gradL) - info%L0 * Y - info%c0*Af*Y_env
484 :
485 304426 : end subroutine compute_Q
486 :
487 : !! Calculates the coefficients of the TDC velocity equation.
488 : !! The velocity equation is
489 : !!
490 : !! 2 dw/dt = xi0 + w * xi1 + w**2 * xi2 - Lambda
491 : !!
492 : !! where Lambda (currently set to zero) captures coupling between cells.
493 : !!
494 : !! The coefficients xi0/1/2 are given by
495 : !!
496 : !! xi0 = (epsilon_q + S * Y) / w
497 : !! xi1 = (-D_r + p_turb dV/dt) / w^2 [here V = 1/rho is the volume]
498 : !! xi2 = -D / w^3
499 : !!
500 : !! Note that these terms are evaluated on faces, because the convection speed
501 : !! is evaluated on faces. As a result all the inputs must either natively
502 : !! live on faces or be interpolated to faces.
503 : !!
504 : !! This function also has some side effects in terms of storing some of the terms
505 : !! it calculates for plotting purposes.
506 : !!
507 : !! @param info tdc_info type storing various quantities that are independent of Y.
508 : !! @param Y superadiabaticity
509 : !! @param xi0 Output, the constant term in the convective velocity equation.
510 : !! @param xi1 Output, the prefactor of the linear term in the convective velocity equation.
511 : !! @param xi2 Output, the prefactor of the quadratic term in the convective velocity equation.
512 608852 : subroutine eval_xis(info, Y, xi0, xi1, xi2)
513 : ! eval_xis sets up Y with partial wrt Z
514 : ! so results come back with partials wrt Z
515 : type(tdc_info), intent(in) :: info
516 : type(auto_diff_real_tdc), intent(in) :: Y
517 : type(auto_diff_real_tdc), intent(out) :: xi0, xi1, xi2
518 : type(auto_diff_real_tdc) :: S0, D0, DR0
519 : type(auto_diff_real_star_order1) :: gammar_div_alfa, Pt0, dVdt
520 : real(dp), parameter :: x_ALFAS = (1.d0/2.d0)*sqrt_2_div_3
521 : real(dp), parameter :: x_CEDE = (8.d0/3.d0)*sqrt_2_div_3
522 : real(dp), parameter :: x_ALFAP = 2.d0/3.d0
523 : real(dp), parameter :: x_GAMMAR = 2.d0*sqrt(3.d0)
524 :
525 304426 : S0 = convert(x_ALFAS*info%mixing_length_alpha*info%Cp*info%T/info%Hp)*info%grada
526 304426 : S0 = S0*Y
527 304426 : D0 = convert(info%alpha_TDC_DAMP*x_CEDE/(info%mixing_length_alpha*info%Hp))
528 304426 : gammar_div_alfa = info%alpha_TDC_DAMPR*x_GAMMAR/(info%mixing_length_alpha*info%Hp)
529 304426 : DR0 = convert(4d0*boltz_sigma*pow2(gammar_div_alfa)*pow3(info%T)/(pow2(info%rho)*info%Cp*info%kap))
530 304426 : Pt0 = info%alpha_TDC_PtdVdt*x_ALFAP*info%rho
531 304426 : dVdt = info%dV/info%dt
532 :
533 304426 : xi0 = S0
534 304426 : xi1 = -(DR0 + convert(Pt0*dVdt))
535 304426 : xi2 = -D0
536 304426 : end subroutine eval_xis
537 :
538 : !! Calculates the solution to the TDC velocity equation.
539 : !! The velocity equation is
540 : !!
541 : !! 2 dw/dt = xi0 + w * xi1 + w**2 * xi2 - Lambda
542 : !!
543 : !! where Lambda (currently set to zero) captures coupling between cells.
544 : !! The xi0/1/2 variables are constants for purposes of solving this equation.
545 : !!
546 : !! An important related parameter is J:
547 : !!
548 : !! J^2 = xi1^2 - 4 * xi0 * xi2
549 : !!
550 : !! When J^2 > 0 the solution for w is hyperbolic in time.
551 : !! When J^2 < 0 the solution is trigonometric, behaving like tan(J * t / 4 + c).
552 : !!
553 : !! In the trigonometric branch note that once the solution passes through its first
554 : !! root the solution for all later times is w(t) = 0. This is not a solution to the
555 : !! convective velocity equation as written above, but *is* a solution to the convective
556 : !! energy equation (multiply the velocity equation by w), which is the one we actually
557 : !! want to solve (per Radek Smolec's thesis). We just solve the velocity form
558 : !! because it's more convenient.
559 : !!
560 : !! @param dt Time-step
561 : !! @param A0 convection speed from the start of the step (cm/s)
562 : !! @param xi0 The constant term in the convective velocity equation.
563 : !! @param xi1 The prefactor of the linear term in the convective velocity equation.
564 : !! @param xi2 The prefactor of the quadratic term in the convective velocity equation.
565 : !! @param Af Output, the convection speed at the end of the step (cm/s)
566 304426 : function eval_Af(dt, A0, xi0, xi1, xi2) result(Af)
567 : real(dp), intent(in) :: dt
568 : type(auto_diff_real_tdc), intent(in) :: A0, xi0, xi1, xi2
569 : type(auto_diff_real_tdc) :: Af ! output
570 : type(auto_diff_real_tdc) :: J2, J, Jt4, num, den, y_for_atan, root
571 :
572 304426 : J2 = pow2(xi1) - 4d0 * xi0 * xi2
573 :
574 304426 : if (J2 > 0d0) then ! Hyperbolic branch
575 180913 : J = sqrt(abs(J2)) ! Only compute once we know J2 is not 0
576 180913 : Jt4 = 0.25d0 * dt * J
577 180913 : num = safe_tanh(Jt4) * (2d0 * xi0 + A0 * xi1) + A0 * J
578 180913 : den = safe_tanh(Jt4) * (xi1 + 2d0 * A0 * xi2) - J
579 180913 : Af = num / den
580 180913 : if (Af < 0d0) then
581 180913 : Af = -Af
582 : end if
583 123513 : else if (J2 < 0d0) then ! Trigonometric branch
584 58271 : J = sqrt(abs(J2)) ! Only compute once we know J2 is not 0
585 58271 : Jt4 = 0.25d0 * dt * J
586 :
587 : ! This branch contains decaying solutions that reach A = 0, at which point
588 : ! they switch onto the 'zero' branch. So we have to calculate the position of
589 : ! the first root to check it against dt.
590 58271 : y_for_atan = xi1 + 2d0 * A0 * xi2
591 58271 : root = atan(xi1 / J) - atan(y_for_atan / J)
592 :
593 : ! The root enters into a tangent, so we can freely shift it by pi and
594 : ! get another root. We care about the first positive root, and the above prescription
595 : ! is guaranteed to give an answer between (-2*pi,2*pi) because atan produces an answer in [-pi,pi],
596 : ! so we add/subtract a multiple of pi to get the root into [0,pi).
597 58271 : if (root > pi) then
598 0 : root = root - pi
599 58271 : else if (root < -pi) then
600 0 : root = root + 2d0*pi
601 58271 : else if (root < 0d0) then
602 0 : root = root + pi
603 : end if
604 :
605 58271 : if (Jt4 < root) then
606 1576 : num = -xi1 + J * tan(Jt4 + atan(y_for_atan / J))
607 1576 : den = 2d0 * xi2
608 1576 : Af = num / den
609 : else
610 56695 : Af = 0d0
611 : end if
612 : else ! if (J2 == 0d0) then
613 65242 : Af = A0
614 : end if
615 :
616 304426 : end function eval_Af
617 :
618 0 : end module tdc_support
|
