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 : module eps_grav
21 :
22 : use star_private_def
23 : use const_def, only: dp, ln10
24 : use chem_def, only: chem_isos
25 : use utils_lib, only: mesa_error, is_bad
26 : use auto_diff
27 :
28 : implicit none
29 :
30 : private
31 : public :: eval_eps_grav_and_partials, zero_eps_grav_and_partials
32 :
33 : contains
34 :
35 0 : subroutine eval_eps_grav_and_partials(s, k, ierr)
36 : type (star_info), pointer :: s
37 : integer, intent(in) :: k
38 : integer, intent(out) :: ierr
39 : type(auto_diff_real_star_order1) :: eps_grav
40 : include 'formats'
41 0 : ierr = 0
42 :
43 0 : if (s% dt <= 0) then
44 0 : ierr = -1
45 0 : return
46 : end if
47 :
48 0 : call zero_eps_grav_and_partials(s, k)
49 :
50 : ! zero composition derivatives
51 : ! if include_composition_in_eps_grav is true
52 : ! then these are set in the call to eval_eps_grav_composition
53 0 : s% d_eps_grav_dx(:,k) = 0
54 :
55 0 : call eval1_eps_grav_and_partials(s, k, eps_grav, ierr)
56 0 : if (ierr /= 0) return
57 :
58 0 : s% eps_grav_ad(k) = eps_grav
59 :
60 0 : if (s% use_other_eps_grav) then
61 : ! note: call this after 1st doing the standard calculation
62 0 : call s% other_eps_grav(s% id, k, s% dt, ierr)
63 0 : if (ierr /= 0) return
64 : end if
65 :
66 : ! apply user-specified scaling factor after hook
67 0 : s% eps_grav_ad(k) = s% eps_grav_factor * s% eps_grav_ad(k)
68 :
69 : end subroutine eval_eps_grav_and_partials
70 :
71 :
72 0 : subroutine eval1_eps_grav_and_partials(s, k, eps_grav, ierr)
73 : type (star_info), pointer :: s
74 : integer, intent(in) :: k
75 : type(auto_diff_real_star_order1), intent(out) :: eps_grav
76 : integer, intent(out) :: ierr
77 :
78 : type(auto_diff_real_star_order1) :: eps_grav_lnS, eps_grav_std
79 0 : real(dp) :: alfa, Gamma
80 : logical :: using_PC
81 :
82 : include 'formats'
83 : ierr = 0
84 :
85 0 : using_PC = (s% eos_frac_PC(k) > 0)
86 :
87 0 : if (using_PC .and. s% gam_start(k) >= s% Gamma_lnS_eps_grav_full_on) then
88 0 : call do_lnS_eps_grav(s, k, eps_grav, ierr)
89 0 : else if (using_PC .and. s% gam_start(k) > s% Gamma_lnS_eps_grav_full_off) then
90 0 : Gamma = s% gam_start(k)
91 : alfa = (Gamma - s% Gamma_lnS_eps_grav_full_off) / &
92 0 : (s% Gamma_lnS_eps_grav_full_on - s% Gamma_lnS_eps_grav_full_off)
93 0 : call do_lnS_eps_grav(s, k, eps_grav_lnS, ierr)
94 0 : if (ierr /= 0) return
95 0 : call do_std_eps_grav(s, k, eps_grav_std, ierr)
96 0 : if (ierr /= 0) return
97 : ! the derivative of the blending function is missing
98 : ! but historically we've been able to get away with that
99 : ! because the two forms should match in the blend region
100 0 : eps_grav = alfa * eps_grav_lnS + (1d0 - alfa) * eps_grav_std
101 : else
102 0 : call do_std_eps_grav(s, k, eps_grav, ierr)
103 : end if
104 :
105 :
106 0 : if (ierr /= 0 .or. is_bad(eps_grav% val)) then
107 0 : ierr = -1
108 0 : s% retry_message = 'failed in eval_eps_grav_and_partials'
109 0 : if (s% report_ierr) then
110 : write(*,2) &
111 0 : 'failed in eval_eps_grav_and_partials', k, eps_grav% val
112 : end if
113 0 : if (s% stop_for_bad_nums) then
114 0 : call mesa_error(__FILE__,__LINE__,'eval1_eps_grav_and_partials')
115 : end if
116 0 : return
117 : end if
118 :
119 : end subroutine eval1_eps_grav_and_partials
120 :
121 :
122 : ! this uses the given args to calculate -T*ds/dt
123 0 : subroutine do_std_eps_grav(s, k, eps_grav, ierr)
124 : use auto_diff_support
125 : use eos_def, only: i_latent_ddlnT, i_latent_ddlnRho
126 : type (star_info), pointer :: s
127 : integer, intent(in) :: k
128 : type(auto_diff_real_star_order1), intent(out) :: eps_grav
129 : integer, intent(out) :: ierr
130 :
131 : type(auto_diff_real_star_order1) :: T, cp, grada, chiT, chiRho, dlnd_dt, dlnT_dt, &
132 : latent_ddlnT, latent_ddlnRho, eps_grav_start, eps_grav_composition_term
133 :
134 : logical :: test_partials
135 :
136 : real(dp) :: theta
137 :
138 : include 'formats'
139 0 : ierr = 0
140 :
141 : !test_partials = (k == s% solver_test_partials_k)
142 0 : test_partials = .false.
143 :
144 : ! select time-centering
145 0 : if (s% use_time_centered_eps_grav .and. .not. s% doing_relax) then
146 0 : theta = 0.5_dp
147 : else
148 0 : theta = 1.0_dp
149 : end if
150 :
151 : ! total time derivatives of (Rho, T) basis
152 0 : dlnd_dt = wrap_dxh_lnd(s,k)/s% dt
153 0 : dlnT_dt = wrap_dxh_lnT(s,k)/s% dt
154 :
155 : ! thermo quantities
156 0 : T = wrap_T_00(s, k)
157 0 : Cp = wrap_Cp_00(s, k)
158 0 : grada = wrap_grad_ad_00(s, k)
159 0 : chiT = wrap_chiT_00(s, k)
160 0 : chiRho = wrap_chiRho_00(s, k)
161 :
162 : ! MESA I, Equation (12)
163 0 : eps_grav = -T*Cp * ((1d0 - grada*chiT)*dlnT_dt - grada*chiRho*dlnd_dt)
164 :
165 : ! phase transition latent heat
166 : ! Jermyn et al. (2021), Equation (47)
167 0 : latent_ddlnRho = wrap_latent_ddlnRho_00(s,k)
168 0 : latent_ddlnT = wrap_latent_ddlnT_00(s,k)
169 0 : eps_grav = eps_grav - (dlnd_dt * latent_ddlnRho + dlnT_dt * latent_ddlnT)
170 :
171 :
172 : ! for time centered version
173 0 : if (s% use_time_centered_eps_grav .and. .not. s% doing_relax) then
174 :
175 : ! start values are constants during Newton iters
176 : eps_grav_start = -s% T_start(k)*s% cp_start(k) &
177 0 : * ((1d0 - s% grada_start(k)*s% chiT_start(k))*dlnT_dt - s% grada_start(k)*s% chiRho_start(k)*dlnd_dt)
178 :
179 : ! phase transition latent heat
180 0 : eps_grav_start = eps_grav_start - (dlnd_dt * s% latent_ddlnRho_start(k) + dlnT_dt * s% latent_ddlnT_start(k))
181 :
182 0 : eps_grav = theta * eps_grav + (1d0-theta) * eps_grav_start
183 :
184 : end if
185 :
186 :
187 0 : if (s% include_composition_in_eps_grav) then
188 0 : call eval_eps_grav_composition(s, k, eps_grav_composition_term, ierr)
189 0 : if (ierr /= 0) return
190 0 : eps_grav = eps_grav + eps_grav_composition_term
191 : end if
192 :
193 0 : if (is_bad(eps_grav% val)) then
194 0 : ierr = -1
195 0 : s% retry_message = 'do_lnd_eps_grav -- bad value for eps_grav'
196 0 : if (s% report_ierr) &
197 0 : write(*,2) 'do_lnd_eps_grav -- bad value for eps_grav', k, eps_grav% val
198 0 : if (s% stop_for_bad_nums) call mesa_error(__FILE__,__LINE__,'do_lnd_eps_grav')
199 0 : return
200 : end if
201 :
202 : if (test_partials) then
203 : s% solver_test_partials_val = 0
204 : s% solver_test_partials_var = 0
205 : s% solver_test_partials_dval_dx = 0
206 : write(*,*) 'do_std_eps_grav chiT', s% solver_test_partials_var
207 : end if
208 :
209 0 : end subroutine do_std_eps_grav
210 :
211 :
212 0 : subroutine do_lnS_eps_grav(s, k, eps_grav, ierr)
213 0 : use auto_diff_support
214 : type (star_info), pointer :: s
215 : integer, intent(in) :: k
216 : type(auto_diff_real_star_order1), intent(out) :: eps_grav
217 : integer, intent(out) :: ierr
218 :
219 0 : real(dp) :: entropy_start
220 : type(auto_diff_real_star_order1) :: entropy, T
221 :
222 : include 'formats'
223 0 : ierr = 0
224 :
225 0 : T = wrap_T_00(s,k)
226 0 : entropy = wrap_s_00(s,k)
227 0 : entropy_start = exp(s% lnS_start(k))
228 :
229 0 : eps_grav = -T*(entropy - entropy_start)/s% dt
230 :
231 : ! NOTE: the correct version of the composition term to go with TdS is
232 : ! -sum_i (\partial e/\partial Y_i)_{s,\rho} dY_i
233 : ! see for example, MESA IV, equation (59)
234 : !
235 : ! When on PC near crystallization (the only place the lnS form is used)
236 : ! there are typically no composition changes, so we drop this term
237 : !
238 : ! The term we normally use comes from expanding in the (rho,T,X) basis
239 : ! -sum_ (\partial e/\partial X_i)_{T,\rho} dX_i
240 : ! and in practice can be a reasonable approximation to the above.
241 : !
242 : ! If an such approximation is desired, one could use the following code:
243 :
244 : ! if (s% include_composition_in_eps_grav) then
245 : ! call eval_eps_grav_composition(s, k, eps_grav_composition_term, ierr)
246 : ! if (ierr /= 0) return
247 : ! eps_grav = eps_grav + eps_grav_composition_term
248 : ! end if
249 :
250 0 : if (is_bad(eps_grav% val)) then
251 0 : ierr = -1
252 0 : s% retry_message = 'do_lnS_eps_grav -- bad value for eps_grav'
253 0 : if (s% report_ierr) &
254 0 : write(*,2) 'do_lnS_eps_grav -- bad value for eps_grav', k, eps_grav% val
255 0 : if (s% stop_for_bad_nums) call mesa_error(__FILE__,__LINE__,'do_lnS_eps_grav')
256 : return
257 : end if
258 :
259 0 : end subroutine do_lnS_eps_grav
260 :
261 :
262 0 : subroutine eval_eps_grav_composition(s, k, eps_grav_composition_term, ierr)
263 0 : use auto_diff_support, only: wrap
264 : use eos_support, only: get_eos
265 : use eos_def, only: num_eos_basic_results, num_eos_d_dxa_results, i_lnE, i_lnPgas
266 :
267 : type (star_info), pointer :: s
268 : integer, intent(in) :: k
269 : type(auto_diff_real_star_order1), intent(out) :: eps_grav_composition_term
270 : integer, intent(out) :: ierr
271 : real(dp) :: &
272 0 : e, e_start, de, d_de_dlnd, d_de_dlnT, &
273 0 : e_with_xa_start, d_e_with_xa_start_dlnd, d_e_with_xa_start_dlnT, &
274 0 : e_with_DT_start, Pgas_with_DT_start
275 : real(dp), dimension(num_eos_basic_results) :: &
276 0 : res, dres_dlnd, dres_dlnT
277 0 : real(dp) :: dres_dxa(num_eos_d_dxa_results,s% species)
278 : integer :: j
279 : logical :: test_partials
280 :
281 0 : real(dp) :: theta
282 :
283 : include 'formats'
284 0 : ierr = 0
285 :
286 0 : eps_grav_composition_term = 0
287 :
288 0 : if (s% use_time_centered_eps_grav .and. .not. s% doing_relax) then
289 : theta = 0.5_dp
290 : else
291 0 : theta = 1.0_dp
292 : end if
293 :
294 : ! some EOSes have composition partials and some do not
295 : ! those currently without dx partials are PC & Skye
296 : ! however, composition partials of lnE & lnP were revised by the call to
297 : ! fix_d_eos_dxa_partials at the beginning of eval_equ_for_solver
298 :
299 : ! directly fill array with composition derivatives
300 : ! we only have lnE and lnP derivs, so use
301 : ! eps_grav = -de/dt + (P/rho) * dlnd/dt
302 0 : do j=1, s% species
303 : s% d_eps_grav_dx(j,k) = -s% energy(k) * s% dlnE_dxa_for_partials(j,k)/s% dt + &
304 0 : (s% Peos(k) / s% Rho(k)) * s% dlnPeos_dxa_for_partials(j,k) * s% dxh_lnd(k)/s% dt
305 : end do
306 :
307 0 : e = s% energy(k)
308 : call get_eos( &
309 : s, k, s% xa_start(:,k), &
310 : s% rho(k), s% lnd(k)/ln10, s% T(k), s% lnT(k)/ln10, &
311 : res, dres_dlnd, dres_dlnT, &
312 0 : dres_dxa, ierr)
313 0 : if (ierr /= 0) then
314 0 : if (s% report_ierr) write(*,2) 'failed in get_eos with xa_start', k
315 0 : return
316 : end if
317 :
318 0 : e_with_xa_start = exp(res(i_lnE))
319 0 : d_e_with_xa_start_dlnd = dres_dlnd(i_lnE)*e_with_xa_start
320 0 : d_e_with_xa_start_dlnT = dres_dlnT(i_lnE)*e_with_xa_start
321 :
322 0 : de = (e - e_with_xa_start)
323 0 : d_de_dlnd = (s% dE_dRho_for_partials(k)*s% Rho(k) - d_e_with_xa_start_dlnd)
324 0 : d_de_dlnT = (s% Cv_for_partials(k)*s% T(k) - d_e_with_xa_start_dlnT)
325 :
326 0 : if (s% use_time_centered_eps_grav .and. .not. s% doing_relax) then
327 :
328 0 : e_start = s% energy_start(k)
329 :
330 : call get_eos( &
331 : s, k, s% xa(:,k), &
332 : s% rho_start(k), s% lnd_start(k)/ln10, s% T_start(k), s% lnT_start(k)/ln10, &
333 : res, dres_dlnd, dres_dlnT, &
334 0 : dres_dxa, ierr)
335 0 : if (ierr /= 0) then
336 0 : if (s% report_ierr) write(*,2) 'failed in get_eos with xa_start', k
337 0 : return
338 : end if
339 :
340 0 : e_with_DT_start = exp(res(i_lnE))
341 0 : de = theta * de + (1d0 - theta) * (e_with_DT_start - e_start)
342 0 : d_de_dlnd = theta * d_de_dlnd
343 0 : d_de_dlnT = theta * d_de_dlnT
344 :
345 : Pgas_with_DT_start = exp(res(i_lnPgas))
346 :
347 : ! combine with end-of-step composition derivatives
348 :
349 : ! until EOSes always provide composition derivatives, ignore this
350 : ! the end-of-step term should generally be similar enough to be OK
351 :
352 : ! do j=1, s% species
353 : ! s% d_eps_grav_dx(j,k) = theta * s% d_eps_grav_dx(j,k) + (1d0 - theta) * &
354 : ! (-e_with_DT_start*d_dxa(i_lnE,j)/s% dt + &
355 : ! Pgas_with_DT_start*d_dxa(i_lnPgas,j)/s% rho_start(k)*s% dxh_lnd(k)/s% dt)
356 : ! end do
357 :
358 : end if
359 :
360 : call wrap(eps_grav_composition_term, -de/s% dt, &
361 : 0d0, -d_de_dlnd/s% dt, 0d0, &
362 : 0d0, -d_de_dlnT/s% dt, 0d0, &
363 : 0d0, 0d0, 0d0, &
364 : 0d0, 0d0, 0d0, &
365 : 0d0, 0d0, 0d0, &
366 : 0d0, 0d0, 0d0, &
367 : 0d0, 0d0, 0d0, &
368 : 0d0, 0d0, 0d0, &
369 : 0d0, 0d0, 0d0, &
370 : 0d0, 0d0, 0d0, &
371 0 : 0d0, 0d0, 0d0)
372 :
373 : ! add easy access to this quantity in star
374 0 : s% eps_grav_composition_term(k) = eps_grav_composition_term% val
375 :
376 : !test_partials = (k == s% solver_test_partials_k)
377 0 : test_partials = .false.
378 :
379 : if (test_partials) then
380 : s% solver_test_partials_val = de
381 : s% solver_test_partials_var = s% i_lnT
382 : s% solver_test_partials_dval_dx = d_de_dlnT
383 : write(*,*) 'get_dedt', s% solver_test_partials_var
384 : end if
385 :
386 0 : if (is_bad(eps_grav_composition_term% val)) then
387 0 : if (s% report_ierr) then
388 0 : write(*, *) s% retry_message
389 0 : write(*,2) 'eps_grav_composition_term', k, eps_grav_composition_term% val
390 : end if
391 0 : if (s% stop_for_bad_nums) then
392 0 : write(*,2) 'include_composition_in_eps_grav -- bad value for eps_grav_composition_term', &
393 0 : k, eps_grav_composition_term% val
394 0 : call mesa_error(__FILE__,__LINE__,'eval_eps_grav_composition')
395 : end if
396 0 : return
397 : end if
398 :
399 0 : end subroutine eval_eps_grav_composition
400 :
401 :
402 52348 : subroutine zero_eps_grav_and_partials(s, k)
403 : type (star_info), pointer :: s
404 : integer, intent(in) :: k
405 52348 : s% eps_grav_ad(k) = 0
406 0 : end subroutine zero_eps_grav_and_partials
407 :
408 : end module eps_grav
|
