Line data Source code
1 : ! ***********************************************************************
2 : !
3 : ! Copyright (C) 2014-2019 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 adjust_net
21 :
22 : use star_private_def
23 :
24 : implicit none
25 :
26 : private
27 : public :: check_adjust_net
28 :
29 : contains
30 :
31 0 : subroutine check_adjust_net(s, species, &
32 : min_x_for_keep, min_x_for_n, min_x_for_add, &
33 : max_Z_for_add, max_N_for_add, max_A_for_add, ierr)
34 : use adjust_xyz, only: change_net
35 : use utils_lib, only: mkdir
36 : use chem_def, only: &
37 : chem_isos, ineut, max_el_z, el_name, element_min_N, element_max_N
38 : use chem_lib, only: chem_get_iso_id
39 : use star_utils, only: get_string_for_model_number
40 : use net_lib, only: show_net_reactions
41 :
42 : type (star_info), pointer :: s
43 : integer, intent(in) :: species
44 : real(dp), intent(in) :: &
45 : min_x_for_keep, min_x_for_n, min_x_for_add, max_Z_for_add, max_N_for_add, max_A_for_add
46 : integer, intent(out) :: ierr
47 :
48 : character (len=strlen) :: &
49 : net_name, fname, temp_fname, z_plus_n_str, cname, line_buf, species_str
50 : integer :: num_digits, i, j, io, io_new, nz, cid, Z, N, A, &
51 : new_species, next_model_number
52 : integer, dimension(0:max_el_z) :: min_N, max_N
53 0 : real(dp) :: max_x_for_species(species), max_x
54 0 : logical :: have_new_isos, still_in_net(species)
55 : logical, parameter :: adjust_abundances_for_new_isos = .false.
56 :
57 : include 'formats'
58 :
59 0 : ierr = 0
60 0 : nz = s% nz
61 0 : min_N = -1
62 0 : max_N = -1
63 :
64 0 : next_model_number = s% model_number + 1
65 :
66 0 : call include_iso(0,1) ! neut
67 0 : call include_iso(1,0) ! h1
68 0 : call include_iso(2,2) ! he4
69 0 : call include_iso(6,6) ! c12
70 0 : call include_iso(6,7) ! c13
71 0 : call include_iso(7,7) ! n14
72 0 : call include_iso(8,8) ! o16
73 0 : call include_iso(10,10) ! ne20
74 0 : call include_iso(10,12) ! ne22
75 0 : call include_iso(12,12) ! mg24
76 0 : call include_iso(12,14) ! mg26
77 0 : call include_iso(14,14) ! si28
78 0 : call include_iso(16,16) ! s32
79 0 : call include_iso(26,30) ! fe56
80 :
81 0 : do j=1,species
82 0 : max_x = maxval(s% xa(j,1:nz))
83 0 : max_x_for_species(j) = max_x
84 0 : if (max_x < min_x_for_keep) cycle
85 0 : cid = s% chem_id(j)
86 0 : Z = chem_isos% Z(cid)
87 0 : N = chem_isos% N(cid)
88 0 : A=Z+N
89 0 : call include_iso(Z,N)
90 0 : if (max_x >= min_x_for_n .and. N+1 <= max_N_for_add .and. A+1 <= max_A_for_add) then
91 0 : call include_iso(Z,N+1) ! (n,g)
92 0 : call include_iso(Z,N-1) ! (g,n)
93 : end if
94 0 : if (max_x >= min_x_for_add) then
95 0 : if (Z+1 <= max_Z_for_add.and. A+1 <= max_A_for_add) then
96 0 : call include_iso(Z+1,N) ! (p,g)
97 0 : call include_iso(Z-1,N) ! (g,p)
98 : end if
99 0 : if (Z+2 <= max_Z_for_add .and. N+2 <= max_N_for_add .and. A+4 <= max_A_for_add) then
100 0 : call include_iso(Z+2,N+2) ! (a,g)
101 0 : call include_iso(Z-2,N-2) ! (g,a)
102 : end if
103 0 : if (Z+2 <= max_Z_for_add .and. N+1 <= max_N_for_add .and. A+3 <= max_A_for_add) then
104 0 : call include_iso(Z+2,N+1) ! (a,n)
105 0 : call include_iso(Z-2,N-1) ! (n,a)
106 : end if
107 0 : if (Z+1 <= max_Z_for_add .and. N+2 <= max_N_for_add .and. A+3 <= max_A_for_add) then
108 0 : call include_iso(Z+1,N+2) ! (a,p)
109 0 : call include_iso(Z-1,N-2) ! (p,a)
110 : end if
111 0 : if (Z+1 <= max_Z_for_add .and. N+1 <= max_N_for_add .and. A+2 <= max_A_for_add) then
112 0 : call include_iso(Z+1,N-1) ! (p,n)
113 0 : call include_iso(Z-1,N+1) ! (n,p)
114 : end if
115 0 : if (Z+4 <= max_Z_for_add .and. N+4 <= max_N_for_add .and. A+8 <= max_A_for_add) then
116 0 : call include_iso(Z+4,N+4) ! (2a,g) ! extend alpha chain by 2
117 0 : call include_iso(Z+3,N+4) ! (2a,p) ! extend alpha chain by 2
118 : end if
119 : end if
120 : end do
121 :
122 0 : temp_fname = '.temp_net'
123 0 : open(newunit=io, file=trim(temp_fname), action='write', iostat=ierr)
124 0 : if (ierr /= 0) then
125 0 : write(*,*) 'failed to open ' // trim(temp_fname)
126 0 : return
127 : end if
128 :
129 0 : write(io,'(a)') 'add_isos_and_reactions('
130 0 : write(io,'(a)') '! iso log10 max x'
131 :
132 0 : new_species = 0
133 0 : have_new_isos = .false.
134 0 : still_in_net(:) = .false.
135 0 : do Z=0,max_el_z
136 0 : if (min_N(Z) < 0 .or. max_N(Z) < min_N(Z)) cycle
137 0 : if (Z == 0) then
138 0 : i = s% net_iso(ineut)
139 0 : if (i == 0) then
140 0 : write(*,*) ' add ' // trim(el_name(Z))
141 0 : have_new_isos = .true.
142 0 : write(io,'(3x,a8,a)') trim(el_name(Z)), ' ! newly added'
143 : else
144 0 : still_in_net(i) = .true.
145 0 : write(io,'(3x,a8,a,f10.5)') trim(el_name(Z)), ' ! ', &
146 0 : log10(max(1d-199,max_x_for_species(i)))
147 : end if
148 0 : new_species = new_species + 1
149 0 : cycle
150 : end if
151 0 : write(io,'(A)')
152 0 : do N = min_N(Z), max_N(Z)
153 0 : write(z_plus_n_str,'(i4)') Z+N
154 0 : write(cname,'(a)') trim(el_name(Z)) // trim(adjustl(z_plus_n_str))
155 0 : cid = chem_get_iso_id(cname)
156 0 : if (cid <= 0) cycle
157 0 : i = s% net_iso(cid)
158 0 : if (i == 0) then
159 0 : write(*,*) ' add ' // trim(cname)
160 0 : have_new_isos = .true.
161 0 : write(io,'(3x,a8,a)') trim(cname), ' ! newly added'
162 : else
163 0 : still_in_net(i) = .true.
164 0 : write(io,'(3x,a8,a,f10.5)') trim(cname), ' ! ', &
165 0 : log10(max(1d-199,max_x_for_species(i)))
166 : end if
167 0 : new_species = new_species + 1
168 : end do
169 : end do
170 :
171 0 : do i=1,species
172 0 : if (still_in_net(i)) cycle
173 0 : write(*,*) ' drop ' // trim(chem_isos% name(s% chem_id(i)))
174 : end do
175 :
176 0 : write(io,'(a)') ')'
177 0 : write(io,'(A)')
178 :
179 0 : close(io)
180 :
181 0 : if (new_species == species .and. .not. have_new_isos) return
182 :
183 0 : open(newunit=io, file=trim(temp_fname), action='read', status='old', iostat=ierr)
184 0 : if (ierr /= 0) then
185 0 : write(*,*) 'failed to open ' // trim(temp_fname)
186 0 : return
187 : end if
188 :
189 0 : num_digits = 5
190 : call get_string_for_model_number( &
191 0 : '', next_model_number, num_digits, net_name)
192 0 : write(species_str,'(i4)') new_species
193 0 : net_name = trim(net_name) // '_' // trim(adjustl(species_str)) // '.net'
194 0 : fname = 'nets/' // trim(net_name)
195 0 : call mkdir('nets/')
196 0 : open(newunit=io_new, file=trim(fname), action='write', iostat=ierr)
197 0 : if (ierr /= 0) then
198 0 : write(*,*) 'failed to open ' // trim(fname)
199 0 : return
200 : end if
201 :
202 0 : write(io_new,'(a,i8,a)') '! ' // trim(net_name), new_species, ' species'
203 0 : write(io_new,*)
204 :
205 0 : do i=1,10000
206 0 : read(io, fmt='(a)', iostat=ierr) line_buf
207 0 : if (ierr /= 0) exit
208 0 : write(io_new, fmt='(a)') trim(line_buf)
209 : end do
210 0 : ierr = 0
211 :
212 0 : close(io)
213 0 : close(io_new)
214 :
215 0 : write(*,'(i11,a,i8,a)') next_model_number, &
216 0 : ' change to ' // trim(fname), new_species, ' species'
217 :
218 : call change_net( &
219 0 : s% id, adjust_abundances_for_new_isos, net_name, ierr)
220 0 : if (ierr /= 0) then
221 0 : write(*,*) 'failed in change_net ' // trim(net_name)
222 0 : call mesa_error(__FILE__,__LINE__,'check_adjust_net')
223 0 : return
224 : end if
225 :
226 0 : if (net_name /= s% net_name) then
227 0 : write(*,*) ' new net_name ', trim(net_name)
228 0 : write(*,*) 'old s% net_name ', trim(s% net_name)
229 0 : write(*,*) 'failed to change'
230 0 : call mesa_error(__FILE__,__LINE__,'check_adjust_net')
231 : end if
232 :
233 0 : s% using_revised_net_name = .true.
234 0 : s% revised_net_name = s% net_name
235 0 : s% need_to_setvars = .true.
236 :
237 : contains
238 :
239 0 : subroutine include_iso(Z,N)
240 : integer, intent(in) :: Z,N
241 0 : if (Z < 0 .or. Z > max_el_z .or. N < 0) return
242 0 : if (N < element_min_N(Z) .or. N > element_max_N(Z)) return
243 0 : if (N < min_N(Z) .or. min_N(Z) < 0) min_N(Z) = N
244 0 : if (N > max_N(Z)) max_N(Z) = N
245 0 : end subroutine include_iso
246 :
247 : end subroutine check_adjust_net
248 :
249 : end module adjust_net
|
