Line data Source code
1 : ! ***********************************************************************
2 : !
3 : ! Copyright (C) 2012-2019 Phil Arras & 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 create_initial_model
21 :
22 : use star_private_def
23 : use const_def, only: dp, pi, pi4, mp, lsun, msun, standard_cgrav, boltzm, boltz_sigma, &
24 : arg_not_provided, two_thirds, four_thirds_pi
25 : use chem_def
26 :
27 : implicit none
28 :
29 : private
30 : public :: build_initial_model
31 :
32 : integer :: eos_handle, kap_handle, species
33 :
34 : real(dp) :: X, Z, Y, G, abar, zbar, z2bar, z53bar, ye
35 : real (dp) :: Tmin,eps,R_try
36 :
37 : integer, parameter :: max_species=1000
38 : real(dp) :: xa(max_species)
39 : integer, pointer, dimension(:) :: net_iso, chem_id
40 :
41 : integer, parameter :: nzmax=100000
42 :
43 : type create_star_info
44 : integer :: nz
45 : real(dp), dimension(nzmax) :: rg, mg, Pg, Tg, taug, Lg, rhog, intdmTg
46 : real(dp) :: mass, radius, Teff, luminosity
47 : end type create_star_info
48 :
49 : integer, parameter :: max_create_star_handles = 3
50 : type (create_star_info), target, save :: &
51 : create_star_handles(max_create_star_handles)
52 :
53 : contains
54 :
55 0 : subroutine get_create_star_ptr(id,cs,ierr)
56 : integer, intent(in) :: id
57 : type (create_star_info), pointer :: cs
58 : integer, intent(out) :: ierr
59 : include 'formats'
60 0 : if (id < 1 .or. id > max_create_star_handles) then
61 0 : ierr = -1
62 0 : write(*,2) 'bad id for get_create_star_ptr', id
63 0 : return
64 : end if
65 0 : cs => create_star_handles(id)
66 0 : ierr = 0
67 : end subroutine get_create_star_ptr
68 :
69 :
70 0 : subroutine build_initial_model(s, ierr)
71 : use chem_lib, only: basic_composition_info, chem_Xsol
72 : use adjust_xyz, only: get_xa_for_standard_metals
73 : use alloc, only: allocate_star_info_arrays
74 : use star_utils, only: store_r_in_xh, store_rho_in_xh, store_T_in_xh
75 :
76 : type (star_info), pointer :: s
77 : integer, intent(out) :: ierr
78 :
79 : integer :: initial_zfracs, i_lum, i, j, k, itry, max_try, id0, id1, id2
80 0 : real(dp) :: M, R, initial_y, initial_h1, initial_h2, initial_he3, initial_he4, &
81 0 : S0, Pc0, e0(2), S1, Pc1, e1(2), S2, Pc2, e2(2), det, dPc, dS, safefac, &
82 0 : initial_z, xsol_he3, xsol_he4, mass_correction, mat(2,2), minv(2,2), sumx
83 : type (create_star_info), pointer :: cs
84 :
85 : include 'formats'
86 :
87 : ierr = 0
88 :
89 0 : if (s% use_other_build_initial_model) then
90 0 : call s% other_build_initial_model(s% id, ierr)
91 0 : return
92 : end if
93 :
94 0 : id0 = 1; id1 = 2; id2 = 3
95 0 : call get_create_star_ptr(id0, cs, ierr)
96 0 : if (ierr /= 0) call mesa_error(__FILE__,__LINE__,'build_initial_model')
97 :
98 0 : if (s% nvar_hydro > 4) then
99 0 : write(*,*) 'sorry, current build_initial_model only supports the basic 4 vars.'
100 0 : ierr = -1
101 0 : return
102 : end if
103 :
104 0 : M = s% mass_in_gm_for_create_initial_model
105 0 : R = s% radius_in_cm_for_create_initial_model
106 :
107 0 : if (M <= 0) then
108 0 : write(*,1) 'must set mass_in_gm_for_create_initial_model', M
109 0 : ierr = -1
110 0 : return
111 : end if
112 :
113 0 : if (R <= 0) then
114 0 : write(*,1) 'must set radius_in_cm_for_create_initial_model', R
115 0 : ierr = -1
116 0 : return
117 : end if
118 :
119 0 : initial_z = s% initial_z
120 0 : initial_y = s% initial_y
121 0 : initial_h1 = max(0d0, min(1d0, 1d0 - (initial_z + initial_y)))
122 0 : initial_h2 = 0d0
123 0 : xsol_he3 = chem_Xsol('he3')
124 0 : xsol_he4 = chem_Xsol('he4')
125 0 : initial_he3 = initial_y*xsol_he3/(xsol_he3 + xsol_he4)
126 0 : initial_he4 = initial_y*xsol_he4/(xsol_he3 + xsol_he4)
127 0 : initial_zfracs = s% initial_zfracs_for_create_initial_model
128 :
129 0 : chem_id => s% chem_id
130 0 : net_iso => s% net_iso
131 0 : eos_handle = s% eos_handle
132 0 : kap_handle = s% kap_handle
133 0 : species = s% species
134 :
135 : call get_xa_for_standard_metals(s, &
136 : species, s% chem_id, s% net_iso, &
137 : initial_h1, initial_h2, initial_he3, initial_he4, &
138 : initial_zfracs, &
139 0 : s% initial_dump_missing_heaviest, xa, ierr)
140 0 : if (ierr /= 0) then
141 0 : write(*,*) 'create initial model failed in get_xa_for_standard_metals'
142 0 : return
143 : end if
144 :
145 : call basic_composition_info( &
146 : species, s% chem_id, xa(:), x, y, z, abar, zbar, z2bar, z53bar, ye, &
147 0 : mass_correction, sumx)
148 :
149 0 : G = standard_cgrav
150 0 : Tmin = 0.d0 ! sets surface isotherm
151 0 : eps = s% initial_model_eps ! integration accuracy
152 0 : R_try = R ! used to set grid spacing near the center
153 :
154 : ! init guess
155 0 : Pc0=exp10(s% center_logP_1st_try_for_create_initial_model)
156 0 : S0=boltzm/mp * s% entropy_1st_try_for_create_initial_model
157 :
158 0 : max_try = s% max_tries_for_create_initial_model
159 0 : do itry = 1, max_try+1
160 :
161 0 : if (itry > max_try) then
162 : write(*,'(a,i10)') &
163 0 : 'create initial model failed to converge in allowed number of tries', max_try
164 0 : ierr = -1
165 0 : return
166 : end if
167 :
168 0 : write(*,2) 'create initial model iteration', itry
169 :
170 0 : Pc1=Pc0*1.01d0
171 0 : S1=S0
172 :
173 0 : Pc2=Pc0
174 0 : S2=S0*1.01d0
175 :
176 0 : !$OMP PARALLEL DO PRIVATE(i) SCHEDULE(dynamic,2)
177 : do i=0,2
178 : select case(i)
179 : case(0)
180 : call PSerrfunc(id0,Pc0,S0,M,R,e0)
181 : case(1)
182 : call PSerrfunc(id1,Pc1,S1,M,R,e1)
183 : case(2)
184 : call PSerrfunc(id2,Pc2,S2,M,R,e2)
185 : end select
186 : end do
187 : !$OMP END PARALLEL DO
188 :
189 0 : if (abs(e0(1)) < s% abs_e01_tolerance_for_create_initial_model .and. &
190 : abs(e0(2)) < s% abs_e02_tolerance_for_create_initial_model) exit
191 :
192 0 : mat(1,1)=(e1(1)-e0(1))/(0.01d0*Pc0)
193 0 : mat(1,2)=(e2(1)-e0(1))/(0.01d0*S0)
194 0 : mat(2,1)=(e1(2)-e0(2))/(0.01d0*Pc0)
195 0 : mat(2,2)=(e2(2)-e0(2))/(0.01d0*S0)
196 0 : det = mat(1,1)*mat(2,2)-mat(1,2)*mat(2,1)
197 0 : minv(1,1)=mat(2,2)/det
198 0 : minv(2,2)=mat(1,1)/det
199 0 : minv(1,2)=-mat(1,2)/det
200 0 : minv(2,1)=-mat(2,1)/det
201 0 : dPc = - (minv(1,1)*e0(1)+minv(1,2)*e0(2))
202 0 : dS = - (minv(2,1)*e0(1)+minv(2,2)*e0(2))
203 :
204 0 : safefac = 1.d0 + 10.d0*max( abs(dPc/Pc0), abs(dS/S0) )
205 0 : Pc0 = Pc0 + dPc/safefac
206 0 : S0 = S0 + dS/safefac
207 :
208 : end do
209 :
210 0 : s% nz = cs% nz - 1 ! skip center point
211 0 : call allocate_star_info_arrays(s, ierr)
212 0 : if (ierr /= 0) then
213 : return
214 : end if
215 :
216 0 : s% M_center = 0d0
217 0 : s% L_center = 0d0
218 0 : s% R_center = 0d0
219 0 : s% v_center = 0d0
220 :
221 0 : s% mstar = cs% mass
222 0 : s% star_mass = cs% mass/Msun
223 0 : s% xmstar = cs% mass
224 :
225 0 : i_lum = s% i_lum
226 :
227 0 : do k=1, s% nz
228 0 : i = s% nz - k + 2 ! skip center point
229 0 : call store_rho_in_xh(s, k, cs% rhog(i))
230 0 : call store_T_in_xh(s, k, cs% Tg(i))
231 0 : call store_r_in_xh(s, k, cs% rg(i))
232 0 : if (i_lum /= 0) s% xh(i_lum, k) = cs% Lg(i)
233 0 : do j=1,species
234 0 : s% xa(j,k) = xa(j)
235 : end do
236 0 : s% q(k) = cs% mg(i)/s% xmstar
237 : end do
238 0 : s% dq(s% nz) = s% q(s% nz)
239 0 : do k=1, s% nz - 1
240 0 : s% dq(k) = s% q(k) - s% q(k+1)
241 : end do
242 :
243 0 : write(*,*) 'done build_initial_model'
244 :
245 0 : end subroutine build_initial_model
246 :
247 :
248 : ! output
249 : !errvec(1)=(mass-M)/M
250 : !errvec(2)=(radius-R)/R
251 0 : subroutine PSerrfunc(id,Pcguess,Sguess,M,R,errvec)
252 : integer, intent(in) :: id
253 : real(dp), intent(in) :: Pcguess,Sguess,M,R
254 : real(dp), intent(out) :: errvec(2)
255 :
256 : integer :: ierr
257 : real(dp) :: rhoc,Tc,Pc,S
258 0 : real(dp) :: P,T,rho,y(3),dydP(3),r1,m1,intdmT1,dy1(3),dy2(3),dy3(3),dy4(3),d_P
259 0 : real(dp) :: kap,tau,grav
260 : logical :: exitnow
261 : integer :: i, nz, k
262 : type (create_star_info), pointer :: cs
263 :
264 0 : call get_create_star_ptr(id, cs, ierr)
265 0 : if (ierr /= 0) call mesa_error(__FILE__,__LINE__,'PSerrfunc')
266 :
267 0 : Pc=Pcguess
268 0 : S=Sguess
269 :
270 :
271 : ! PA. given S and Pc, integrate adiabat outward
272 :
273 : !S = 13.0 * boltzm / mp
274 : !Pc = 1.0d13
275 :
276 0 : call get_TRho_from_PS(cs,Pc,S,Tc,rhoc)
277 :
278 : ! initialize grids to zero
279 0 : cs% rg=0.d0
280 0 : cs% mg=0.d0
281 0 : cs% Pg=0.d0
282 0 : cs% Tg=0.d0
283 0 : cs% rhog=0.d0
284 0 : cs% intdmTg=0.d0
285 0 : cs% taug=0.d0
286 0 : cs% Lg=0.d0
287 :
288 : ! record central point
289 0 : i=1
290 0 : cs% rg(i)=0.d0
291 0 : cs% mg(i)=0.d0
292 0 : cs% Pg(i) = Pc
293 0 : cs% Tg(i) = Tc
294 0 : cs% rhog(i) = rhoc
295 0 : cs% intdmTg(i)=0.d0
296 0 : cs% taug(i) = 1.d20
297 :
298 : ! compute first point off center
299 0 : r1 = 1.d-2 * R_try * min(0.1d0,eps)
300 :
301 0 : m1 = four_thirds_pi * rhoc * r1*r1*r1
302 0 : P = Pc - two_thirds*pi * G*rhoc*rhoc*r1*r1
303 0 : intdmT1=m1*Tc
304 0 : y=[r1,m1,intdmT1]
305 0 : call get_TRho_from_PS(cs,P,S,T,rho)
306 :
307 : ! record first point off center
308 0 : i=2
309 0 : cs% rg(i)=r1
310 0 : cs% mg(i)=m1
311 0 : cs% Pg(i) = P
312 0 : cs% Tg(i) = T
313 0 : cs% rhog(i) = rho
314 0 : cs% intdmTg(i) = intdmT1
315 0 : cs% taug(i) = 1.d20
316 0 : cs% Lg(i)=0.d0
317 :
318 0 : exitnow=.false.
319 0 : do
320 :
321 : ! what should stepsize be
322 0 : call derivs(cs,P,S,y,dydP)
323 0 : d_P=-eps*min( P , abs(y(1)/dydP(1)) , abs(y(2)/dydP(2)) )
324 :
325 : ! 4th order rk step
326 0 : dy1=dydP*d_P
327 0 : call derivs(cs,P+d_P/2d0,S,y+dy1/2d0,dydP)
328 0 : dy2=dydP*d_P
329 0 : call derivs(cs,P+d_P/2d0,S,y+dy2/2d0,dydP)
330 0 : dy3=dydP*d_P
331 0 : call derivs(cs,P+d_P,S,y+dy3,dydP)
332 0 : dy4=dydP*d_P
333 0 : y=y+dy1/6.d0+dy2/3.d0+dy3/3.d0+dy4/6.d0
334 0 : P=P+d_P
335 :
336 0 : call get_TRho_from_PS(cs,P,S,T,rho)
337 :
338 0 : call get_kap_from_rhoT(cs,log10(rho),log10(T),kap)
339 0 : grav = G*y(2)/(y(1)*y(1))
340 0 : tau = kap * P / grav
341 :
342 0 : i=i+1
343 0 : cs% rg(i)=y(1)
344 0 : cs% mg(i)=y(2)
345 0 : cs% Pg(i) = P
346 0 : cs% Tg(i) = T
347 0 : cs% rhog(i) = rho
348 0 : cs% taug(i) = tau
349 0 : cs% intdmTg(i) = y(3)
350 :
351 0 : if (T<150.d0) then
352 0 : print *,"temp too low in integration"
353 0 : stop
354 : end if
355 :
356 0 : if (tau < 2.d0/3.d0) then
357 0 : exitnow=.true.
358 : end if
359 :
360 0 : if (exitnow) exit
361 :
362 : end do
363 :
364 : ! interpolate last point to tau=2/3
365 0 : nz=i
366 0 : cs% nz = nz
367 0 : P = cs% Pg(i-1) + (cs% Pg(i)-cs% Pg(i-1))*(2.d0/3.d0-cs% taug(i-1))/(cs% taug(i)-cs% taug(i-1))
368 0 : cs% rg(i) = cs% rg(i-1) + (cs% rg(i)-cs% rg(i-1))*(P-cs% Pg(i-1))/(cs% Pg(i)-cs% Pg(i-1))
369 0 : cs% mg(i) = cs% mg(i-1) + (cs% mg(i)-cs% mg(i-1))*(P-cs% Pg(i-1))/(cs% Pg(i)-cs% Pg(i-1))
370 0 : cs% Tg(i) = cs% Tg(i-1) + (cs% Tg(i)-cs% Tg(i-1))*(P-cs% Pg(i-1))/(cs% Pg(i)-cs% Pg(i-1))
371 0 : cs% rhog(i) = cs% rhog(i-1) + (cs% rhog(i)-cs% rhog(i-1))*(P-cs% Pg(i-1))/(cs% Pg(i)-cs% Pg(i-1))
372 0 : cs% intdmTg(i) = cs% intdmTg(i-1) + (cs% intdmTg(i)-cs% intdmTg(i-1))*(P-cs% Pg(i-1))/(cs% Pg(i)-cs% Pg(i-1))
373 0 : cs% Pg(i)=P
374 0 : cs% taug(i)=2.d0/3.d0
375 :
376 0 : cs% mass = cs% mg(nz)
377 0 : cs% radius = cs% rg(nz)
378 0 : cs% Teff = cs% Tg(nz)
379 0 : cs% luminosity = pi4*pow2(cs% radius)*boltz_sigma*pow4(cs% Teff)
380 0 : do k=1,nz
381 0 : cs% Lg(k)=cs% luminosity*cs% intdmTg(k)/cs% intdmTg(nz)
382 : end do
383 0 : write(*,"(a20,2x,es15.8)") "log10(L/Lsun)=", log10(cs%luminosity/Lsun)
384 0 : write(*,"(a20,2x,es15.8)") "Mass=", cs% mass
385 0 : write(*,"(a20,2x,es15.8)") "Radius=", cs% radius
386 0 : write(*,*) ''
387 :
388 0 : errvec(1)=(cs% mass-M)/M
389 0 : errvec(2)=(cs% radius-R)/R
390 :
391 0 : end subroutine PSerrfunc
392 :
393 :
394 : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
395 :
396 :
397 : ! y(1)=r
398 : ! y(2)=m
399 : ! y(3)=\int_0^m dm T(m)
400 :
401 0 : subroutine derivs(cs,P,S,y,dydP)
402 : type (create_star_info), pointer :: cs
403 : real(dp), intent(in) :: P,S,y(3)
404 : real(dp), intent(out) :: dydP(3)
405 0 : real(dp) :: r,m,T,rho,intdmT
406 :
407 0 : r=y(1)
408 0 : m=y(2)
409 0 : intdmT=y(3)
410 0 : call get_TRho_from_PS(cs,P,S,T,rho)
411 : !write(*,"(a24,10(2x,es15.8))") "S*mp/kb,lgPc,lgTc,rhoc=",S*mp/boltzm,log10(P),log10(T),rho
412 0 : dydP(1)=-r*r/(G*m*rho)
413 0 : dydP(2)=-pi4*r*r*r*r/(G*m)
414 0 : dydP(3)=dydP(2)*T
415 :
416 0 : end subroutine derivs
417 :
418 :
419 : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
420 :
421 :
422 0 : subroutine get_kap_from_rhoT(cs,logrho,logT,kap)
423 : use kap_def, only: num_kap_fracs
424 : use kap_lib
425 : type (create_star_info), pointer :: cs
426 : real(dp), intent(in) :: logrho,logT
427 : real(dp), intent(out) :: kap
428 : real(dp) :: lnfree_e, d_lnfree_e_dlnRho, d_lnfree_e_dlnT
429 : real(dp) :: eta, d_eta_dlnRho, d_eta_dlnT
430 : real(dp) :: dlnkap_dlnRho, dlnkap_dlnT
431 0 : real(dp) :: kap_fracs(num_kap_fracs), dlnkap_dxa(species)
432 : integer :: ierr
433 :
434 : ierr=0
435 :
436 : ! this ignores lnfree and eta
437 0 : lnfree_e=0; d_lnfree_e_dlnRho=0; d_lnfree_e_dlnT=0
438 0 : eta=0; d_eta_dlnRho=0; d_eta_dlnT=0
439 :
440 : call kap_get( &
441 : kap_handle, species, chem_id, net_iso, xa, &
442 : logRho, logT, &
443 : lnfree_e, d_lnfree_e_dlnRho, d_lnfree_e_dlnT, &
444 : eta, d_eta_dlnRho, d_eta_dlnT, &
445 0 : kap_fracs, kap, dlnkap_dlnRho, dlnkap_dlnT, dlnkap_dxa, ierr)
446 0 : if (ierr /= 0) then
447 0 : write(*,*) 'failed in kap_get get_kap_from_rhoT'
448 : return
449 : end if
450 :
451 : if(ierr/=0) then
452 : write(*,*) 'kap_get failed'
453 : stop
454 : end if
455 :
456 0 : end subroutine get_kap_from_rhoT
457 :
458 :
459 : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
460 :
461 : ! ignore radiation pressure: P=Pgas here
462 :
463 0 : subroutine get_TRho_from_PS(cs,P,S,T,rho)
464 0 : use eos_lib
465 : use eos_def
466 : type (create_star_info), pointer :: cs
467 : real(dp), intent(in) :: P,S
468 : real(dp), intent(out) :: T,rho
469 :
470 0 : real(dp) :: logT_result,log10Rho,dlnRho_dlnPgas_const_T,dlnRho_dlnT_const_Pgas
471 : real(dp), dimension(num_eos_basic_results) :: &
472 0 : res, d_dlnRho_const_T, d_dlnT_const_Rho
473 0 : real(dp), dimension(num_eos_d_dxa_results, species) :: d_dxa_const_TRho
474 : integer, parameter :: max_iter = 100
475 : integer :: eos_calls,ierr
476 : real(dp), parameter :: logT_tol = 1.d-6, other_tol = 1.d-6, logT_guess = 4.d0, &
477 : logT_bnd1= arg_not_provided, logT_bnd2= arg_not_provided, &
478 : other_at_bnd1= arg_not_provided, other_at_bnd2= arg_not_provided
479 :
480 : call eosPT_get_T( &
481 : eos_handle, &
482 : species, chem_id, net_iso, xa, &
483 : log10(P), i_lnS, log(S), &
484 : logT_tol, other_tol, max_iter, logT_guess, &
485 : logT_bnd1, logT_bnd2, other_at_bnd1, other_at_bnd2, &
486 : logT_result, rho, log10Rho, dlnRho_dlnPgas_const_T, dlnRho_dlnT_const_Pgas, &
487 : res, d_dlnRho_const_T, d_dlnT_const_Rho, &
488 : d_dxa_const_TRho, &
489 0 : eos_calls, ierr)
490 0 : if (ierr /=0) then
491 0 : print *,"failure in eosPT_get_T"
492 0 : stop
493 : end if
494 0 : T = exp10(logT_result)
495 :
496 : ! don't let T get below min temp. use to make a surface isotherm.
497 0 : if (T < Tmin) T=Tmin
498 :
499 0 : end subroutine get_TRho_from_PS
500 :
501 0 : end module create_initial_model
|
