Line data Source code
1 : ! ***********************************************************************
2 : !
3 : ! Copyright (C) 2010-2025 Ebraheem Farag & 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_hydro_support
21 :
22 : use star_private_def
23 : use const_def, only: dp, ln10, pi, lsun, rsun
24 : use utils_lib, only: is_bad
25 : use auto_diff
26 : use auto_diff_support
27 : use star_utils
28 :
29 : implicit none
30 :
31 : private
32 : public :: remesh_for_TDC_pulsations
33 :
34 : contains
35 :
36 0 : subroutine remesh_for_TDC_pulsations(s, ierr)
37 : ! uses these controls
38 : ! TDC_hydro_nz = 190
39 : ! TDC_hydro_nz_outer = 40
40 : ! TDC_hydro_nz_inner = 20
41 : ! TDC_hydro_nz_T_gradient = 20
42 : ! TDC_hydro_T_anchor = 11d3
43 : ! TDC_hydro_dq_1_factor = 2d0
44 : use interp_1d_def, only: pm_work_size
45 : use interp_1d_lib, only: interpolate_vector_pm
46 : use hydro_vars, only: set_cgrav
47 : type(star_info), pointer :: s
48 : integer, intent(out) :: ierr
49 : integer :: k, j, nz_old, nz
50 : real(dp) :: xm_anchor, P_surf, T_surf, old_L1, old_r1, old_J, old_abs_J
51 : real(dp), allocatable, dimension(:) :: &
52 0 : xm_old, xm, xm_mid_old, xm_mid, v_old, v_new
53 0 : real(dp), pointer :: work1(:) ! =(nz_old+1, pm_work_size)
54 : include 'formats'
55 0 : ierr = 0
56 0 : nz_old = s%nz
57 0 : nz = s%TDC_hydro_nz
58 0 : call validate_controls2
59 0 : if (ierr /= 0) return
60 0 : call setvars2(ierr)
61 0 : if (ierr /= 0) call mesa_error(__FILE__, __LINE__, 'remesh_for_TDC failed in setvars')
62 0 : old_J = 0d0
63 0 : old_abs_J = 0d0
64 0 : if (s%rotation_flag) then
65 0 : old_J = dot_product(s%dm_bar(1:nz_old), s%j_rot(1:nz_old))
66 0 : old_abs_J = dot_product(s%dm_bar(1:nz_old), abs(s%j_rot(1:nz_old)))
67 : end if
68 0 : old_L1 = s%L(1)
69 0 : old_r1 = s%r(1)
70 0 : call set_phot_info(s) ! sets Teff
71 0 : call get_PT_surf2(P_surf, T_surf, ierr)
72 0 : if (ierr /= 0) call mesa_error(__FILE__, __LINE__, 'remesh_for_TDC failed in get_PT_surf')
73 : allocate ( &
74 0 : xm_old(nz_old + 1), xm_mid_old(nz_old), v_old(nz_old + 1), &
75 0 : xm(nz + 1), xm_mid(nz), v_new(nz + 1), work1((nz_old + 1)*pm_work_size))
76 0 : call set_xm_old2
77 0 : call find_xm_anchor2
78 0 : if (ierr /= 0) then
79 0 : deallocate(work1)
80 0 : return
81 : end if
82 0 : call set_xm_new2
83 0 : if (ierr /= 0) then
84 0 : deallocate(work1)
85 0 : return
86 : end if
87 0 : call interpolate1_face_val2(s%i_lnR, log(max(1d0, s%r_center)))
88 0 : call check_new_lnR2
89 0 : call interpolate1_face_val2(s%i_lum, s%L_center)
90 0 : if (s%i_v /= 0) call interpolate1_face_val2(s%i_v, s%v_center)
91 0 : call set_new_lnd2
92 0 : call interpolate1_cell_val2(s%i_lnT)
93 0 : if (s%i_u /= 0) call interpolate1_cell_val2(s%i_u)
94 0 : do j = 1, s%species
95 0 : call remap1_xa2(j)
96 : end do
97 0 : if (s%rotation_flag) call remap_rotation2(old_J, old_abs_J)
98 0 : s%nz = nz
99 0 : call update_composition_info2
100 0 : call set_cgrav(s, ierr)
101 0 : if (ierr /= 0) call mesa_error(__FILE__, __LINE__, 'remesh_for_TDC failed in set_cgrav')
102 0 : call revise_lnT_for_QHSE2(P_surf, ierr)
103 0 : if (ierr /= 0) call mesa_error(__FILE__, __LINE__, 'remesh_for_TDC failed in revise_lnT_for_QHSE')
104 0 : if (s%rotation_flag) call set_rotation_seed2
105 0 : deallocate (work1)
106 0 : write (*, 1) 'new old L_surf/Lsun', s%xh(s%i_lum, 1)/Lsun, old_L1/Lsun
107 0 : write (*, 1) 'new old R_surf/Rsun', exp(s%xh(s%i_lnR, 1))/Rsun, old_r1/Rsun
108 0 : write (*, '(A)')
109 :
110 : contains
111 :
112 0 : subroutine validate_controls2
113 : integer :: nz_base
114 : include 'formats'
115 0 : nz_base = nz - s%TDC_hydro_nz_T_gradient
116 0 : if (s%RSP_flag .or. s%RSP2_flag) then
117 0 : write(*,'(A)') 'TDC remesh cannot be applied after enabling RSP or RSP2'
118 0 : ierr = -1
119 0 : else if (nz > nz_old) then
120 0 : write(*,3) 'TDC remesh cannot increase the number of zones', nz, nz_old
121 0 : ierr = -1
122 0 : else if (s%TDC_hydro_nz_T_gradient < 0) then
123 0 : write(*,1) 'TDC_hydro_nz_T_gradient must be nonnegative', &
124 0 : real(s%TDC_hydro_nz_T_gradient, dp)
125 0 : ierr = -1
126 0 : else if (s%TDC_hydro_nz_inner < 0) then
127 0 : write(*,1) 'TDC_hydro_nz_inner must be nonnegative', &
128 0 : real(s%TDC_hydro_nz_inner, dp)
129 0 : ierr = -1
130 0 : else if (s%TDC_hydro_nz_outer < 2) then
131 0 : write(*,1) 'TDC_hydro_nz_outer must be at least 2', &
132 0 : real(s%TDC_hydro_nz_outer, dp)
133 0 : ierr = -1
134 0 : else if (s%remesh_for_TDC_pulsations_log_core_zoning .and. &
135 : nz_base - s%TDC_hydro_nz_outer < 2) then
136 0 : write(*,3) 'TDC remesh needs at least two interior zones', &
137 0 : nz_base, s%TDC_hydro_nz_outer
138 0 : ierr = -1
139 0 : else if (.not. s%remesh_for_TDC_pulsations_log_core_zoning .and. &
140 : nz_base - s%TDC_hydro_nz_outer - s%TDC_hydro_nz_inner < 2) then
141 0 : write(*,3) 'TDC remesh needs at least two middle zones', &
142 0 : nz_base, s%TDC_hydro_nz_outer, s%TDC_hydro_nz_inner
143 0 : ierr = -1
144 : else if (.not. s%remesh_for_TDC_pulsations_log_core_zoning .and. &
145 0 : s%TDC_hydro_nz_inner > 0 .and. s%max_center_cell_dq <= 0d0) then
146 0 : write(*,1) 'max_center_cell_dq must be positive for inner TDC zoning', &
147 0 : s%max_center_cell_dq
148 0 : ierr = -1
149 0 : else if (s%TDC_hydro_dq_1_factor <= 0d0) then
150 0 : write(*,1) 'TDC_hydro_dq_1_factor must be positive', &
151 0 : s%TDC_hydro_dq_1_factor
152 0 : ierr = -1
153 0 : else if (s%TDC_hydro_T_anchor <= 0d0) then
154 0 : write(*,1) 'TDC_hydro_T_anchor must be positive', s%TDC_hydro_T_anchor
155 0 : ierr = -1
156 : end if
157 0 : end subroutine validate_controls2
158 :
159 0 : subroutine setvars2(ierr)
160 : use hydro_vars, only: unpack_xh, set_hydro_vars
161 : integer, intent(out) :: ierr
162 : logical, parameter :: &
163 : skip_basic_vars = .false., &
164 : skip_micro_vars = .false., &
165 : skip_m_grav_and_grav = .false., &
166 : skip_net = .true., &
167 : skip_neu = .true., &
168 : skip_kap = .false., &
169 : skip_grads = .true., &
170 : skip_rotation = .true., &
171 : skip_brunt = .true., &
172 : skip_other_cgrav = .true., &
173 : skip_mixing_info = .true., &
174 : skip_set_cz_bdy_mass = .true., &
175 : skip_mlt = .true., &
176 : skip_eos = .false.
177 : ierr = 0
178 0 : call unpack_xh(s, ierr)
179 0 : if (ierr /= 0) call mesa_error(__FILE__, __LINE__, 'remesh_for_TDC failed in unpack_xh')
180 : call set_hydro_vars( &
181 : s, 1, nz_old, skip_basic_vars, &
182 : skip_micro_vars, skip_m_grav_and_grav, skip_eos, skip_net, skip_neu, &
183 : skip_kap, skip_grads, skip_rotation, skip_brunt, skip_other_cgrav, &
184 0 : skip_mixing_info, skip_set_cz_bdy_mass, skip_mlt, ierr)
185 0 : if (ierr /= 0) call mesa_error(__FILE__, __LINE__, 'remesh_for_TDC failed in set_hydro_vars')
186 0 : end subroutine setvars2
187 :
188 0 : subroutine get_PT_surf2(P_surf, T_surf, ierr)
189 : use atm_support, only: get_atm_PT
190 : real(dp), intent(out) :: P_surf, T_surf
191 : integer, intent(out) :: ierr
192 : real(dp) :: &
193 : Teff, lnT_surf, dlnT_dL, dlnT_dlnR, dlnT_dlnM, dlnT_dlnkap, &
194 : lnP_surf, dlnP_dL, dlnP_dlnR, dlnP_dlnM, dlnP_dlnkap
195 : logical, parameter :: skip_partials = .true.
196 0 : ierr = 0
197 0 : call set_phot_info(s) ! sets s% Teff
198 0 : Teff = s%Teff
199 : call get_atm_PT( & ! this uses s% opacity(1)
200 : s, s%tau_factor*s%tau_base, s%L(1), s%r(1), s%m(1), s%cgrav(1), skip_partials, &
201 : Teff, lnT_surf, dlnT_dL, dlnT_dlnR, dlnT_dlnM, dlnT_dlnkap, &
202 0 : lnP_surf, dlnP_dL, dlnP_dlnR, dlnP_dlnM, dlnP_dlnkap, ierr)
203 0 : if (ierr /= 0) call mesa_error(__FILE__, __LINE__, 'get_P_surf failed in get_atm_PT')
204 0 : P_surf = exp(lnP_surf)
205 0 : T_surf = exp(lnT_surf)
206 0 : end subroutine get_PT_surf2
207 :
208 0 : subroutine set_xm_old2
209 0 : xm_old(1) = 0d0
210 0 : do k = 2, nz_old
211 0 : xm_old(k) = xm_old(k - 1) + s%dm(k - 1)
212 : end do
213 0 : xm_old(nz_old + 1) = s%xmstar
214 0 : do k = 1, nz_old
215 0 : xm_mid_old(k) = xm_old(k) + 0.5d0*s%dm(k)
216 : end do
217 0 : end subroutine set_xm_old2
218 :
219 0 : subroutine find_xm_anchor2
220 : real(dp) :: lnT_anchor, xmm1, xm00, lnTm1, lnT00
221 : include 'formats'
222 0 : lnT_anchor = log(s%TDC_hydro_T_anchor)
223 0 : if (lnT_anchor <= s%xh(s%i_lnT, 1)) then
224 0 : write (*, 1) 'T_anchor < T_surf', s%TDC_hydro_T_anchor, exp(s%xh(s%i_lnT, 1))
225 0 : call mesa_error(__FILE__, __LINE__, 'find_xm_anchor')
226 : end if
227 0 : xm_anchor = xm_old(nz_old)
228 0 : do k = 2, nz_old
229 0 : if (s%xh(s%i_lnT, k) >= lnT_anchor) then
230 0 : xmm1 = xm_mid_old(k - 1)
231 0 : xm00 = xm_mid_old(k)
232 0 : lnTm1 = s%xh(s%i_lnT, k - 1)
233 0 : lnT00 = s%xh(s%i_lnT, k)
234 : xm_anchor = xmm1 + &
235 0 : (xm00 - xmm1)*(lnT_anchor - lnTm1)/(lnT00 - lnTm1)
236 0 : if (is_bad(xm_anchor) .or. xm_anchor <= 0d0) then
237 0 : write (*, 2) 'bad xm_anchor', k, xm_anchor, xmm1, xm00, lnTm1, lnT00, lnT_anchor, s%lnT(1)
238 0 : call mesa_error(__FILE__, __LINE__, 'find_xm_anchor')
239 : end if
240 0 : return
241 : end if
242 : end do
243 0 : write(*,1) 'T_anchor exceeds the model temperature', s%TDC_hydro_T_anchor
244 0 : ierr = -1
245 : end subroutine find_xm_anchor2
246 :
247 0 : subroutine set_xm_new2 ! sets xm, dm, m, dq, q
248 : integer :: nz_outer, nz_inner, nz_T_gradient, nz_base, n_middle, k
249 : real(dp) :: dq_1_factor, dxm_outer, lnx, dlnx, base_dm, &
250 : rem_mass, center_dm, peak_dm, correction, ratio, H, H_inner
251 : real(dp) :: H_low, H_high, H_mid, f_low, f_high, f_mid
252 : integer :: iter
253 : include 'formats'
254 0 : nz_outer = s%TDC_hydro_nz_outer
255 0 : nz_inner = s%TDC_hydro_nz_inner
256 0 : nz_T_gradient = s%TDC_hydro_nz_T_gradient
257 0 : nz_base = nz - nz_T_gradient
258 0 : dq_1_factor = s%TDC_hydro_dq_1_factor
259 0 : ratio = max(2.5d0, s%mesh_max_allowed_ratio)
260 0 : dxm_outer = xm_anchor/(nz_outer - 1d0 + dq_1_factor)
261 0 : xm(1) = 0d0
262 0 : xm(2) = dxm_outer*dq_1_factor
263 0 : s%dm(1) = xm(2)
264 0 : do k = 3, nz_outer + 1
265 0 : xm(k) = xm(k - 1) + dxm_outer
266 0 : s%dm(k - 1) = dxm_outer
267 : end do
268 :
269 0 : if (.not. s%remesh_for_TDC_pulsations_log_core_zoning) then
270 : ! do rsp style core zoning with a power law on dq
271 0 : n_middle = nz_base - nz_outer - nz_inner
272 0 : rem_mass = s%xmstar - xm(nz_outer + 1)
273 0 : base_dm = dxm_outer
274 :
275 0 : if (nz_inner == 0) then
276 : ! Original single inward-increasing power law.
277 0 : H_low = 1.001d0
278 0 : H_high = 1.40d0
279 0 : f_low = base_dm*(1d0 - pow(H_low, real(n_middle, dp)))/(1d0 - H_low) - rem_mass
280 0 : f_high = base_dm*(1d0 - pow(H_high, real(n_middle, dp)))/(1d0 - H_high) - rem_mass
281 0 : if (f_low*f_high > 0d0) then
282 0 : write(*,2) 'failed to bracket TDC core zoning ramp', &
283 0 : n_middle, f_low, f_high
284 0 : ierr = -1
285 0 : return
286 : end if
287 0 : do iter = 1, 1000
288 0 : H_mid = 0.5d0*(H_low + H_high)
289 0 : f_mid = base_dm*(1d0 - pow(H_mid, real(n_middle, dp)))/(1d0 - H_mid) - rem_mass
290 0 : if (abs(f_mid) < 1d-12*rem_mass) exit
291 0 : if (f_low*f_mid <= 0d0) then
292 0 : H_high = H_mid
293 0 : f_high = f_mid
294 : else
295 0 : H_low = H_mid
296 0 : f_low = f_mid
297 : end if
298 : end do
299 0 : H = H_mid
300 :
301 0 : s%dm(nz_outer + 1) = base_dm
302 0 : do k = nz_outer + 2, nz_base - 1
303 0 : s%dm(k) = H*s%dm(k - 1)
304 : end do
305 0 : s%dm(nz_base) = s%xmstar - sum(s%dm(1:nz_base - 1))
306 :
307 : else
308 : ! Add an inward-decreasing ramp that ends at max_center_cell_dq.
309 : H_low = 1d0
310 0 : H_high = ratio
311 0 : center_dm = s%max_center_cell_dq*s%xmstar
312 : H_low = max(H_low, &
313 0 : pow(center_dm/base_dm, 1d0/real(n_middle - 1, dp)))
314 0 : if (H_low > H_high) then
315 0 : write(*,2) 'TDC inner zoning requires a larger mesh ratio', &
316 0 : H_low, H_high
317 0 : ierr = -1
318 0 : return
319 : end if
320 : f_low = double_ramp_mass2( &
321 0 : H_low, base_dm, center_dm, n_middle, nz_inner) - rem_mass
322 : f_high = double_ramp_mass2( &
323 0 : H_high, base_dm, center_dm, n_middle, nz_inner) - rem_mass
324 0 : if (f_low*f_high > 0d0) then
325 0 : write(*,2) 'failed to bracket TDC core zoning ramp', &
326 0 : n_middle, f_low, f_high
327 0 : ierr = -1
328 0 : return
329 : end if
330 0 : do iter = 1, 1000
331 0 : H_mid = 0.5d0*(H_low + H_high)
332 : f_mid = double_ramp_mass2( &
333 0 : H_mid, base_dm, center_dm, n_middle, nz_inner) - rem_mass
334 0 : if (abs(f_mid) < 1d-12*rem_mass) exit
335 0 : if (f_low*f_mid <= 0d0) then
336 0 : H_high = H_mid
337 0 : f_high = f_mid
338 : else
339 0 : H_low = H_mid
340 0 : f_low = f_mid
341 : end if
342 : end do
343 0 : H = H_mid
344 :
345 0 : do k = nz_outer + 1, nz_outer + n_middle
346 0 : s%dm(k) = base_dm*pow(H, real(k - nz_outer - 1, dp))
347 : end do
348 :
349 0 : peak_dm = s%dm(nz_outer + n_middle)
350 0 : H_inner = pow(peak_dm/center_dm, 1d0/real(nz_inner, dp))
351 0 : do k = nz_outer + n_middle + 1, nz_base
352 0 : s%dm(k) = center_dm*pow(H_inner, real(nz_base - k, dp))
353 : end do
354 :
355 0 : correction = s%xmstar - sum(s%dm(1:nz_base))
356 : s%dm(nz_outer + n_middle) = &
357 0 : s%dm(nz_outer + n_middle) + correction
358 : end if
359 :
360 : else ! use log zoning inward from anchor to core.
361 0 : lnx = log(xm(nz_outer + 1))
362 0 : if (is_bad(lnx)) then
363 0 : write (*, 2) 'bad lnx', nz_outer + 1, lnx, xm(nz_outer + 1)
364 0 : call mesa_error(__FILE__, __LINE__, 'set_xm_new')
365 : end if
366 0 : dlnx = (log(s%xmstar) - lnx)/(nz_base - nz_outer)
367 0 : do k = nz_outer + 2, nz_base
368 0 : lnx = lnx + dlnx
369 0 : xm(k) = exp(lnx)
370 0 : s%dm(k - 1) = xm(k) - xm(k - 1)
371 : end do
372 0 : s%dm(nz_base) = s%xmstar - xm(nz_base)
373 :
374 : ! enforce the last boundary at total mass
375 0 : xm(nz_base + 1) = s%xmstar
376 :
377 : ! recompute cell masses
378 0 : do k = nz_outer + 1, nz_base
379 0 : s%dm(k) = xm(k + 1) - xm(k)
380 : end do
381 :
382 : end if
383 :
384 0 : xm(1) = 0d0
385 0 : do k = 1, nz_base
386 0 : xm(k + 1) = xm(k) + s%dm(k)
387 : end do
388 0 : xm(nz_base + 1) = s%xmstar
389 :
390 0 : if (nz_T_gradient > 0) then
391 0 : call add_T_gradient_zones2(nz_outer, nz_base, nz_T_gradient)
392 0 : if (ierr /= 0) return
393 : end if
394 :
395 0 : do k = 1, nz
396 0 : if (s%dm(k) <= 0d0 .or. is_bad(s%dm(k))) then
397 0 : write(*,2) 'bad cell mass in TDC remesh', k, s%dm(k)
398 0 : ierr = -1
399 0 : return
400 : end if
401 : end do
402 :
403 0 : if ((.not. s%remesh_for_TDC_pulsations_log_core_zoning .and. nz_inner > 0) .or. &
404 : nz_T_gradient > 0) then
405 0 : do k = 2, nz
406 0 : if (s%dm(k) > ratio*(1d0 + 1d-10)*s%dm(k - 1) .or. &
407 0 : s%dm(k - 1) > ratio*(1d0 + 1d-10)*s%dm(k)) then
408 0 : write(*,2) 'bad adjacent cell mass ratio in TDC mesh', &
409 0 : k, s%dm(k)/s%dm(k - 1), ratio
410 0 : ierr = -1
411 0 : return
412 : end if
413 : end do
414 : end if
415 :
416 0 : xm(1) = 0d0
417 0 : do k = 1, nz
418 0 : xm(k + 1) = xm(k) + s%dm(k)
419 : end do
420 0 : xm(nz + 1) = s%xmstar
421 :
422 0 : do k = 1, nz - 1
423 0 : xm_mid(k) = 0.5d0*(xm(k) + xm(k + 1))
424 : end do
425 0 : xm_mid(nz) = 0.5d0*(xm(nz) + s%xmstar)
426 0 : s%m(1) = s%mstar
427 0 : s%q(1) = 1d0
428 0 : s%dq(1) = s%dm(1)/s%xmstar
429 0 : do k = 2, nz
430 0 : s%m(k) = s%m(k - 1) - s%dm(k - 1)
431 0 : s%dq(k) = s%dm(k)/s%xmstar
432 0 : s%q(k) = s%q(k - 1) - s%dq(k - 1)
433 : end do
434 0 : call set_dm_bar(s, nz, s%dm, s%dm_bar)
435 : end subroutine set_xm_new2
436 :
437 0 : real(dp) function geometric_sum2(H, n) result(sum_H)
438 : real(dp), intent(in) :: H
439 : integer, intent(in) :: n
440 :
441 0 : if (abs(H - 1d0) < 1d-12) then
442 0 : sum_H = real(n, dp)
443 : else
444 0 : sum_H = (pow(H, real(n, dp)) - 1d0)/(H - 1d0)
445 : end if
446 0 : end function geometric_sum2
447 :
448 0 : real(dp) function double_ramp_mass2( &
449 : H, base_dm, center_dm, n_middle, n_inner) result(total_mass)
450 : real(dp), intent(in) :: H, base_dm, center_dm
451 : integer, intent(in) :: n_middle, n_inner
452 : real(dp) :: H_inner, peak_dm
453 :
454 0 : peak_dm = base_dm*pow(H, real(n_middle - 1, dp))
455 0 : H_inner = pow(peak_dm/center_dm, 1d0/real(n_inner, dp))
456 : total_mass = base_dm*geometric_sum2(H, n_middle) + &
457 0 : center_dm*geometric_sum2(H_inner, n_inner)
458 0 : end function double_ramp_mass2
459 :
460 0 : subroutine add_T_gradient_zones2(nz_outer, nz_base, nz_T_gradient)
461 : integer, intent(in) :: nz_outer, nz_base, nz_T_gradient
462 : integer :: i, i_base, i_old, j, k_new, max_points, npts, nz_core_base
463 : real(dp) :: dlnT_total, eps_xm, fraction, next_base, next_old, &
464 : next_xm, target
465 0 : real(dp), allocatable :: base_xm(:), lnT_monitor(:), &
466 0 : monitor(:), monitor_xm(:)
467 : include 'formats'
468 :
469 0 : nz_core_base = nz_base - nz_outer
470 0 : max_points = nz_old + nz_core_base + 2
471 : allocate(base_xm(nz_base + 1), lnT_monitor(max_points), &
472 0 : monitor(max_points), monitor_xm(max_points))
473 0 : base_xm = xm(1:nz_base + 1)
474 0 : eps_xm = 16d0*epsilon(1d0)*max(abs(s%xmstar), 1d0)
475 :
476 : ! Include both old temperature points and base-grid boundaries in the monitor.
477 0 : npts = 1
478 0 : monitor_xm(npts) = xm_anchor
479 0 : i_base = nz_outer + 2
480 0 : i_old = 1
481 0 : do while (i_old <= nz_old .and. xm_mid_old(i_old) <= xm_anchor + eps_xm)
482 0 : i_old = i_old + 1
483 : end do
484 :
485 : do
486 0 : next_base = huge(1d0)
487 0 : if (i_base <= nz_base) next_base = base_xm(i_base)
488 0 : next_old = huge(1d0)
489 0 : if (i_old <= nz_old) next_old = xm_mid_old(i_old)
490 0 : next_xm = min(next_base, next_old)
491 0 : if (next_xm >= s%xmstar - eps_xm) exit
492 :
493 0 : if (next_xm > monitor_xm(npts) + eps_xm) then
494 0 : npts = npts + 1
495 0 : monitor_xm(npts) = next_xm
496 : end if
497 0 : if (next_base <= next_xm + eps_xm) i_base = i_base + 1
498 0 : if (next_old <= next_xm + eps_xm) i_old = i_old + 1
499 : end do
500 0 : npts = npts + 1
501 0 : monitor_xm(npts) = s%xmstar
502 :
503 0 : do i = 1, npts
504 0 : lnT_monitor(i) = lnT_at_xm2(monitor_xm(i))
505 : end do
506 0 : monitor(1) = 0d0
507 0 : do i = 2, npts
508 : monitor(i) = monitor(i - 1) + &
509 0 : abs(lnT_monitor(i) - lnT_monitor(i - 1))
510 : end do
511 0 : dlnT_total = monitor(npts)
512 :
513 0 : if (dlnT_total > 0d0) then
514 0 : do i = 1, npts
515 : monitor(i) = base_mesh_coordinate2( &
516 : monitor_xm(i), base_xm, nz_outer, nz_base) + &
517 0 : nz_T_gradient*monitor(i)/dlnT_total
518 : end do
519 : else
520 0 : do i = 1, npts
521 : monitor(i) = base_mesh_coordinate2( &
522 : monitor_xm(i), base_xm, nz_outer, nz_base)* &
523 0 : real(nz_core_base + nz_T_gradient, dp)/real(nz_core_base, dp)
524 : end do
525 : end if
526 0 : monitor(1) = 0d0
527 0 : monitor(npts) = real(nz_core_base + nz_T_gradient, dp)
528 :
529 0 : xm(1:nz_outer + 1) = base_xm(1:nz_outer + 1)
530 0 : j = 2
531 0 : do k_new = nz_outer + 2, nz
532 0 : target = real(k_new - nz_outer - 1, dp)
533 0 : do while (j < npts .and. monitor(j) < target)
534 0 : j = j + 1
535 : end do
536 0 : if (monitor(j) <= monitor(j - 1)) then
537 0 : write(*,2) 'bad TDC temperature-gradient monitor', &
538 0 : j, monitor(j - 1), monitor(j)
539 0 : ierr = -1
540 0 : deallocate(base_xm, lnT_monitor, monitor, monitor_xm)
541 : return
542 : end if
543 0 : fraction = (target - monitor(j - 1))/(monitor(j) - monitor(j - 1))
544 : xm(k_new) = monitor_xm(j - 1) + &
545 0 : fraction*(monitor_xm(j) - monitor_xm(j - 1))
546 : end do
547 0 : xm(nz + 1) = s%xmstar
548 0 : do i = nz_outer + 1, nz
549 0 : s%dm(i) = xm(i + 1) - xm(i)
550 : end do
551 :
552 0 : deallocate(base_xm, lnT_monitor, monitor, monitor_xm)
553 : end subroutine add_T_gradient_zones2
554 :
555 0 : real(dp) function base_mesh_coordinate2( &
556 0 : x, base_xm, nz_outer, nz_base) result(coordinate)
557 : real(dp), intent(in) :: x, base_xm(:)
558 : integer, intent(in) :: nz_outer, nz_base
559 : integer :: i
560 :
561 0 : if (x <= base_xm(nz_outer + 1)) then
562 0 : coordinate = 0d0
563 : return
564 : end if
565 0 : do i = nz_outer + 1, nz_base
566 0 : if (x <= base_xm(i + 1)) then
567 : coordinate = real(i - nz_outer - 1, dp) + &
568 0 : (x - base_xm(i))/(base_xm(i + 1) - base_xm(i))
569 0 : return
570 : end if
571 : end do
572 0 : coordinate = real(nz_base - nz_outer, dp)
573 0 : end function base_mesh_coordinate2
574 :
575 0 : real(dp) function lnT_at_xm2(x) result(lnT_at_xm)
576 : real(dp), intent(in) :: x
577 : integer :: i
578 : real(dp) :: fraction
579 :
580 0 : if (nz_old == 1) then
581 0 : lnT_at_xm = s%xh(s%i_lnT, 1)
582 0 : return
583 : end if
584 0 : if (x <= xm_mid_old(1)) then
585 : i = 2
586 : else
587 0 : do i = 2, nz_old
588 0 : if (x <= xm_mid_old(i)) exit
589 : end do
590 0 : if (i > nz_old) i = nz_old
591 : end if
592 : fraction = (x - xm_mid_old(i - 1))/ &
593 0 : (xm_mid_old(i) - xm_mid_old(i - 1))
594 : lnT_at_xm = s%xh(s%i_lnT, i - 1) + fraction* &
595 0 : (s%xh(s%i_lnT, i) - s%xh(s%i_lnT, i - 1))
596 0 : end function lnT_at_xm2
597 :
598 0 : subroutine interpolate1_face_val2(i, cntr_val)
599 : integer, intent(in) :: i
600 : real(dp), intent(in) :: cntr_val
601 0 : do k = 1, nz_old
602 0 : v_old(k) = s%xh(i, k)
603 : end do
604 0 : v_old(nz_old + 1) = cntr_val
605 : call interpolate_vector_pm( &
606 0 : nz_old + 1, xm_old, nz + 1, xm, v_old, v_new, work1, 'remesh_for_TDC', ierr)
607 0 : if (ierr /= 0) call mesa_error(__FILE__, __LINE__, 'TDC remesh face interpolation failed')
608 0 : do k = 1, nz
609 0 : s%xh(i, k) = v_new(k)
610 : end do
611 0 : end subroutine interpolate1_face_val2
612 :
613 0 : subroutine check_new_lnR2
614 : include 'formats'
615 0 : do k = 1, nz
616 0 : s%lnR(k) = s%xh(s%i_lnR, k)
617 0 : s%r(k) = exp(s%lnR(k))
618 : end do
619 0 : do k = 1, nz - 1
620 0 : if (s%r(k) <= s%r(k + 1)) then
621 0 : write (*, 2) 'bad r', k, s%r(k), s%r(k + 1)
622 0 : call mesa_error(__FILE__, __LINE__, 'check_new_lnR remesh TDC')
623 : end if
624 : end do
625 0 : if (s%r(nz) <= s%r_center) then
626 0 : write (*, 2) 'bad r center', nz, s%r(nz), s%r_center
627 0 : call mesa_error(__FILE__, __LINE__, 'check_new_lnR remesh TDC')
628 : end if
629 0 : end subroutine check_new_lnR2
630 :
631 0 : subroutine set_new_lnd2
632 : real(dp) :: vol, r300, r3p1
633 : include 'formats'
634 0 : do k = 1, nz
635 0 : r300 = pow3(s%r(k))
636 0 : if (k < nz) then
637 0 : r3p1 = pow3(s%r(k + 1))
638 : else
639 0 : r3p1 = pow3(s%r_center)
640 : end if
641 0 : vol = (4d0*pi/3d0)*(r300 - r3p1)
642 0 : s%rho(k) = s%dm(k)/vol
643 0 : s%lnd(k) = log(s%rho(k))
644 0 : s%xh(s%i_lnd, k) = s%lnd(k)
645 0 : if (is_bad(s%lnd(k))) then
646 0 : write (*, 2) 'bad lnd vol dm r300 r3p1', k, s%lnd(k), vol, s%dm(k), r300, r3p1
647 0 : call mesa_error(__FILE__, __LINE__, 'remesh for TDC')
648 : end if
649 : end do
650 0 : end subroutine set_new_lnd2
651 :
652 0 : subroutine interpolate1_cell_val2(i)
653 : integer, intent(in) :: i
654 0 : do k = 1, nz_old
655 0 : v_old(k) = s%xh(i, k)
656 : end do
657 : call interpolate_vector_pm( &
658 0 : nz_old, xm_mid_old, nz, xm_mid, v_old, v_new, work1, 'remesh_for_TDC', ierr)
659 0 : if (ierr /= 0) call mesa_error(__FILE__, __LINE__, 'TDC remesh cell interpolation failed')
660 0 : do k = 1, nz
661 0 : s%xh(i, k) = v_new(k)
662 : end do
663 0 : end subroutine interpolate1_cell_val2
664 :
665 0 : subroutine remap1_xa2(j)
666 : integer, intent(in) :: j
667 : integer :: k_old, k_scan, k_new
668 : real(dp) :: overlap, species_mass
669 :
670 0 : v_old(1:nz_old) = s%xa(j, 1:nz_old)
671 0 : k_old = 1
672 0 : do k_new = 1, nz
673 0 : do while (k_old < nz_old .and. xm_old(k_old + 1) <= xm(k_new))
674 0 : k_old = k_old + 1
675 : end do
676 : species_mass = 0d0
677 : k_scan = k_old
678 0 : do while (k_scan <= nz_old .and. xm_old(k_scan) < xm(k_new + 1))
679 : ! A mass-overlap average conserves the total mass of each species.
680 : overlap = min(xm(k_new + 1), xm_old(k_scan + 1)) - &
681 0 : max(xm(k_new), xm_old(k_scan))
682 0 : if (overlap > 0d0) species_mass = species_mass + overlap*v_old(k_scan)
683 0 : k_scan = k_scan + 1
684 : end do
685 0 : s%xa(j, k_new) = species_mass/s%dm(k_new)
686 : end do
687 0 : end subroutine remap1_xa2
688 :
689 0 : subroutine remap_rotation2(old_J, old_abs_J)
690 : real(dp), intent(in) :: old_J, old_abs_J
691 : real(dp) :: new_J
692 : include 'formats'
693 :
694 0 : v_old(1:nz_old) = s%j_rot(1:nz_old)
695 0 : call remap1_face_average2
696 0 : s%j_rot(1:nz) = v_new(1:nz)
697 0 : if (s%i_j_rot /= 0) s%xh(s%i_j_rot,1:nz) = s%j_rot(1:nz)
698 :
699 0 : v_old(1:nz_old) = s%omega(1:nz_old)
700 0 : call remap1_face_average2
701 0 : s%omega(1:nz) = v_new(1:nz)
702 :
703 0 : new_J = dot_product(s%dm_bar(1:nz), s%j_rot(1:nz))
704 0 : if (abs(new_J - old_J) > 1d-12*max(old_abs_J, tiny(1d0))) then
705 0 : write(*,1) 'relative angular momentum error in TDC remesh', &
706 0 : (new_J - old_J)/max(abs(old_J), old_abs_J, tiny(1d0))
707 0 : call mesa_error(__FILE__, __LINE__, 'TDC remesh failed to conserve angular momentum')
708 : end if
709 0 : s%total_angular_momentum = new_J
710 : s%total_abs_angular_momentum = &
711 0 : dot_product(s%dm_bar(1:nz), abs(s%j_rot(1:nz)))
712 0 : end subroutine remap_rotation2
713 :
714 0 : subroutine remap1_face_average2
715 : integer :: k_old, k_scan, k_new
716 : real(dp) :: new_outer, new_inner, old_outer, old_inner, &
717 : overlap, integral
718 :
719 0 : k_old = 1
720 0 : do k_new = 1, nz
721 0 : if (k_new == 1) then
722 : new_outer = 0d0
723 : else
724 0 : new_outer = xm_mid(k_new - 1)
725 : end if
726 0 : if (k_new == nz) then
727 0 : new_inner = s%xmstar
728 : else
729 0 : new_inner = xm_mid(k_new)
730 : end if
731 :
732 0 : do while (k_old < nz_old .and. xm_mid_old(k_old) <= new_outer)
733 0 : k_old = k_old + 1
734 : end do
735 : integral = 0d0
736 : k_scan = k_old
737 0 : do while (k_scan <= nz_old)
738 0 : if (k_scan == 1) then
739 : old_outer = 0d0
740 : else
741 0 : old_outer = xm_mid_old(k_scan - 1)
742 : end if
743 0 : if (k_scan == nz_old) then
744 0 : old_inner = s%xmstar
745 : else
746 0 : old_inner = xm_mid_old(k_scan)
747 : end if
748 0 : if (old_outer >= new_inner) exit
749 0 : overlap = min(new_inner, old_inner) - max(new_outer, old_outer)
750 0 : if (overlap > 0d0) integral = integral + overlap*v_old(k_scan)
751 0 : k_scan = k_scan + 1
752 : end do
753 0 : v_new(k_new) = integral/(new_inner - new_outer)
754 : end do
755 0 : end subroutine remap1_face_average2
756 :
757 0 : subroutine update_composition_info2
758 : use chem_lib, only: basic_composition_info
759 : real(dp) :: sumx
760 :
761 0 : do k = 1, nz
762 : call basic_composition_info( &
763 : s%species, s%chem_id, s%xa(1:s%species,k), &
764 : s%X(k), s%Y(k), s%Z(k), s%abar(k), s%zbar(k), s%z2bar(k), &
765 0 : s%z53bar(k), s%ye(k), s%mass_correction(k), sumx)
766 : end do
767 0 : end subroutine update_composition_info2
768 :
769 0 : subroutine set_rotation_seed2
770 : use hydro_rotation, only: w_div_w_roche_jrot
771 :
772 0 : do k = 1, nz
773 : s%w_div_w_crit_roche(k) = &
774 : w_div_w_roche_jrot(s%r(k), s%m(k), s%j_rot(k), s%cgrav(k), &
775 0 : s%w_div_wcrit_max, s%w_div_wcrit_max2, s%w_div_wc_flag)
776 : end do
777 0 : if (s%i_w_div_wc /= 0) &
778 0 : s%xh(s%i_w_div_wc,1:nz) = s%w_div_w_crit_roche(1:nz)
779 0 : end subroutine set_rotation_seed2
780 :
781 0 : subroutine revise_lnT_for_QHSE2(P_surf, ierr)
782 : use eos_def, only: num_eos_basic_results, num_eos_d_dxa_results
783 : use chem_def, only: chem_isos
784 : use eos_support, only: solve_eos_given_DP
785 : use eos_def, only: i_eta, i_lnfree_e
786 : use kap_def, only: num_kap_fracs
787 : use kap_support, only: get_kap
788 : real(dp), intent(in) :: P_surf
789 : integer, intent(out) :: ierr
790 : real(dp) :: logRho, logP, logT_guess, &
791 : logT_tol, logP_tol, logT, P_m1, P_00, dm_face, &
792 : kap_fracs(num_kap_fracs), kap, dlnkap_dlnRho, dlnkap_dlnT, &
793 : old_kap, new_P_surf, new_T_surf
794 : real(dp), dimension(num_eos_basic_results) :: &
795 : res, d_dlnd, d_dlnT
796 0 : real(dp) :: dres_dxa(num_eos_d_dxa_results, s%species)
797 : include 'formats'
798 0 : ierr = 0
799 0 : P_m1 = P_surf
800 0 : do k = 1, nz
801 0 : s%lnT(k) = s%xh(s%i_lnT, k)
802 0 : s%lnR(k) = s%xh(s%i_lnR, k)
803 0 : s%r(k) = exp(s%lnR(k))
804 : end do
805 0 : do k = 1, nz
806 0 : if (k < nz) then
807 0 : dm_face = s%dm_bar(k)
808 : else
809 0 : dm_face = 0.5d0*(s%dm(k - 1) + s%dm(k))
810 : end if
811 0 : P_00 = P_m1 + s%cgrav(k)*s%m(k)*dm_face/(4d0*pi*pow4(s%r(k)))
812 0 : logP = log10(P_00) ! value for QHSE
813 0 : s%lnPeos(k) = logP*ln10
814 0 : s%Peos(k) = P_00
815 0 : logRho = s%lnd(k)/ln10
816 0 : logT_guess = s%lnT(k)/ln10
817 0 : logT_tol = 1d-11
818 0 : logP_tol = 1d-11
819 : call solve_eos_given_DP( &
820 : s, k, s%xa(:, k), &
821 : logRho, logP, logT_guess, logT_tol, logP_tol, &
822 0 : logT, res, d_dlnd, d_dlnT, dres_dxa, ierr)
823 0 : if (ierr /= 0) then
824 0 : write (*, 2) 'solve_eos_given_DP failed', k
825 0 : write (*, '(A)')
826 0 : write (*, 1) 'sum(xa)', sum(s%xa(:, k))
827 0 : do j = 1, s%species
828 0 : write (*, 4) 'xa(j,k) '//trim(chem_isos%name(s%chem_id(j))), j, j + s%nvar_hydro, k, s%xa(j, k)
829 : end do
830 0 : write (*, 1) 'logRho', logRho
831 0 : write (*, 1) 'logP', logP
832 0 : write (*, 1) 'logT_guess', logT_guess
833 0 : write (*, 1) 'logT_tol', logT_tol
834 0 : write (*, 1) 'logP_tol', logP_tol
835 0 : write (*, '(A)')
836 0 : call mesa_error(__FILE__, __LINE__, 'revise_lnT_for_QHSE')
837 : end if
838 0 : s%lnT(k) = logT*ln10
839 0 : s%xh(s%i_lnT, k) = s%lnT(k)
840 0 : P_m1 = P_00
841 :
842 0 : if (k == 1) then ! get opacity and recheck surf BCs
843 : call get_kap( &
844 : s, k, s%zbar(k), s%xa(:, k), logRho, logT, &
845 : res(i_lnfree_e), d_dlnd(i_lnfree_e), d_dlnT(i_lnfree_e), &
846 : res(i_eta), d_dlnd(i_eta), d_dlnT(i_eta), &
847 : kap_fracs, kap, dlnkap_dlnRho, dlnkap_dlnT, &
848 0 : ierr)
849 0 : if (ierr /= 0) then
850 0 : write (*, 2) 'get_kap failed', k
851 0 : call mesa_error(__FILE__, __LINE__, 'revise_lnT_for_QHSE')
852 : end if
853 0 : old_kap = s%opacity(1)
854 0 : s%opacity(1) = kap ! for use by atm surf PT
855 0 : call get_PT_surf2(new_P_surf, new_T_surf, ierr)
856 0 : if (ierr /= 0) then
857 0 : write (*, 2) 'get_PT_surf failed', k
858 0 : call mesa_error(__FILE__, __LINE__, 'revise_lnT_for_QHSE')
859 : end if
860 0 : write (*, 1) 'new old T_surf', new_T_surf, T_surf
861 0 : write (*, 1) 'new old P_surf', new_P_surf, P_surf
862 0 : write (*, 1) 'new old kap(1)', kap, old_kap
863 : end if
864 :
865 : end do
866 0 : end subroutine revise_lnT_for_QHSE2
867 :
868 : end subroutine remesh_for_TDC_pulsations
869 :
870 : end module tdc_hydro_support
|
