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 :
21 : module tdc
22 :
23 : use const_def, only: dp, sqrt_2_div_3
24 : use num_lib
25 : use utils_lib
26 : use auto_diff
27 : use tdc_support
28 :
29 : implicit none
30 :
31 : private
32 : public :: get_TDC_solution
33 :
34 : contains
35 :
36 : !> Computes the outputs of time-dependent convection theory following the model specified in
37 : !! Radek Smolec's thesis [https://users.camk.edu.pl/smolec/phd_smolec.pdf], which in turn
38 : !! follows the model of Kuhfuss 1986.
39 : !!
40 : !! Internally this solves the equation L = L_conv + L_rad.
41 : !!
42 : !! @param info tdc_info type storing various quantities that are independent of Y.
43 : !! @param scale The scale for computing residuals to the luminosity equation (erg/s).
44 : !! @param Zlb Lower bound on Z.
45 : !! @param Zub Upper bound on Z.
46 : !! @param conv_vel The convection speed (cm/s).
47 : !! @param Y_face The superadiabaticity (dlnT/dlnP - grada, output).
48 : !! @param tdc_num_iters Number of iterations taken in the TDC solver.
49 : !! @param ierr Tracks errors (output).
50 : !! @param Y_face_guess Candidate superadiabaticity for the local solve. Non-positive values disable seeding.
51 195726 : subroutine get_TDC_solution(info, scale, Zlb, Zub, conv_vel, Y_face, tdc_num_iters, Y_face_guess, ierr)
52 : type(tdc_info), intent(in) :: info
53 : real(dp), intent(in) :: scale
54 : type(auto_diff_real_tdc), intent(in) :: Zlb, Zub
55 : type(auto_diff_real_star_order1),intent(out) :: conv_vel, Y_face
56 : integer, intent(out) :: tdc_num_iters, ierr
57 : real(dp), intent(in) :: Y_face_guess
58 :
59 : logical :: Y_is_positive
60 : type(auto_diff_real_tdc) :: Af, Y, Y0, Y1, Z0, Z1, radY
61 : type(auto_diff_real_tdc) :: Q, Q0
62 : logical :: has_root
63 : integer :: iter
64 : include 'formats'
65 :
66 65242 : ierr = 0
67 65242 : if (info%mixing_length_alpha == 0d0 .or. info%dt <= 0d0) then
68 0 : call mesa_error(__FILE__,__LINE__,'bad call to TDC get_TDC_solution')
69 : end if
70 :
71 : ! Determine the sign of the solution.
72 : !
73 : ! If Q(Y=0) is positive then the luminosity is too great to be carried radiatively, so
74 : ! we'll necessarily have Y > 0.
75 : !
76 : ! If Q(Y=0) is negative then the luminosity can be carried by radiation alone, so we'll
77 : ! necessarily have Y < 0.
78 65242 : Y = 0d0
79 65242 : call compute_Q(info, Y, Q0, Af)
80 65242 : if (Q0 > 0d0) then
81 8981 : Y_is_positive = .true.
82 : else
83 56261 : Y_is_positive = .false.
84 : end if
85 :
86 65242 : if (info%report) then
87 0 : open(unit=4,file='out.data')
88 0 : do iter=1,1000
89 0 : Z0 = (Zlb + (Zub-Zlb)*(iter-1)/1000)
90 0 : Y = set_Y(Y_is_positive,Z0)
91 0 : call compute_Q(info, Y, Q, Af)
92 0 : write(4,*) Y%val, Q%val
93 : end do
94 0 : write(*,*) 'Wrote Q(Y) to out.data'
95 : end if
96 :
97 : ! Start down the chain of logic...
98 65242 : if (Y_is_positive) then
99 : ! If Y > 0 then Q(Y) is monotone and we can jump straight to the search.
100 : call bracket_plus_Newton_search(info, scale, Y_is_positive, Zlb, Zub, Y_face, Af, tdc_num_iters, &
101 8981 : Y_face_guess, ierr)
102 8981 : if (ierr /= 0) return
103 8981 : Y = convert(Y_face)
104 8981 : if (info%report) write(*,*) 'Y is positive, Y=',Y_face%val
105 : else
106 56261 : if (info%report) write(*,*) 'Y is negative.'
107 : ! If Y < 0 then Q(Y) is not guaranteed to be monotone, so we have to be more careful.
108 : ! One root we could have is the radiative solution (with Af==0), given by
109 56261 : radY = (info%L - info%L0 * info%gradL) / info%L0
110 :
111 : ! If A0 == 0 then, because Af(Y) is monotone-increasing in Y, we know that for Y < 0, Af(Y) = 0.
112 : ! As a result we can directly write down the solution and get just radY.
113 : ! For numerical reasons we impose a cutoff of 1d-10 cm/s (much smaller than this and the point where
114 : ! A0 reaches zero could be at Z < lower_bound_Z).
115 56261 : if (info%A0 < 1d-10) then
116 56227 : Y = radY
117 56227 : if (info%report) write(*,*) 'A0 == 0, Y=',Y%val
118 : else
119 34 : if (info%report) write(*,*) 'A0 > 0.'
120 : ! Otherwise, we keep going.
121 :
122 : ! We divide the possible functions Af(Z) into three classes:
123 : ! 1. Af(lower_bound_Z) == 0. In this case Af is zero over the whole interval, because d(Af)/dZ <= 0.
124 : ! This means that the system becomes radiative instantly, so we just return the radiative answer.
125 : ! 2. Af(upper_bound_Z) > 0. In this case Af is non-zero over the whole interval, so there is no discontinuity
126 : ! in dQ/dZ, and we can proceed to bisect to search for a root in dQ/dZ.
127 : ! 3. Af(upper_bound_Z) == 0 and Af(lower_bound_Z) > 0. In this case Af hits zero somewhere in our
128 : ! search interval and then sticks there for the rest of the interval. This requires more care to handle,
129 : ! as it produces a discontinuity in dQ/dZ.
130 : !
131 : ! We identify and handle these three cases below.
132 : ! In cases 2 and 3 our goal is to approximate the greatest Z0 such that dQ/dZ is continuous on [Zlb,Z0] and
133 : ! Z0 <= Zub. The continuity requirement is equivalent to Af(Z0) > 0, and the 'approximate the greatest'
134 : ! requirement means that Z0 is near the point where Af(Z) first equals 0.
135 :
136 34 : Y0 = set_Y(.false., Zlb)
137 34 : call compute_Q(info, Y0, Q, Af)
138 34 : if (Af == 0) then
139 : ! Means we're in case 1. Return the radiative Y.
140 0 : Y = radY
141 0 : if (info%report) write(*,*) 'Case 1: Af(Zlb) == 0, Y=radY=',Y%val
142 : else
143 34 : Y0 = set_Y(.false., Zub)
144 34 : call compute_Q(info, Y0, Q, Af)
145 34 : if (Af > 0) then
146 : ! Means we're in case 2. Af > 0 on the whole interval, so dQ/dZ is continuous on the whole interval,
147 : ! so return Z0 = Zub.
148 0 : Z0 = Zub
149 0 : if (info%report) write(*,*) 'Case 1: Af(Zub) > 0, Z0=',Z0%val
150 : else
151 : ! Means we're in case 3.
152 :
153 : ! We now identify the point where Af(Y) = 0.
154 : ! Once we find that, we return a Y that is just slightly less negative so that Af > 0.
155 : ! Call this Y0, corresponding to Z0.
156 : ! This is important for our later bisection of dQ/dZ, because we want the upper-Z end of
157 : ! the domain we bisect to have dQ/dZ < 0, which means it has to capture the fact that d(Af)/dZ < 0,
158 : ! and so the Z we return has to be on the positive-Af side of the discontinuity in dQdZ.
159 34 : call Af_bisection_search(info, Zlb, Zub, Z0, Af, ierr)
160 :
161 : ! ierr /= 0 should be impossible, because we checked the necessary conditions
162 : ! for the bisection search above. Nonetheless, bugs can crop up, so we leave this
163 : ! check in here and leave the checks in Af_bisection_search.
164 34 : if (ierr /= 0) return
165 34 : Y0 = set_Y(.false., Z0)
166 34 : call compute_Q(info, Y0, Q, Af)
167 34 : if (info%report) write(*,*) 'Bisected Af. Y0=',Y0%val,'Af(Y0)=',Af%val
168 : end if
169 :
170 : ! If we're still here it means we were in either case 2 or case 3.
171 : ! In either case, we now need to do a search for where dQ/dZ == 0
172 : ! over the interval [Y0,0] (equivalently from Z=lower_bound to Z=Z0).
173 34 : call dQdZ_bisection_search(info, Zlb, Z0, Z1, has_root)
174 34 : if (has_root) then
175 34 : Y1 = set_Y(.false., Z1)
176 34 : if (info%report) write(*,*) 'Bisected dQdZ, found root, ',Y1%val
177 34 : call compute_Q(info, Y1, Q, Af)
178 34 : if (Q < 0) then ! Means there are no roots with Af > 0.
179 18 : if (info%report) write(*,*) 'Root has Q<0, Q=',Q%val,'Y=',radY%val
180 18 : Y = radY
181 : else
182 16 : if (info%report) write(*,*) 'Root has Q>0. Q(Y1)=',Q%val
183 : ! Do a search over [lower_bound, Z1]. If we find a root, that's the root closest to zero so call it done.
184 16 : if (info%report) write(*,*) 'Searching from Y=',-exp(Zlb%val),'to Y=',-exp(Z1%val)
185 : call bracket_plus_Newton_search(info, scale, Y_is_positive, Zlb, Z1, Y_face, Af, tdc_num_iters, &
186 16 : 0d0, ierr)
187 16 : Y = convert(Y_face)
188 16 : if (info%report) write(*,*) 'ierr',ierr, tdc_num_iters
189 16 : if (ierr /= 0) then
190 0 : if (info%report) write(*,*) 'No root found. Searching from Y=',-exp(Z1%val),'to Y=',-exp(Z0%val)
191 : ! Do a search over [Z1, Z0]. If we find a root, that's the root closest to zero so call it done.
192 : ! Note that if we get to this stage there is (mathematically) guaranteed to be a root, modulo precision issues.
193 : call bracket_plus_Newton_search(info, scale, Y_is_positive, Z1, Z0, Y_face, Af, tdc_num_iters, &
194 0 : 0d0, ierr)
195 0 : Y = convert(Y_face)
196 : end if
197 16 : if (info%report) write(*,*) 'Y=',Y%val
198 : end if
199 : else
200 0 : if (info%report) write(*,*) 'Bisected dQdZ, no root found.'
201 0 : call compute_Q(info, Y0, Q, Af)
202 0 : if (Q > 0) then ! Means there's a root in [Y0,0] so we bracket search from [lower_bound,Z0]
203 : call bracket_plus_Newton_search(info, scale, Y_is_positive, Zlb, Z0, Y_face, Af, tdc_num_iters, &
204 0 : 0d0, ierr)
205 0 : Y = convert(Y_face)
206 0 : if (info%report) write(*,*) 'Q(Y0) > 0, bisected and found Y=',Y%val
207 : else ! Means there's no root in [Y0,0] so the only root is radY.
208 0 : if (info%report) write(*,*) 'Q(Y0) < 0, Y=',radY%val
209 0 : Y = radY
210 : end if
211 : end if
212 : end if
213 : end if
214 : end if
215 :
216 :
217 : ! Process Y into the various outputs.
218 65242 : call compute_Q(info, Y, Q, Af)
219 65242 : Y_face = unconvert(Y)
220 65242 : conv_vel = sqrt_2_div_3*unconvert(Af)
221 :
222 : end subroutine get_TDC_solution
223 :
224 : !> Performs a bracket search for the solution to Q=0 over
225 : !! a domain in which Q is guaranteed to be monotone. Then
226 : !! refines the result with a Newton solver to both get a better
227 : !! solution and endow the solution with partial derivatives.
228 : !!
229 : !! @param info tdc_info type storing various quantities that are independent of Y.
230 : !! @param scale The scale for computing residuals to the luminosity equation (erg/s).
231 : !! @param Zlb Lower bound on Z.
232 : !! @param Zub Upper bound on Z.
233 : !! @param Y_face The superadiabaticity (dlnT/dlnP - grada, output).
234 : !! @param tdc_num_iters Number of iterations taken in the TDC solver.
235 : !! @param ierr Tracks errors (output).
236 : !! @param Y_face_guess Candidate superadiabaticity for the local solve. Non-positive values disable seeding.
237 26991 : subroutine bracket_plus_Newton_search(info, scale, Y_is_positive, Zlb, Zub, Y_face, Af, tdc_num_iters, &
238 : Y_face_guess, ierr)
239 : type(tdc_info), intent(in) :: info
240 : logical, intent(in) :: Y_is_positive
241 : real(dp), intent(in) :: scale
242 : type(auto_diff_real_tdc), intent(in) :: Zlb, Zub
243 : type(auto_diff_real_star_order1),intent(out) :: Y_face
244 : type(auto_diff_real_tdc), intent(out) :: Af
245 : integer, intent(out) :: tdc_num_iters
246 : integer, intent(out) :: ierr
247 : real(dp), intent(in) :: Y_face_guess
248 :
249 : type(auto_diff_real_tdc) :: Y, Z, Q, Qc, Z_new, correction, lower_bound_Z, upper_bound_Z
250 : type(auto_diff_real_tdc) :: dQdZ
251 : integer :: iter, line_iter
252 : logical :: converged, have_derivatives, corr_has_derivatives
253 : real(dp), parameter :: correction_tolerance = 1d-13
254 : real(dp), parameter :: residual_tolerance = 1d-8
255 : integer, parameter :: max_iter = 200
256 : integer, parameter :: max_line_search_iter = 5
257 : include 'formats'
258 :
259 8997 : ierr = 0
260 :
261 : ! We start by bisecting to find a narrow interval around the root.
262 8997 : lower_bound_Z = Zlb
263 8997 : upper_bound_Z = Zub
264 :
265 8997 : call try_seeded_positive_Y_bracket(info, Y_is_positive, Zlb, Zub, lower_bound_Z, upper_bound_Z, Y_face_guess)
266 :
267 : ! Perform bisection search.
268 8997 : call Q_bisection_search(info, Y_is_positive, lower_bound_Z, upper_bound_Z, Z, ierr)
269 8997 : if (ierr /= 0) return
270 :
271 : ! Set up Z from bisection search
272 8997 : Z%d1val1 = 1d0 ! Set derivative dZ/dZ=1 for Newton iterations.
273 8997 : if (info%report) write(*,*) 'Z from bisection search', Z%val
274 8997 : if (info%report) write(*,*) 'lower_bound_Z, upper_bound_Z',lower_bound_Z%val,upper_bound_Z%val
275 :
276 : ! Now we refine the solution with a Newton solve.
277 : ! This also let's us pick up the derivative of the solution with respect to input parameters.
278 :
279 : ! Initialize starting values for TDC Newton iterations.
280 8997 : dQdz = 0d0
281 8997 : converged = .false.
282 8997 : have_derivatives = .false. ! Tracks if we've done at least one Newton iteration.
283 : ! Need to do this before returning to endow Y with partials
284 : ! with respect to the structure variables.
285 41127 : do iter = 1, max_iter
286 41127 : Y = set_Y(Y_is_positive, Z)
287 41127 : call compute_Q(info, Y, Q, Af)
288 41127 : if (is_bad(Q%val)) then
289 0 : ierr = 1
290 0 : exit
291 : end if
292 :
293 41127 : if (abs(Q%val)/scale <= residual_tolerance .and. have_derivatives) then
294 : ! Can't exit on the first iteration, otherwise we have no derivative information.
295 8997 : if (info%report) write(*,2) 'converged', iter, abs(Q%val)/scale, residual_tolerance
296 : converged = .true.
297 : exit
298 : end if
299 :
300 : ! We use the fact that Q(Y) is monotonic to iteratively refined bounds on Q.
301 32130 : dQdZ = differentiate_1(Q)
302 32130 : if (Q > 0d0 .and. dQdZ < 0d0) then
303 5838 : lower_bound_Z = Z
304 26292 : else if (Q > 0d0 .and. dQdZ > 0d0) then
305 29 : upper_bound_Z = Z
306 26263 : else if (Q < 0d0 .and. dQdZ < 0d0) then
307 26242 : upper_bound_Z = Z
308 : else
309 21 : lower_bound_Z = Z
310 : end if
311 :
312 32130 : if (is_bad(dQdZ%val) .or. abs(dQdZ%val) < 1d-99) then
313 0 : ierr = 1
314 0 : exit
315 : end if
316 :
317 32130 : correction = -Q/dQdz
318 32130 : corr_has_derivatives = .true.
319 :
320 : ! Do a line search to avoid steps that are too big.
321 32130 : do line_iter=1,max_line_search_iter
322 :
323 32130 : if (abs(correction) < correction_tolerance .and. have_derivatives) then
324 : ! Can't get much more precision than this.
325 0 : converged = .true.
326 0 : exit
327 : end if
328 :
329 32130 : Z_new = Z + correction
330 : if (corr_has_derivatives) then
331 : have_derivatives = .true.
332 : end if
333 :
334 : ! If the correction pushes the solution out of bounds then we know
335 : ! that was a bad step. Bad steps are still in the same direction, they just
336 : ! go too far, so we replace that result with one that's halfway to the relevant bound.
337 32130 : if (Z_new > upper_bound_Z) then
338 980 : Z_new = (Z + upper_bound_Z) / 2d0
339 31150 : else if (Z_new < lower_bound_Z) then
340 0 : Z_new = (Z + lower_bound_Z) / 2d0
341 : end if
342 :
343 32130 : Y = set_Y(Y_is_positive,Z_new)
344 :
345 32130 : call compute_Q(info, Y, Qc, Af)
346 :
347 32130 : if (is_bad(Qc%val)) then
348 0 : correction = 0.5d0 * correction
349 32130 : else if (abs(Qc) < abs(Q)) then
350 : exit
351 : else
352 0 : correction = 0.5d0 * correction
353 : end if
354 : end do
355 :
356 32130 : if (info%report) write(*,3) 'i, li, Z_new, Z, low_bnd, upr_bnd, Q, dQdZ, corr', iter, line_iter, &
357 0 : Z_new%val, Z%val, lower_bound_Z%val, upper_bound_Z%val, Q%val, dQdZ%val, correction%val
358 32130 : Z_new%d1val1 = 1d0 ! Ensures that dZ/dZ = 1.
359 32130 : Z = Z_new
360 :
361 32130 : Y = set_Y(Y_is_positive,Z)
362 105387 : if (converged) exit
363 :
364 : end do
365 :
366 0 : if (.not. converged) then
367 0 : ierr = 1
368 0 : if (info%report) then
369 0 : !$OMP critical (tdc_crit0)
370 0 : write(*,*) 'failed get_TDC_solution TDC_iter', &
371 0 : iter
372 0 : write(*,*) 'scale', scale
373 0 : write(*,*) 'Q/scale', Q%val/scale
374 0 : write(*,*) 'tolerance', residual_tolerance
375 0 : write(*,*) 'dQdZ', dQdZ%val
376 0 : write(*,*) 'Y', Y%val
377 0 : write(*,*) 'dYdZ', Y%d1val1
378 0 : write(*,*) 'exp(Z)', exp(Z%val)
379 0 : write(*,*) 'Z', Z%val
380 0 : write(*,*) 'Af', Af%val
381 0 : write(*,*) 'A0', info%A0%val
382 0 : write(*,*) 'c0', info%c0%val
383 0 : write(*,*) 'L', info%L%val
384 0 : write(*,*) 'L0', info%L0%val
385 0 : write(*,*) 'grada', info%grada%val
386 0 : write(*,*) 'gradL', info%gradL%val
387 0 : write(*,'(A)')
388 : !$OMP end critical (tdc_crit0)
389 : end if
390 0 : return
391 : end if
392 :
393 : ! Unpack output
394 8997 : Y_face = unconvert(Y)
395 8997 : tdc_num_iters = iter
396 : end subroutine bracket_plus_Newton_search
397 :
398 :
399 8997 : subroutine try_seeded_positive_Y_bracket(info, Y_is_positive, Zlb, Zub, lower_bound_Z, upper_bound_Z, Y_face_guess)
400 : type(tdc_info), intent(in) :: info
401 : logical, intent(in) :: Y_is_positive
402 : type(auto_diff_real_tdc), intent(in) :: Zlb, Zub
403 : type(auto_diff_real_tdc), intent(inout) :: lower_bound_Z, upper_bound_Z
404 : real(dp), intent(in) :: Y_face_guess
405 :
406 : type(auto_diff_real_tdc) :: Af, Q_lb, Q_ub, Y, Z_seed, Z_seed_lb, Z_seed_ub
407 : ! Try a local bracket within a factor exp(seed_bracket_half_width) of the seed.
408 : real(dp), parameter :: seed_bracket_half_width = 3d0
409 :
410 17978 : if (.not. Y_is_positive) return
411 8981 : if (is_bad(Y_face_guess) .or. Y_face_guess <= 0d0) return
412 :
413 0 : Z_seed = log(Y_face_guess)
414 0 : if (is_bad(Z_seed%val)) return
415 0 : if (Z_seed%val < Zlb%val .or. Z_seed%val > Zub%val) return
416 :
417 0 : Z_seed_lb = max(Zlb%val, Z_seed%val - seed_bracket_half_width)
418 0 : Z_seed_ub = min(Zub%val, Z_seed%val + seed_bracket_half_width)
419 :
420 0 : Y = set_Y(Y_is_positive, Z_seed_lb)
421 0 : call compute_Q(info, Y, Q_lb, Af)
422 0 : Y = set_Y(Y_is_positive, Z_seed_ub)
423 0 : call compute_Q(info, Y, Q_ub, Af)
424 :
425 0 : if (is_bad(Q_lb%val) .or. is_bad(Q_ub%val)) return
426 0 : if (Q_lb * Q_ub > 0d0) return
427 :
428 0 : lower_bound_Z = Z_seed_lb
429 0 : upper_bound_Z = Z_seed_ub
430 : end subroutine try_seeded_positive_Y_bracket
431 :
432 : end module tdc
|
