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