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