Line data Source code
1 : ! ***********************************************************************
2 : !
3 : ! Copyright (C) 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 eoscms_eval
21 :
22 : use eos_def
23 : use num_lib
24 : use const_def, only: avo, crad, ln10, arg_not_provided, mp, kerg, dp, qp, mesa_data_dir, one_third
25 : use utils_lib, only: is_bad, mesa_error
26 : use math_lib
27 : use interp_2d_lib_db
28 :
29 : implicit none
30 :
31 : logical, parameter :: CMS_cubic_in_X = .false.
32 :
33 : integer, parameter :: CMS_num_Xs = 11
34 : integer, parameter :: min_for_cubic = 2
35 : integer, parameter :: max_for_cubic = CMS_num_Xs - 2
36 :
37 : character(len=3) :: CMS_Xstr(CMS_num_Xs) = ['000','010','020','030','040','050','060','070','080','090','100']
38 :
39 : real(dp) :: CMS_Xvals(CMS_num_Xs) = [ 0.0_dp, 0.1_dp, 0.2_dp, 0.3_dp, 0.4_dp, 0.5_dp, 0.6_dp, 0.7_dp, 0.8_dp, 0.9_dp, 1.0_dp]
40 :
41 : !modeled on EosDT_XZ_Info from eos_def
42 : type eosCMS_X_Info
43 : real(dp) :: logRho_min
44 : real(dp) :: logRho_max
45 : real(dp) :: delta_logRho
46 : integer :: num_logRhos
47 : real(dp) :: logT_min
48 : real(dp) :: logT_max
49 : real(dp) :: delta_logT
50 : integer :: num_logTs
51 : real(dp), pointer :: logRhos(:), logTs(:)
52 : real(dp), pointer :: tbl1(:)
53 : integer :: version
54 : end type eosCMS_X_Info
55 :
56 : type(eosCMS_X_Info), target :: eosCMS_X_data(CMS_num_Xs)
57 : logical :: eosCMS_X_loaded(CMS_num_Xs) = .false.
58 :
59 : contains
60 :
61 0 : subroutine eosCMS_init(ierr)
62 : integer, intent(out) :: ierr
63 0 : eosCMS_X_loaded = .false.
64 0 : ierr=0
65 0 : end subroutine eosCMS_init
66 :
67 :
68 0 : subroutine Get_CMS_alfa( &
69 : rq, logRho, logT, Z, abar, zbar, &
70 : alfa, d_alfa_dlogT, d_alfa_dlogRho, &
71 : ierr)
72 : use auto_diff
73 : type (EoS_General_Info), pointer :: rq
74 : real(dp), intent(in) :: logRho, logT, Z, abar, zbar
75 : real(dp), intent(out) :: alfa, d_alfa_dlogT, d_alfa_dlogRho
76 : integer, intent(out) :: ierr
77 :
78 : type(auto_diff_real_2var_order1) :: logT_auto, logRho_auto, logQ_auto
79 : type(auto_diff_real_2var_order1) :: blend, blend_logT, blend_logRho, blend_logQ
80 :
81 : include 'formats'
82 :
83 0 : ierr = 0
84 :
85 : ! logRho is val1
86 0 : logRho_auto% val = logRho
87 0 : logRho_auto% d1val1 = 1d0
88 0 : logRho_auto% d1val2 = 0d0
89 :
90 : ! logT is val2
91 0 : logT_auto% val = logT
92 0 : logT_auto% d1val1 = 0d0
93 0 : logT_auto% d1val2 = 1d0
94 :
95 0 : logQ_auto = logRho_auto - 2d0*logT_auto + 12d0
96 :
97 : ! for blend variables 1 is CMS, 0 is other
98 : ! (this is the opposite of the final alfa)
99 :
100 : ! logT blend
101 0 : if (logT_auto < rq% logT_min_for_any_CMS) then
102 0 : blend_logT = 0d0
103 0 : else if (logT_auto < rq% logT_min_for_all_CMS) then
104 0 : blend_logT = (logT_auto - rQ% logT_min_for_any_CMS) / (rq% logT_min_for_all_CMS - rq% logT_min_for_any_CMS)
105 0 : else if (logT_auto < rq% logT_max_for_all_CMS) then
106 0 : blend_logT = 1d0
107 0 : else if (logT_auto < rq% logT_max_for_any_CMS) then
108 0 : blend_logT = (logT_auto - rQ% logT_max_for_any_CMS) / (rq% logT_max_for_all_CMS - rq% logT_max_for_any_CMS)
109 : else
110 0 : blend_logT = 0
111 : end if
112 :
113 :
114 : ! logRho blend
115 0 : if (logRho_auto < rq% logRho_min_for_any_CMS) then
116 0 : blend_logRho = 0d0
117 0 : else if (logRho_auto < rq% logRho_min_for_all_CMS) then
118 0 : blend_logRho = (logRho_auto - rQ% logRho_min_for_any_CMS) / (rq% logRho_min_for_all_CMS - rq% logRho_min_for_any_CMS)
119 0 : else if (logRho_auto < rq% logRho_max_for_all_CMS) then
120 0 : blend_logRho = 1d0
121 0 : else if (logRho_auto < rq% logRho_max_for_any_CMS) then
122 0 : blend_logRho = (logRho_auto - rQ% logRho_max_for_any_CMS) / (rq% logRho_max_for_all_CMS - rq% logRho_max_for_any_CMS)
123 : else
124 0 : blend_logRho = 0
125 : end if
126 :
127 :
128 : ! logQ blend
129 0 : if (logQ_auto < rq% logQ_min_for_any_CMS) then
130 0 : blend_logQ = 0d0
131 0 : else if (logQ_auto < rq% logQ_min_for_all_CMS) then
132 0 : blend_logQ = (logQ_auto - rQ% logQ_min_for_any_CMS) / (rq% logQ_min_for_all_CMS - rq% logQ_min_for_any_CMS)
133 0 : else if (logQ_auto < rq% logQ_max_for_all_CMS) then
134 0 : blend_logQ = 1d0
135 0 : else if (logQ_auto < rq% logQ_max_for_any_CMS) then
136 0 : blend_logQ = (logQ_auto - rQ% logQ_max_for_any_CMS) / (rq% logQ_max_for_all_CMS - rq% logQ_max_for_any_CMS)
137 : else
138 0 : blend_logQ = 0
139 : end if
140 :
141 : ! combine blends
142 0 : blend = blend_logRho * blend_logT * blend_logQ
143 :
144 0 : alfa = 1d0 - blend% val
145 0 : d_alfa_dlogRho = -blend% d1val1
146 0 : d_alfa_dlogT = -blend% d1val2
147 :
148 0 : end subroutine Get_CMS_alfa
149 :
150 :
151 0 : subroutine get_CMS_for_eosdt( &
152 : handle, dbg, Z, X, abar, zbar, &
153 0 : species, chem_id, net_iso, xa, &
154 : rho, logRho, T, logT, remaining_fraction, &
155 0 : res, d_dlnd, d_dlnT, d_dxa, &
156 : skip, ierr)
157 0 : use chem_def, only: chem_isos
158 : use interp_1d_lib, only: interp_4pt_pm
159 : integer, intent(in) :: handle
160 : logical, intent(in) :: dbg
161 : real(dp), intent(in) :: Z, X, abar, zbar, remaining_fraction
162 : integer, intent(in) :: species
163 : integer, pointer :: chem_id(:), net_iso(:)
164 : real(dp), intent(in) :: xa(:)
165 : real(dp), intent(in) :: rho, logRho, T, logT
166 : real(dp), intent(inout), dimension(nv) :: &
167 : res, d_dlnd, d_dlnT
168 : real(dp), intent(inout), dimension(nv, species) :: d_dxa
169 : logical, intent(out) :: skip
170 : integer, intent(out) :: ierr
171 : integer :: iX, i
172 : type (EoS_General_Info), pointer :: rq
173 0 : real(dp) :: res1(nv), res2(nv), res3(nv), res4(nv)
174 0 : real(dp) :: dres1_dlnT(nv), dres2_dlnT(nv), dres3_dlnT(nv), dres4_dlnT(nv)
175 0 : real(dp) :: dres1_dlnRho(nv), dres2_dlnRho(nv), dres3_dlnRho(nv), dres4_dlnRho(nv)
176 0 : real(dp) :: d_dX(nv)
177 0 : real(dp) :: alfa, beta, dbeta_dX, dalfa_dX, xx(4), y(4), a(3), dX
178 0 : rq => eos_handles(handle)
179 :
180 0 : if(rq% CMS_use_fixed_composition)then !do fixed composition (one table only)
181 :
182 0 : if(rq% CMS_fixed_composition_index < 0 .or. rq% CMS_fixed_composition_index > 10)then
183 0 : write(*,*) 'invalid value for CMS_fixed_composition_index. See eos.defaults.'
184 0 : ierr=-1
185 0 : return
186 : end if
187 :
188 0 : iX = rq% CMS_fixed_composition_index + 1
189 :
190 0 : call eval_eosCMS_fixed_X(iX,logRho,logT,res,d_dlnT,d_dlnd,ierr)
191 :
192 0 : if(ierr/=0) then
193 0 : write(*,*) 'failed in get_CMS_for_eosdt'
194 0 : return
195 : end if
196 :
197 : ! composition derivatives; here composition is constant so no change
198 0 : d_dxa = 0
199 :
200 :
201 : else !do full composition
202 : !locate X values in the tables such that Xvals(iX) <= X < Xvals(iX+1)
203 0 : if (X <= CMS_Xvals(1)) then
204 0 : iX = 1
205 0 : if (X < 0) write(*,*) 'warning: X < 0 in eosCMS'
206 0 : else if (X >= CMS_Xvals(CMS_num_Xs-1)) then
207 0 : iX = CMS_num_Xs-1
208 0 : if (X > 1) write(*,*) 'warning: X > 1 in eosCMS'
209 : else
210 0 : do i = 2, CMS_num_Xs-1
211 0 : if (X < CMS_Xvals(i) )then
212 0 : iX = i-1; exit
213 : end if
214 : end do
215 : end if
216 :
217 : !interpolation always bicubic in logRho and logT
218 0 : if(CMS_cubic_in_X .and. iX > 2 .and. iX < CMS_num_Xs - 2)then
219 : !do cubic interpolation in X
220 : call eval_eosCMS_fixed_X(iX-1,logRho,logT,res1,dres1_dlnT,dres1_dlnRho,ierr)
221 : call eval_eosCMS_fixed_X(iX ,logRho,logT,res2,dres2_dlnT,dres2_dlnRho,ierr)
222 : call eval_eosCMS_fixed_X(iX+1,logRho,logT,res3,dres3_dlnT,dres3_dlnRho,ierr)
223 : call eval_eosCMS_fixed_X(iX+2,logRho,logT,res4,dres4_dlnT,dres4_dlnRho,ierr)
224 : if(ierr/=0) then
225 : write(*,*) 'failed in get_CMS_for_eosdt'
226 : return
227 : end if
228 : XX(1:4) = CMS_Xvals(iX-1:iX+2)
229 : dX=X-CMS_Xvals(iX) ! assumes fixed dX spacing
230 : do i=1,nv
231 : !result
232 : y(1:4) = [res1(i), res2(i), res3(i), res4(i)]
233 : call interp_4pt_pm(XX,y,a)
234 : res(i) = y(2) + dX*(a(1) + dX*(a(2) + dX*a(3)))
235 : !dres/dX
236 : d_dX(i) = a(1) + dX*(2*a(2) + 3*dX*a(3))
237 : !dres/dlnT
238 : y(1:4) = [dres1_dlnT(i), dres2_dlnT(i), dres3_dlnT(i), dres4_dlnT(i)]
239 : call interp_4pt_pm(XX,y,a)
240 : d_dlnT(i) = y(2) + dX*(a(1) + dX*(a(2) + dX*a(3)))
241 : !dres/dlnRho
242 : y(1:4) = [dres1_dlnRho(i), dres2_dlnRho(i), dres3_dlnRho(i), dres4_dlnRho(i)]
243 : call interp_4pt_pm(XX,y,a)
244 : d_dlnd(i) = y(2) + dX*(a(1) + dX*(a(2) + dX*a(3)))
245 : end do
246 : else !linear interpolation in X
247 0 : call eval_eosCMS_fixed_X(iX ,logRho,logT,res1,dres1_dlnT,dres1_dlnRho,ierr)
248 0 : call eval_eosCMS_fixed_X(iX+1,logRho,logT,res2,dres2_dlnT,dres2_dlnRho,ierr)
249 0 : if(ierr/=0) then
250 0 : write(*,*) 'failed in get_CMS_for_eosdt'
251 0 : return
252 : end if
253 0 : beta = (X-CMS_Xvals(iX))/(CMS_Xvals(iX+1)-CMS_Xvals(iX))
254 0 : alfa = 1._dp - beta
255 0 : dbeta_dX = 1d0/(CMS_Xvals(iX+1)-CMS_Xvals(iX))
256 0 : dalfa_dX = -dbeta_dX
257 0 : res = alfa*res1 + beta*res2
258 0 : d_dX = dalfa_dX*res1 + dbeta_dX*res2
259 0 : d_dlnT = alfa*dres1_dlnT + beta*dres2_dlnT
260 0 : d_dlnd = alfa*dres1_dlnRho + beta*dres2_dlnRho
261 : end if
262 :
263 : ! composition derivatives
264 0 : do i=1,species
265 0 : select case(chem_isos% Z(chem_id(i))) ! charge
266 : case (1) ! X
267 0 : d_dxa(:,i) = d_dX
268 : case (2) ! Y
269 0 : d_dxa(:,i) = 0
270 : case default ! Z
271 0 : d_dxa(:,i) = 0
272 : end select
273 : end do
274 :
275 : end if
276 :
277 0 : skip = .false.
278 :
279 : ! CMS tables do not include radiation. Add it.
280 0 : if (rq% include_radiation) call include_radiation(Z, X, abar, zbar, &
281 : species, chem_id, net_iso, xa, &
282 : rho, logRho, T, logT, &
283 : res, d_dlnd, d_dlnT, d_dxa, &
284 0 : ierr)
285 :
286 : ! zero phase information
287 0 : res(i_phase:i_latent_ddlnRho) = 0d0
288 0 : d_dlnT(i_phase:i_latent_ddlnRho) = 0d0
289 0 : d_dlnd(i_phase:i_latent_ddlnRho) = 0d0
290 :
291 : ! zero all components
292 0 : res(i_frac:i_frac+num_eos_frac_results-1) = 0.0d0
293 0 : d_dlnd(i_frac:i_frac+num_eos_frac_results-1) = 0.0d0
294 0 : d_dlnT(i_frac:i_frac+num_eos_frac_results-1) = 0.0d0
295 :
296 : ! mark this one
297 0 : res(i_frac_CMS) = 1.0d0
298 :
299 0 : end subroutine get_CMS_for_eosdt
300 :
301 :
302 0 : subroutine include_radiation(Z, X, abar, zbar, &
303 : species, chem_id, net_iso, xa, &
304 : rho_in, logRho, T_in, logT, &
305 : res, d_dlnd, d_dlnT, d_dxa, &
306 : ierr)
307 0 : use const_def, only: dp
308 : use auto_diff
309 : real(dp), intent(in) :: Z, X, abar, zbar
310 : integer, intent(in) :: species
311 : integer, pointer :: chem_id(:), net_iso(:)
312 : real(dp), intent(in) :: xa(:)
313 : real(dp), intent(in) :: rho_in, logRho, T_in, logT
314 : real(dp), intent(inout), dimension(nv) :: &
315 : res, d_dlnd, d_dlnT
316 : real(dp), intent(inout), dimension(nv, species) :: d_dxa
317 : integer, intent(out) :: ierr
318 :
319 : type(auto_diff_real_2var_order1) :: T, Rho, P, Pgas, Prad, lnPgas, lnE, lnS, &
320 : grad_ad, chiRho, chiT, Cp, Cv, dE_dRho, dS_dT, dS_dRho, mu, lnfree_e, gamma1, gamma3, eta
321 :
322 0 : ierr = 0
323 :
324 0 : T% val = T_in
325 0 : T% d1val1 = 1d0
326 0 : T% d1val2 = 0d0
327 :
328 0 : Rho% val = Rho_in
329 0 : Rho% d1val1 = 0d0
330 0 : Rho% d1val2 = 1d0
331 :
332 :
333 : ! unpack results into auto_diff types
334 0 : lnPgas % val = res(i_lnPgas)
335 0 : lnE % val = res(i_lnE)
336 0 : lnS % val = res(i_lnS)
337 : grad_ad % val = res(i_grad_ad)
338 0 : chiRho % val = res(i_chiRho)
339 0 : chiT % val = res(i_chiT)
340 : Cp % val = res(i_Cp)
341 0 : Cv % val = res(i_Cv)
342 : dE_dRho % val = res(i_dE_dRho)
343 : dS_dT % val = res(i_dS_dT)
344 : dS_dRho % val = res(i_dS_dRho)
345 0 : mu % val = res(i_mu)
346 0 : lnfree_e % val = res(i_lnfree_e)
347 : gamma1 % val = res(i_gamma1)
348 : gamma3 % val = res(i_gamma3)
349 0 : eta % val = res(i_eta)
350 :
351 0 : lnPgas % d1val1 = d_dlnT(i_lnPgas) / T% val
352 0 : lnE % d1val1 = d_dlnT(i_lnE) / T% val
353 0 : lnS % d1val1 = d_dlnT(i_lnS) / T% val
354 : grad_ad % d1val1 = d_dlnT(i_grad_ad) / T% val
355 0 : chiRho % d1val1 = d_dlnT(i_chiRho) / T% val
356 0 : chiT % d1val1 = d_dlnT(i_chiT) / T% val
357 : Cp % d1val1 = d_dlnT(i_Cp) / T% val
358 0 : Cv % d1val1 = d_dlnT(i_Cv) / T% val
359 : dE_dRho % d1val1 = d_dlnT(i_dE_dRho) / T% val
360 : dS_dT % d1val1 = d_dlnT(i_dS_dT) / T% val
361 : dS_dRho % d1val1 = d_dlnT(i_dS_dRho) / T% val
362 0 : mu % d1val1 = d_dlnT(i_mu) / T% val
363 0 : lnfree_e % d1val1 = d_dlnT(i_lnfree_e) / T% val
364 : gamma1 % d1val1 = d_dlnT(i_gamma1) / T% val
365 : gamma3 % d1val1 = d_dlnT(i_gamma3) / T% val
366 0 : eta % d1val1 = d_dlnT(i_eta) / T% val
367 :
368 0 : lnPgas % d1val2 = d_dlnd(i_lnPgas) / Rho% val
369 0 : lnE % d1val2 = d_dlnd(i_lnE) / Rho% val
370 0 : lnS % d1val2 = d_dlnd(i_lnS) / Rho% val
371 : grad_ad % d1val2 = d_dlnd(i_grad_ad) / Rho% val
372 0 : chiRho % d1val2 = d_dlnd(i_chiRho) / Rho% val
373 0 : chiT % d1val2 = d_dlnd(i_chiT) / Rho% val
374 : Cp % d1val2 = d_dlnd(i_Cp) / Rho% val
375 0 : Cv % d1val2 = d_dlnd(i_Cv) / Rho% val
376 : dE_dRho % d1val2 = d_dlnd(i_dE_dRho) / Rho% val
377 : dS_dT % d1val2 = d_dlnd(i_dS_dT) / Rho% val
378 : dS_dRho % d1val2 = d_dlnd(i_dS_dRho) / Rho% val
379 0 : mu % d1val2 = d_dlnd(i_mu) / Rho% val
380 0 : lnfree_e % d1val2 = d_dlnd(i_lnfree_e) / Rho% val
381 : gamma1 % d1val2 = d_dlnd(i_gamma1) / Rho% val
382 : gamma3 % d1val2 = d_dlnd(i_gamma3) / Rho% val
383 0 : eta % d1val2 = d_dlnd(i_eta) / Rho% val
384 :
385 :
386 : ! add radiation
387 0 : Pgas = exp(lnPgas)
388 0 : Prad = one_third*crad*pow4(T)
389 0 : P = Pgas + Prad
390 0 : lnE = log(exp(lnE) + 3d0*Prad/Rho)
391 0 : lnS = log(exp(lnS) + 4d0*Prad/(Rho*T))
392 0 : Cv = Cv + 12d0*Prad/(Rho*T)
393 0 : chiT = chiT * Pgas/P + 4d0*Prad/P
394 0 : chiRho = chiRho * Pgas/P
395 0 : gamma3 = 1d0 + P/Rho * chiT/(T*Cv)
396 0 : gamma1 = chiT*(gamma3-1d0) + chiRho
397 0 : grad_ad = (gamma3-1d0)/gamma1
398 0 : Cp = Cv * gamma1/chiRho
399 0 : dE_dRho = (1d0-chiT)*P/(Rho*Rho)
400 0 : dS_dT = Cv/T
401 0 : dS_dRho = -P*chiT/(Rho*Rho*T)
402 :
403 :
404 : ! repack results
405 0 : res(i_lnPgas) = lnPgas % val
406 0 : res(i_lnE) = lnE % val
407 0 : res(i_lnS) = lnS % val
408 0 : res(i_grad_ad) = grad_ad % val
409 0 : res(i_chiRho) = chiRho % val
410 0 : res(i_chiT) = chiT % val
411 0 : res(i_Cp) = Cp % val
412 0 : res(i_Cv) = Cv % val
413 0 : res(i_dE_dRho) = dE_dRho % val
414 0 : res(i_dS_dT) = dS_dT % val
415 0 : res(i_dS_dRho) = dS_dRho % val
416 0 : res(i_mu) = mu % val
417 0 : res(i_lnfree_e) = lnfree_e % val
418 0 : res(i_gamma1) = gamma1 % val
419 0 : res(i_gamma3) = gamma3 % val
420 0 : res(i_eta) = eta % val
421 :
422 0 : d_dlnT(i_lnPgas) = lnPgas % d1val1 * T % val
423 0 : d_dlnT(i_lnE) = lnE % d1val1 * T % val
424 0 : d_dlnT(i_lnS) = lnS % d1val1 * T % val
425 0 : d_dlnT(i_grad_ad) = grad_ad % d1val1 * T % val
426 0 : d_dlnT(i_chiRho) = chiRho % d1val1 * T % val
427 0 : d_dlnT(i_chiT) = chiT % d1val1 * T % val
428 0 : d_dlnT(i_Cp) = Cp % d1val1 * T % val
429 0 : d_dlnT(i_Cv) = Cv % d1val1 * T % val
430 0 : d_dlnT(i_dE_dRho) = dE_dRho % d1val1 * T % val
431 0 : d_dlnT(i_dS_dT) = dS_dT % d1val1 * T % val
432 0 : d_dlnT(i_dS_dRho) = dS_dRho % d1val1 * T % val
433 0 : d_dlnT(i_mu) = mu % d1val1 * T % val
434 0 : d_dlnT(i_lnfree_e) = lnfree_e % d1val1 * T % val
435 0 : d_dlnT(i_gamma1) = gamma1 % d1val1 * T % val
436 0 : d_dlnT(i_gamma3) = gamma3 % d1val1 * T % val
437 0 : d_dlnT(i_eta) = eta % d1val1 * T % val
438 :
439 0 : d_dlnd(i_lnPgas) = lnPgas % d1val2 * RHO % val
440 0 : d_dlnd(i_lnE) = lnE % d1val2 * RHO % val
441 0 : d_dlnd(i_lnS) = lnS % d1val2 * RHO % val
442 0 : d_dlnd(i_grad_ad) = grad_ad % d1val2 * RHO % val
443 0 : d_dlnd(i_chiRho) = chiRho % d1val2 * RHO % val
444 0 : d_dlnd(i_chiT) = chiT % d1val2 * RHO % val
445 0 : d_dlnd(i_Cp) = Cp % d1val2 * RHO % val
446 0 : d_dlnd(i_Cv) = Cv % d1val2 * RHO % val
447 0 : d_dlnd(i_dE_dRho) = dE_dRho % d1val2 * RHO % val
448 0 : d_dlnd(i_dS_dT) = dS_dT % d1val2 * RHO % val
449 0 : d_dlnd(i_dS_dRho) = dS_dRho % d1val2 * RHO % val
450 0 : d_dlnd(i_mu) = mu % d1val2 * RHO % val
451 0 : d_dlnd(i_lnfree_e) = lnfree_e % d1val2 * RHO % val
452 0 : d_dlnd(i_gamma1) = gamma1 % d1val2 * RHO % val
453 0 : d_dlnd(i_gamma3) = gamma3 % d1val2 * RHO % val
454 0 : d_dlnd(i_eta) = eta % d1val2 * RHO % val
455 :
456 0 : end subroutine include_radiation
457 :
458 :
459 0 : subroutine eval_eosCMS_fixed_X(iX,logRho,logT,res,dres_dlnT,dres_dlnRho,ierr)
460 0 : use eosdt_support, only: Do_EoS_Interpolations
461 : integer, intent(in) :: iX
462 : real(dp), intent(in) :: logRho, logT
463 : real(dp), intent(out) :: res(nv), dres_dlnT(nv), dres_dlnRho(nv)
464 : integer, intent(out) :: ierr
465 : type(eosCMS_X_info), pointer :: c
466 : integer :: iRho, iT
467 0 : real(dp) :: fval(nv), df_dx(nv), df_dy(nv)
468 0 : real(dp) :: logT0, logRho0, logT1, logRho1, my_logT, my_logRho
469 0 : ierr = 0
470 0 : !$OMP CRITICAL(OMP_CRITICAL_IX)
471 0 : if(.not.eosCMS_X_loaded(iX)) call load_eosCMS_table(iX,ierr)
472 : !$OMP END CRITICAL(OMP_CRITICAL_IX)
473 :
474 0 : my_logT = logT
475 0 : my_logRho = logRho
476 :
477 0 : c => eosCMS_X_data(iX)
478 :
479 0 : call locate_logRho(c, my_logRho, iRho, logRho0, logRho1)
480 0 : call locate_logT (c, my_logT, iT, logT0, logT1)
481 :
482 : call Do_EoS_Interpolations(1, nv, nv, &
483 : c% num_logTs, c% logTs, c% num_logRhos, c% logRhos, &
484 : c% tbl1, iT, iRho, logT0, my_logT, logT1, &
485 0 : logRho0, my_logRho, logRho1, fval, df_dx, df_dy, ierr)
486 :
487 0 : if(ierr/=0) then
488 0 : write(*,*) 'failed in eval_eosCMS_fixed_X'
489 : return
490 : end if
491 :
492 0 : res = fval
493 0 : dres_dlnT = df_dx * iln10
494 0 : dres_dlnRho = df_dy * iln10
495 :
496 : ! specific heats were log
497 0 : res(i_Cp) = exp(fval(i_Cp))
498 0 : dres_dlnT(i_Cp) = dres_dlnT(i_Cp) * res(i_Cp)
499 0 : dres_dlnRho(i_Cp) = dres_dlnRho(i_Cp) * res(i_Cp)
500 :
501 0 : res(i_Cv) = exp(fval(i_Cv))
502 0 : dres_dlnT(i_Cv) = dres_dlnT(i_Cv) * res(i_Cv)
503 0 : dres_dlnRho(i_Cv) = dres_dlnRho(i_Cv) * res(i_Cv)
504 :
505 0 : end subroutine eval_eosCMS_fixed_X
506 :
507 0 : subroutine locate_logT(c,logT, iT, logT0, logT1)
508 : type(eosCMS_X_info), pointer :: c
509 : real(dp), intent(inout) :: logT
510 : integer, intent(out) :: iT
511 : real(dp), intent(out) :: logT0, logT1
512 0 : iT = int((logT - c% logT_min)/c% delta_logT + 0.0001_dp) + 1
513 0 : if(iT < 1 .or. iT >= c% num_logTs)then
514 0 : if(iT < 1) then
515 0 : iT=1
516 0 : logT0=c% logT_min
517 0 : logT1= logT0 + c% delta_logT
518 0 : logT = logT0
519 : else
520 0 : iT = c% num_logTs -1
521 0 : logT0 = c% logT_min + real(iT-1,kind=dp)*c% delta_logT
522 0 : logT1 = logT0 + c% delta_logT
523 0 : logT = logT1
524 : end if
525 : else
526 0 : logT0 = c% logT_min + (iT-1)*c% delta_logT
527 0 : logT1 = logT0 + c% delta_logT
528 : end if
529 0 : end subroutine locate_logT
530 :
531 :
532 0 : subroutine locate_logRho(c,logRho, iRho, logRho0, logRho1)
533 : type(eosCMS_X_info), pointer :: c
534 : real(dp), intent(inout) :: logRho
535 : integer, intent(out) :: iRho
536 : real(dp), intent(out) :: logRho0, logRho1
537 0 : iRho = int((logRho - c% logRho_min)/c% delta_logRho + 0.0001_dp) + 1
538 0 : if(iRho < 1 .or. iRho >= c% num_logRhos)then
539 0 : if(iRho < 1) then
540 0 : iRho=1
541 0 : logRho0=c% logRho_min
542 0 : logRho1= logRho0 + c% delta_logRho
543 0 : logRho = logRho0
544 : else
545 0 : iRho = c% num_logRhos -1
546 0 : logRho0 = c% logRho_min + real(iRho-1,kind=dp)*c% delta_logRho
547 0 : logRho1 = logRho0 + c% delta_logRho
548 0 : logRho = logRho1
549 : end if
550 : else
551 0 : logRho0 = c% logRho_min + (iRho-1)*c% delta_logRho
552 0 : logRho1 = logRho0 + c% delta_logRho
553 : end if
554 0 : end subroutine locate_logRho
555 :
556 0 : subroutine load_eosCMS_table(iX, ierr)
557 : integer, intent(in) :: iX
558 : integer, intent(out) :: ierr
559 : character(len=256) :: filename, data_sub_dir
560 : character(len=1024) :: message
561 : integer :: io, ios, i, j, n, v, ili_logTs, ili_logRhos, ibcxmin, ibcxmax, ibcymin, ibcymax
562 0 : real(dp), allocatable, target :: f1_ary(:)
563 0 : real(dp), pointer :: f1(:), f(:,:,:), vec(:), tbl(:,:,:,:)
564 0 : real(dp), target :: vec_ary(50)
565 : type(eosCMS_X_Info), pointer :: c
566 0 : real(dp), allocatable :: bcxmin(:), bcxmax(:), bcymin(:), bcymax(:)
567 0 : real(dp) :: X_in, Z_in
568 :
569 0 : ierr=0
570 0 : c => eosCMS_X_data(iX)
571 0 : vec => vec_ary
572 :
573 0 : data_sub_dir = '/eosCMS_data/'
574 0 : filename = trim(mesa_data_dir) // trim(data_sub_dir) // 'mesa-CMS_' // CMS_Xstr(iX) // 'x.data'
575 :
576 0 : open(newunit=io,file=trim(filename), action='read', status='old', iostat=ios)
577 0 : if(ios/=0) then
578 0 : write(*,'(a)') 'failed while reading ' // trim(filename)
579 0 : call mesa_error(__FILE__,__LINE__)
580 0 : close(io)
581 0 : ierr=-1
582 0 : return
583 : end if
584 :
585 0 : read(io,*) !header
586 0 : read(io,'(a)') message
587 0 : call str_to_vector(message, vec, n, ierr)
588 0 : if(ierr/=0) return
589 :
590 0 : c% version = int(vec(1))
591 0 : X_in = vec(2)
592 0 : Z_in = vec(3)
593 0 : c% num_logTs = int(vec(4))
594 0 : c% logT_min= vec(5)
595 0 : c% logT_Max = vec(6)
596 0 : c% delta_logT = vec(7)
597 0 : c% num_logRhos = int(vec(8))
598 0 : c% logRho_min = vec(9)
599 0 : c% logRho_max = vec(10)
600 0 : c% delta_logRho = vec(11)
601 :
602 0 : allocate(c% logTs(c% num_logTs))
603 0 : allocate(c% logRhos(c% num_logRhos))
604 :
605 0 : c% logTs(1) = c% logT_min
606 0 : do i = 2, c% num_logTs
607 0 : c% logTs(i) = c% logTs(i-1) + c% delta_logT
608 : end do
609 :
610 0 : c% logRhos(1) = c% logRho_min
611 0 : do i = 2, c% num_logRhos
612 0 : c% logRhos(i) = c% logRhos(i-1) + c% delta_logRho
613 : end do
614 :
615 : !check that the input X value is compatible with the expected value
616 0 : if(abs(X_in - CMS_Xvals(iX)) > 0.01_dp) then
617 0 : write(*,*) ' eosCMS X value does not match table '
618 0 : write(*,*) ' expected X = ', CMS_Xvals(iX)
619 0 : write(*,*) ' received X = ', X_in
620 0 : call mesa_error(__FILE__,__LINE__)
621 0 : ierr=-1
622 0 : return
623 : end if
624 :
625 0 : read(io,*) !header
626 0 : read(io,*) !header
627 :
628 0 : allocate(c% tbl1(sz_per_eos_point * nv * c% num_logRhos * c% num_logTs))
629 : tbl(1:sz_per_eos_point, 1:nv, 1:c% num_logTs, 1:c% num_logRhos) => &
630 0 : c% tbl1(1:sz_per_eos_point*nv*c% num_logTs*c% num_logRhos)
631 :
632 0 : allocate(f1_ary(sz_per_eos_point * c% num_logRhos * c% num_logTs))
633 0 : f1 => f1_ary
634 : f(1:sz_per_eos_point,1:c% num_logTs,1:c% num_logRhos) => &
635 0 : f1(1:sz_per_eos_point*c% num_logTs*c% num_logRhos)
636 :
637 0 : do i=1,c% num_logTs
638 0 : do j=1,c% num_logRhos
639 0 : read(io,'(a)') message
640 0 : call str_to_vector(message, vec, n, ierr)
641 0 : if(ierr/=0)then
642 0 : close(io)
643 0 : return
644 : end if
645 0 : tbl(1,i_lnPgas,i,j) = vec(3)*ln10
646 0 : tbl(1,i_lnE,i,j) = vec(4)*ln10
647 0 : tbl(1,i_lnS,i,j) = vec(5)*ln10
648 0 : tbl(1,i_chiRho,i,j) = vec(6)
649 0 : tbl(1,i_chiT,i,j) = vec(7)
650 0 : tbl(1,i_Cp,i,j) = safe_log(vec(8))
651 0 : tbl(1,i_Cv,i,j) = safe_log(vec(9))
652 0 : tbl(1,i_dE_dRho,i,j) = vec(10)
653 0 : tbl(1,i_dS_dT,i,j) = vec(11)
654 0 : tbl(1,i_dS_dRho,i,j) = vec(12)
655 0 : tbl(1,i_mu,i,j) = vec(13)
656 0 : tbl(1,i_lnfree_e,i,j) = vec(14)
657 0 : tbl(1,i_gamma1,i,j) = vec(15)
658 0 : tbl(1,i_gamma3,i,j) = vec(16)
659 0 : tbl(1,i_grad_ad,i,j) = vec(17)
660 0 : tbl(1,i_eta,i,j) = vec(18)
661 : end do
662 : end do
663 :
664 0 : close(io)
665 :
666 : ! logT is "x"
667 : ! logRho is "y"
668 :
669 : ! "not a knot" bc's at edges of tables
670 0 : allocate(bcxmin(c% num_logRhos), bcxmax(c% num_logRhos))
671 0 : allocate(bcymin(c% num_logTs), bcymax(c% num_logTs))
672 0 : ibcxmin = 0; bcxmin(:) = 0
673 0 : ibcxmax = 0; bcxmax(:) = 0
674 0 : ibcymin = 0; bcymin(:) = 0
675 0 : ibcymax = 0; bcymax(:) = 0
676 :
677 : !create table for bicubic spline
678 0 : do v = 1, nv
679 0 : do j = 1, c% num_logRhos
680 0 : do i = 1, c% num_logTs
681 0 : f(1,i,j) = tbl(1,v,i,j)
682 : end do
683 : end do
684 :
685 : call interp_mkbicub_db( &
686 : c% logTs, c% num_logTs, c% logRhos, c% num_logRhos, f1, c% num_logTs, &
687 : ibcxmin, bcxmin, ibcxmax, bcxmax, ibcymin, bcymin, ibcymax, bcymax, &
688 0 : ili_logTs, ili_logRhos, ierr)
689 :
690 0 : do j = 1, c% num_logRhos
691 0 : do i = 1, c% num_logTs
692 0 : tbl(2,v,i,j) = f(2,i,j)
693 0 : tbl(3,v,i,j) = f(3,i,j)
694 0 : tbl(4,v,i,j) = f(4,i,j)
695 : end do
696 : end do
697 : end do
698 :
699 0 : if(ierr==0) eosCMS_X_loaded(iX) = .true.
700 :
701 0 : end subroutine load_eosCMS_table
702 :
703 0 : end module eoscms_eval
|
