Line data Source code
1 : ! ***********************************************************************
2 : !
3 : ! Copyright (C) 2010-2020 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 hydro_rsp2_support
21 :
22 : use star_private_def
23 : use const_def, only: dp, ln10, pi, lsun, rsun, msun
24 : use utils_lib, only: is_bad
25 : use auto_diff
26 : use auto_diff_support
27 : use accurate_sum_auto_diff_star_order1
28 : use star_utils
29 :
30 : implicit none
31 :
32 : private
33 : public :: remesh_for_RSP2
34 :
35 : contains
36 :
37 0 : subroutine remesh_for_RSP2(s,ierr)
38 : ! uses these controls
39 : ! RSP2_nz = 150
40 : ! RSP2_nz_outer = 40
41 : ! RSP2_T_anchor = 11d3
42 : ! RSP2_dq_1_factor = 2d0
43 : use interp_1d_def, only: pm_work_size
44 : use interp_1d_lib, only: interpolate_vector_pm
45 : type (star_info), pointer :: s
46 : integer, intent(out) :: ierr
47 : integer :: k, j, nz_old, nz
48 0 : real(dp) :: xm_anchor, P_surf, T_surf, old_L1, old_r1
49 : real(dp), allocatable, dimension(:) :: &
50 0 : xm_old, xm, xm_mid_old, xm_mid, v_old, v_new
51 0 : real(dp), pointer :: work1(:) ! =(nz_old+1, pm_work_size)
52 : include 'formats'
53 0 : ierr = 0
54 0 : nz_old = s% nz
55 0 : nz = s% RSP2_nz
56 0 : if (nz == nz_old) return ! assume have already done remesh for RSP2
57 0 : if (nz > nz_old) call mesa_error(__FILE__,__LINE__,'remesh_for_RSP2 cannot increase nz')
58 0 : call setvars(ierr)
59 0 : if (ierr /= 0) call mesa_error(__FILE__,__LINE__,'remesh_for_RSP2 failed in setvars')
60 0 : old_L1 = s% L(1)
61 0 : old_r1 = s% r(1)
62 0 : call set_phot_info(s) ! sets Teff
63 0 : call get_PT_surf(P_surf, T_surf, ierr)
64 0 : if (ierr /= 0) call mesa_error(__FILE__,__LINE__,'remesh_for_RSP2 failed in get_PT_surf')
65 : allocate(&
66 0 : xm_old(nz_old+1), xm_mid_old(nz_old), v_old(nz_old+1), &
67 0 : xm(nz+1), xm_mid(nz), v_new(nz+1), work1((nz_old+1)*pm_work_size))
68 0 : call set_xm_old
69 0 : call find_xm_anchor
70 0 : call set_xm_new
71 0 : call interpolate1_face_val(s% i_lnR, log(max(1d0,s% r_center)))
72 0 : call check_new_lnR
73 0 : call interpolate1_face_val(s% i_lum, s% L_center)
74 0 : if (s% i_v /= 0) call interpolate1_face_val(s% i_v, s% v_center)
75 0 : call set_new_lnd
76 0 : call interpolate1_cell_val(s% i_lnT)
77 0 : call interpolate1_cell_val(s% i_w)
78 0 : do j=1,s% species
79 0 : call interpolate1_xa(j)
80 : end do
81 0 : call rescale_xa
82 0 : call revise_lnT_for_QHSE(P_surf, ierr)
83 0 : if (ierr /= 0) call mesa_error(__FILE__,__LINE__,'remesh_for_RSP2 failed in revise_lnT_for_QHSE')
84 0 : do k=1,nz
85 0 : call set_Hp_face(k)
86 : end do
87 0 : deallocate(work1)
88 0 : s% nz = nz
89 0 : write(*,1) 'new old L_surf/Lsun', s% xh(s% i_lum,1)/Lsun, old_L1/Lsun
90 0 : write(*,1) 'new old R_surf/Rsun', exp(s% xh(s% i_lnR,1))/Rsun, old_r1/Rsun
91 0 : write(*,'(A)')
92 : !call mesa_error(__FILE__,__LINE__,'remesh_for_RSP2')
93 :
94 : contains
95 :
96 0 : subroutine setvars(ierr)
97 0 : use hydro_vars, only: unpack_xh, set_hydro_vars
98 : integer, intent(out) :: ierr
99 : logical, parameter :: &
100 : skip_basic_vars = .false., &
101 : skip_micro_vars = .false., &
102 : skip_m_grav_and_grav = .false., &
103 : skip_net = .true., &
104 : skip_neu = .true., &
105 : skip_kap = .false., &
106 : skip_grads = .true., &
107 : skip_rotation = .true., &
108 : skip_brunt = .true., &
109 : skip_other_cgrav = .true., &
110 : skip_mixing_info = .true., &
111 : skip_set_cz_bdy_mass = .true., &
112 : skip_mlt = .true., &
113 : skip_eos = .false.
114 : ierr = 0
115 0 : call unpack_xh(s,ierr)
116 0 : if (ierr /= 0) call mesa_error(__FILE__,__LINE__,'remesh_for_RSP2 failed in unpack_xh')
117 : call set_hydro_vars( &
118 : s, 1, nz_old, skip_basic_vars, &
119 : skip_micro_vars, skip_m_grav_and_grav, skip_eos, skip_net, skip_neu, &
120 : skip_kap, skip_grads, skip_rotation, skip_brunt, skip_other_cgrav, &
121 0 : skip_mixing_info, skip_set_cz_bdy_mass, skip_mlt, ierr)
122 0 : if (ierr /= 0) call mesa_error(__FILE__,__LINE__,'remesh_for_RSP2 failed in set_hydro_vars')
123 0 : end subroutine setvars
124 :
125 0 : subroutine get_PT_surf(P_surf, T_surf, ierr)
126 0 : use atm_support, only: get_atm_PT
127 : real(dp), intent(out) :: P_surf, T_surf
128 : integer, intent(out) :: ierr
129 : real(dp) :: &
130 0 : Teff, lnT_surf, dlnT_dL, dlnT_dlnR, dlnT_dlnM, dlnT_dlnkap, &
131 0 : lnP_surf, dlnP_dL, dlnP_dlnR, dlnP_dlnM, dlnP_dlnkap
132 : logical, parameter :: skip_partials = .true.
133 : include 'formats'
134 0 : ierr = 0
135 0 : call set_phot_info(s) ! sets s% Teff
136 0 : Teff = s% Teff
137 : call get_atm_PT( & ! this uses s% opacity(1)
138 : s, s% tau_factor*s% tau_base, s% L(1), s% r(1), s% m(1), s% cgrav(1), skip_partials, &
139 : Teff, lnT_surf, dlnT_dL, dlnT_dlnR, dlnT_dlnM, dlnT_dlnkap, &
140 0 : lnP_surf, dlnP_dL, dlnP_dlnR, dlnP_dlnM, dlnP_dlnkap, ierr)
141 0 : if (ierr /= 0) call mesa_error(__FILE__,__LINE__,'get_P_surf failed in get_atm_PT')
142 0 : P_surf = exp(lnP_surf)
143 0 : T_surf = exp(lnT_surf)
144 : return
145 :
146 : write(*,1) 'get_PT_surf P_surf', P_surf
147 : write(*,1) 'get_PT_surf T_surf', T_surf
148 : write(*,1) 'get_PT_surf Teff', Teff
149 : write(*,1) 'get_PT_surf opacity(1)', s% opacity(1)
150 : write(*,1)
151 : !call mesa_error(__FILE__,__LINE__,'get_PT_surf')
152 0 : end subroutine get_PT_surf
153 :
154 0 : subroutine set_xm_old
155 0 : xm_old(1) = 0d0
156 0 : do k=2,nz_old
157 0 : xm_old(k) = xm_old(k-1) + s% dm(k-1)
158 : end do
159 0 : xm_old(nz_old+1) = s% xmstar
160 0 : do k=1,nz_old
161 0 : xm_mid_old(k) = xm_old(k) + 0.5d0*s% dm(k)
162 : end do
163 0 : end subroutine set_xm_old
164 :
165 0 : subroutine find_xm_anchor
166 0 : real(dp) :: lnT_anchor, xmm1, xm00, lnTm1, lnT00
167 : include 'formats'
168 0 : lnT_anchor = log(s% RSP2_T_anchor)
169 0 : if (lnT_anchor <= s% xh(s% i_lnT,1)) then
170 0 : write(*,1) 'T_anchor < T_surf', s% RSP2_T_anchor, exp(s% xh(s% i_lnT,1))
171 0 : call mesa_error(__FILE__,__LINE__,'find_xm_anchor')
172 : end if
173 0 : xm_anchor = xm_old(nz_old)
174 0 : do k=2,nz_old
175 0 : if (s% xh(s% i_lnT,k) >= lnT_anchor) then
176 0 : xmm1 = xm_old(k-1)
177 0 : xm00 = xm_old(k)
178 0 : lnTm1 = s% xh(s% i_lnT,k-1)
179 0 : lnT00 = s% xh(s% i_lnT,k)
180 : xm_anchor = xmm1 + &
181 0 : (xm00 - xmm1)*(lnT_anchor - lnTm1)/(lnT00 - lnTm1)
182 0 : if (is_bad(xm_anchor) .or. xm_anchor <= 0d0) then
183 0 : write(*,2) 'bad xm_anchor', k, xm_anchor, xmm1, xm00, lnTm1, lnT00, lnT_anchor, s% lnT(1)
184 0 : call mesa_error(__FILE__,__LINE__,'find_xm_anchor')
185 : end if
186 0 : return
187 : end if
188 : end do
189 : end subroutine find_xm_anchor
190 :
191 0 : subroutine set_xm_new ! sets xm, dm, m, dq, q
192 : integer :: nz_outer, k
193 0 : real(dp) :: dq_1_factor, dxm_outer, lnx, dlnx
194 : include 'formats'
195 0 : nz_outer = s% RSP2_nz_outer
196 0 : dq_1_factor = s% RSP2_dq_1_factor
197 0 : dxm_outer = xm_anchor/(nz_outer - 1d0 + dq_1_factor)
198 : !write(*,2) 'dxm_outer', nz_outer, dxm_outer, xm_anchor
199 0 : xm(1) = 0d0
200 0 : xm(2) = dxm_outer*dq_1_factor
201 0 : s% dm(1) = xm(2)
202 0 : do k=3,nz_outer+1
203 0 : xm(k) = xm(k-1) + dxm_outer
204 0 : s% dm(k-1) = dxm_outer
205 : end do
206 0 : lnx = log(xm(nz_outer+1))
207 0 : if (is_bad(lnx)) then
208 0 : write(*,2) 'bad lnx', nz_outer+1, lnx, xm(nz_outer+1)
209 0 : call mesa_error(__FILE__,__LINE__,'set_xm_new')
210 : end if
211 0 : dlnx = (log(s% xmstar) - lnx)/(nz - nz_outer)
212 0 : do k=nz_outer+2,nz
213 0 : lnx = lnx + dlnx
214 0 : xm(k) = exp(lnx)
215 0 : s% dm(k-1) = xm(k) - xm(k-1)
216 : end do
217 0 : s% dm(nz) = s% xmstar - xm(nz)
218 0 : do k=1,nz-1
219 0 : xm_mid(k) = 0.5d0*(xm(k) + xm(k+1))
220 : end do
221 0 : xm_mid(nz) = 0.5d0*(xm(nz) + s% xmstar)
222 0 : s% m(1) = s% mstar
223 0 : s% q(1) = 1d0
224 0 : s% dq(1) = s% dm(1)/s% xmstar
225 0 : do k=2,nz
226 0 : s% m(k) = s% m(k-1) - s% dm(k-1)
227 0 : s% dq(k) = s% dm(k)/s% xmstar
228 0 : s% q(k) = s% q(k-1) - s% dq(k-1)
229 : end do
230 0 : call set_dm_bar(s, s% nz, s% dm, s% dm_bar)
231 0 : return
232 :
233 : do k=2,nz
234 : write(*,2) 'dm(k)/dm(k-1) m(k)', k, s%dm(k)/s%dm(k-1), s%m(k)/Msun
235 : end do
236 : write(*,1) 'm_center', s% m_center/msun
237 : call mesa_error(__FILE__,__LINE__,'set_xm_new')
238 : end subroutine set_xm_new
239 :
240 0 : subroutine interpolate1_face_val(i, cntr_val)
241 : integer, intent(in) :: i
242 : real(dp), intent(in) :: cntr_val
243 0 : do k=1,nz_old
244 0 : v_old(k) = s% xh(i,k)
245 : end do
246 0 : v_old(nz_old+1) = cntr_val
247 : call interpolate_vector_pm( &
248 0 : nz_old+1, xm_old, nz+1, xm, v_old, v_new, work1, 'remesh_for_RSP2', ierr)
249 0 : do k=1,nz
250 0 : s% xh(i,k) = v_new(k)
251 : end do
252 0 : end subroutine interpolate1_face_val
253 :
254 0 : subroutine check_new_lnR
255 : include 'formats'
256 0 : do k=1,nz
257 0 : s% lnR(k) = s% xh(s% i_lnR,k)
258 0 : s% r(k) = exp(s% lnR(k))
259 : end do
260 0 : do k=1,nz-1
261 0 : if (s% r(k) <= s% r(k+1)) then
262 0 : write(*,2) 'bad r', k, s% r(k), s% r(k+1)
263 0 : call mesa_error(__FILE__,__LINE__,'check_new_lnR remesh rsp2')
264 : end if
265 : end do
266 0 : if (s% r(nz) <= s% r_center) then
267 0 : write(*,2) 'bad r center', nz, s% r(nz), s% r_center
268 0 : call mesa_error(__FILE__,__LINE__,'check_new_lnR remesh rsp2')
269 : end if
270 0 : end subroutine check_new_lnR
271 :
272 0 : subroutine set_new_lnd
273 0 : real(dp) :: vol, r300, r3p1
274 : include 'formats'
275 0 : do k=1,nz
276 0 : r300 = pow3(s% r(k))
277 0 : if (k < nz) then
278 0 : r3p1 = pow3(s% r(k+1))
279 : else
280 0 : r3p1 = pow3(s% r_center)
281 : end if
282 0 : vol = (4d0*pi/3d0)*(r300 - r3p1)
283 0 : s% rho(k) = s% dm(k)/vol
284 0 : s% lnd(k) = log(s% rho(k))
285 0 : s% xh(s% i_lnd,k) = s% lnd(k)
286 0 : if (is_bad(s% lnd(k))) then
287 0 : write(*,2) 'bad lnd vol dm r300 r3p1', k, s% lnd(k), vol, s% dm(k), r300, r3p1
288 0 : call mesa_error(__FILE__,__LINE__,'remesh for rsp2')
289 : end if
290 : end do
291 0 : end subroutine set_new_lnd
292 :
293 0 : subroutine interpolate1_cell_val(i)
294 : integer, intent(in) :: i
295 0 : do k=1,nz_old
296 0 : v_old(k) = s% xh(i,k)
297 : end do
298 : call interpolate_vector_pm( &
299 0 : nz_old, xm_mid_old, nz, xm_mid, v_old, v_new, work1, 'remesh_for_RSP2', ierr)
300 0 : do k=1,nz
301 0 : s% xh(i,k) = v_new(k)
302 : end do
303 0 : end subroutine interpolate1_cell_val
304 :
305 0 : subroutine interpolate1_xa(j)
306 : integer, intent(in) :: j
307 0 : do k=1,nz_old
308 0 : v_old(k) = s% xa(j,k)
309 : end do
310 : call interpolate_vector_pm( &
311 0 : nz_old, xm_mid_old, nz, xm_mid, v_old, v_new, work1, 'remesh_for_RSP2', ierr)
312 0 : do k=1,nz
313 0 : s% xa(j,k) = v_new(k)
314 : end do
315 0 : end subroutine interpolate1_xa
316 :
317 0 : subroutine rescale_xa
318 : integer :: k, j
319 0 : real(dp) :: sum_xa
320 0 : do k=1,nz
321 0 : sum_xa = sum(s% xa(1:s% species,k))
322 0 : do j=1,s% species
323 0 : s% xa(j,k) = s% xa(j,k)/sum_xa
324 : end do
325 : end do
326 0 : end subroutine rescale_xa
327 :
328 0 : subroutine revise_lnT_for_QHSE(P_surf, ierr)
329 : use eos_def, only: num_eos_basic_results, num_eos_d_dxa_results
330 : use chem_def, only: chem_isos
331 : use eos_support, only: solve_eos_given_DP
332 : use eos_def, only: i_eta, i_lnfree_e
333 : use kap_def, only: num_kap_fracs
334 : use kap_support, only: get_kap
335 : real(dp), intent(in) :: P_surf
336 : integer, intent(out) :: ierr
337 0 : real(dp) :: logRho, logP, logT_guess, &
338 0 : logT_tol, logP_tol, logT, P_m1, P_00, dm_face, &
339 0 : kap_fracs(num_kap_fracs), kap, dlnkap_dlnRho, dlnkap_dlnT, &
340 0 : old_kap, new_P_surf, new_T_surf
341 : real(dp), dimension(num_eos_basic_results) :: &
342 0 : res, d_dlnd, d_dlnT
343 0 : real(dp) :: dres_dxa(num_eos_d_dxa_results,s% species)
344 : include 'formats'
345 0 : ierr = 0
346 0 : P_m1 = P_surf
347 0 : do k=1,nz
348 0 : s% lnT(k) = s% xh(s% i_lnT,k)
349 0 : s% lnR(k) = s% xh(s% i_lnR,k)
350 0 : s% r(k) = exp(s% lnR(k))
351 : end do
352 : !write(*,1) 'before revise_lnT_for_QHSE: logT cntr', s% lnT(nz)/ln10
353 0 : do k=1,nz
354 0 : if (k < nz) then
355 0 : dm_face = s% dm_bar(k)
356 : else
357 0 : dm_face = 0.5d0*(s% dm(k-1) + s% dm(k))
358 : end if
359 0 : P_00 = P_m1 + s% cgrav(k)*s% m(k)*dm_face/(4d0*pi*pow4(s% r(k)))
360 0 : logP = log10(P_00) ! value for QHSE
361 0 : s% lnPeos(k) = logP/ln10
362 0 : s% Peos(k) = P_00
363 0 : logRho = s% lnd(k)/ln10
364 0 : logT_guess = s% lnT(k)/ln10
365 0 : logT_tol = 1d-11
366 0 : logP_tol = 1d-11
367 : call solve_eos_given_DP( &
368 : s, k, s% xa(:,k), &
369 : logRho, logP, logT_guess, logT_tol, logP_tol, &
370 0 : logT, res, d_dlnd, d_dlnT, dres_dxa, ierr)
371 0 : if (ierr /= 0) then
372 0 : write(*,2) 'solve_eos_given_DP failed', k
373 0 : write(*,'(A)')
374 0 : write(*,1) 'sum(xa)', sum(s% xa(:,k))
375 0 : do j=1,s% species
376 0 : write(*,4) 'xa(j,k) ' // trim(chem_isos% name(s% chem_id(j))), j, j+s% nvar_hydro, k, s% xa(j,k)
377 : end do
378 0 : write(*,1) 'logRho', logRho
379 0 : write(*,1) 'logP', logP
380 0 : write(*,1) 'logT_guess', logT_guess
381 0 : write(*,1) 'logT_tol', logT_tol
382 0 : write(*,1) 'logP_tol', logP_tol
383 0 : write(*,'(A)')
384 0 : call mesa_error(__FILE__,__LINE__,'revise_lnT_for_QHSE')
385 : end if
386 0 : s% lnT(k) = logT*ln10
387 0 : s% xh(s% i_lnT,k) = s% lnT(k)
388 : !write(*,2) 'logP dlogT logT logT_guess logRho', k, &
389 : ! logP, logT - logT_guess, logT, logT_guess, logRho
390 0 : P_m1 = P_00
391 :
392 0 : if (k == 1) then ! get opacity and recheck surf BCs
393 : call get_kap( & ! assume zbar is set
394 : s, k, s% zbar(k), s% xa(:,k), logRho, logT, &
395 : res(i_lnfree_e), d_dlnd(i_lnfree_e), d_dlnT(i_lnfree_e), &
396 : res(i_eta), d_dlnd(i_eta), d_dlnT(i_eta), &
397 : kap_fracs, kap, dlnkap_dlnRho, dlnkap_dlnT, &
398 0 : ierr)
399 0 : if (ierr /= 0) then
400 0 : write(*,2) 'get_kap failed', k
401 0 : call mesa_error(__FILE__,__LINE__,'revise_lnT_for_QHSE')
402 : end if
403 0 : old_kap = s% opacity(1)
404 0 : s% opacity(1) = kap ! for use by atm surf PT
405 0 : call get_PT_surf(new_P_surf, new_T_surf, ierr)
406 0 : if (ierr /= 0) then
407 0 : write(*,2) 'get_PT_surf failed', k
408 0 : call mesa_error(__FILE__,__LINE__,'revise_lnT_for_QHSE')
409 : end if
410 0 : write(*,1) 'new old T_surf', new_T_surf, T_surf
411 0 : write(*,1) 'new old P_surf', new_P_surf, P_surf
412 0 : write(*,1) 'new old kap(1)', kap, old_kap
413 : !call mesa_error(__FILE__,__LINE__,'revise_lnT_for_QHSE')
414 : end if
415 :
416 : end do
417 : !write(*,1) 'after revise_lnT_for_QHSE: logT cntr', s% lnT(nz)/ln10
418 : !stop
419 0 : end subroutine revise_lnT_for_QHSE
420 :
421 0 : subroutine set_Hp_face(k)
422 0 : use hydro_rsp2, only: get_RSP2_alfa_beta_face_weights
423 : integer, intent(in) :: k
424 0 : real(dp) :: r_00, d_00, Peos_00, Peos_div_rho, Hp_face, &
425 0 : d_m1, Peos_m1, alfa, beta
426 0 : r_00 = s% r(k)
427 0 : d_00 = s% rho(k)
428 0 : Peos_00 = s% Peos(k)
429 0 : if (k == 1) then
430 0 : Peos_div_rho = Peos_00/d_00
431 0 : Hp_face = pow2(r_00)*Peos_div_rho/(s% cgrav(k)*s% m(k))
432 : else
433 0 : d_m1 = s% rho(k-1)
434 0 : Peos_m1 = s% Peos(k-1)
435 0 : call get_RSP2_alfa_beta_face_weights(s, k, alfa, beta)
436 0 : Peos_div_rho = alfa*Peos_00/d_00 + beta*Peos_m1/d_m1
437 0 : Hp_face = pow2(r_00)*Peos_div_rho/(s% cgrav(k)*s% m(k))
438 : end if
439 0 : s% Hp_face(k) = Hp_face
440 0 : s% xh(s% i_Hp, k) = Hp_face
441 0 : end subroutine set_Hp_face
442 :
443 : end subroutine remesh_for_RSP2
444 :
445 : end module hydro_rsp2_support
|
