Line data Source code
1 : ! ***********************************************************************
2 : !
3 : ! Copyright (C) 2010-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 brunt
21 :
22 : use star_private_def
23 : use const_def, only: dp, pi4, crad
24 : use utils_lib
25 :
26 : implicit none
27 :
28 : private
29 : public :: do_brunt_B
30 : public :: do_brunt_N2
31 :
32 : logical, parameter :: dbg = .false.
33 :
34 : contains
35 :
36 : ! call this before mlt
37 56329 : subroutine do_brunt_B(s,nzlo,nzhi,ierr)
38 : use star_utils, only: get_face_values
39 : use utils_lib, only: set_nan
40 : use interp_1d_def
41 :
42 : type (star_info), pointer :: s
43 : integer, intent(in) :: nzlo, nzhi
44 : integer, intent(out) :: ierr
45 : integer :: nz, k
46 45 : real(dp), allocatable, dimension(:) :: smoothing_array
47 :
48 : include 'formats'
49 :
50 45 : ierr = 0
51 45 : nz = s% nz
52 :
53 45 : if (.not. s% calculate_Brunt_B) then
54 0 : call set_nan(s% brunt_B(1:nz))
55 0 : call set_nan(s% unsmoothed_brunt_B(1:nz))
56 0 : return
57 : end if
58 :
59 45 : if (s% use_other_brunt) then
60 0 : call s% other_brunt(s% id, ierr)
61 0 : if (ierr /= 0) then
62 0 : s% retry_message = 'failed in other_brunt'
63 0 : if (s% report_ierr) write(*, *) s% retry_message
64 : end if
65 : else
66 45 : call do_brunt_B_MHM_form(s,nzlo,nzhi,ierr)
67 45 : if (ierr /= 0) then
68 0 : s% retry_message = 'failed in do_brunt_B_MHM_form'
69 0 : if (s% report_ierr) write(*, *) s% retry_message
70 : end if
71 : end if
72 45 : if (ierr /= 0) return
73 :
74 : ! save unsmoothed brunt_B
75 56329 : do k=nzlo,nzhi
76 56284 : s% unsmoothed_brunt_B(k) = s% brunt_B(k)
77 56329 : if (is_bad(s% unsmoothed_brunt_B(k))) then
78 0 : write(*,2) 'unsmoothed_brunt_B(k)', k, s% unsmoothed_brunt_B(k)
79 0 : call mesa_error(__FILE__,__LINE__,'brunt')
80 : end if
81 : end do
82 :
83 45 : allocate(smoothing_array(nz))
84 45 : call smooth_brunt_B(smoothing_array)
85 45 : if (s% use_other_brunt_smoothing) then
86 0 : call s% other_brunt_smoothing(s% id, ierr)
87 0 : if (ierr /= 0) then
88 0 : s% retry_message = 'failed in other_brunt_smoothing'
89 0 : if (s% report_ierr) write(*, *) s% retry_message
90 0 : return
91 : end if
92 : end if
93 :
94 : contains
95 :
96 45 : subroutine smooth_brunt_B(work)
97 45 : use star_utils, only: threshold_smoothing
98 : real(dp) :: work(:)
99 : logical, parameter :: preserve_sign = .false.
100 45 : if (s% num_cells_for_smooth_brunt_B <= 0) return
101 : call threshold_smoothing( &
102 : s% brunt_B, s% threshold_for_smooth_brunt_B, s% nz, &
103 45 : s% num_cells_for_smooth_brunt_B, preserve_sign, work)
104 45 : end subroutine smooth_brunt_B
105 :
106 : end subroutine do_brunt_B
107 :
108 :
109 : ! call this after mlt
110 45 : subroutine do_brunt_N2(s,nzlo,nzhi,ierr)
111 : use star_utils, only: get_face_values
112 :
113 : type (star_info), pointer :: s
114 : integer, intent(in) :: nzlo, nzhi
115 : integer, intent(out) :: ierr
116 :
117 45 : real(dp) :: f
118 : real(dp), allocatable, dimension(:) :: &
119 45 : rho_P_chiT_chiRho, rho_P_chiT_chiRho_face
120 :
121 : integer :: nz, k
122 :
123 : include 'formats'
124 :
125 45 : ierr = 0
126 45 : nz = s% nz
127 :
128 45 : if (.not. (s% calculate_Brunt_B .and. s% calculate_Brunt_N2)) then
129 0 : call set_nan(s% brunt_N2(1:nz))
130 0 : call set_nan(s% brunt_N2_composition_term(1:nz))
131 0 : return
132 : end if
133 :
134 45 : allocate(rho_P_chiT_chiRho(nz), rho_P_chiT_chiRho_face(nz))
135 :
136 56329 : do k=1,nz
137 56284 : rho_P_chiT_chiRho(k) = (s% rho(k)/s% Peos(k))*(s% chiT(k)/s% chiRho(k))
138 : ! correct for difference between gravitational mass density and baryonic mass density (rho)
139 56329 : if (s% use_mass_corrections) then
140 0 : rho_P_chiT_chiRho(k) = s% mass_correction(k)*rho_P_chiT_chiRho(k)
141 : end if
142 : end do
143 :
144 : call get_face_values( &
145 45 : s, rho_P_chiT_chiRho, rho_P_chiT_chiRho_face, ierr)
146 45 : if (ierr /= 0) then
147 0 : s% retry_message = 'failed in get_face_values'
148 0 : if (s% report_ierr) write(*, *) s% retry_message
149 0 : return
150 : end if
151 :
152 56329 : do k=1,nz ! clip B and calculate N^2 from B
153 : if (abs(s% brunt_B(k)) < s% min_magnitude_brunt_B .or. &
154 56284 : s% gradT(k) == 0 .or. is_bad(s% gradT_sub_grada(k))) then
155 0 : s% brunt_B(k) = 0
156 0 : s% brunt_N2(k) = 0
157 0 : s% brunt_N2_composition_term(k) = 0
158 0 : cycle
159 : end if
160 56284 : f = pow2(s% grav(k))*rho_P_chiT_chiRho_face(k)
161 56284 : if (is_bad(f) .or. is_bad(s% brunt_B(k)) .or. is_bad(s% gradT_sub_grada(k))) then
162 0 : write(*,2) 'f', k, f
163 0 : write(*,2) 's% brunt_B(k)', k, s% brunt_B(k)
164 0 : write(*,2) 's% gradT_sub_grada(k)', k, s% gradT_sub_grada(k)
165 0 : write(*,2) 's% gradT(k)', k, s% gradT(k)
166 0 : write(*,2) 's% grada_face(k)', k, s% grada_face(k)
167 0 : call mesa_error(__FILE__,__LINE__,'brunt')
168 : end if
169 56284 : s% brunt_N2(k) = f*(s% brunt_B(k) - s% gradT_sub_grada(k))
170 56329 : s% brunt_N2_composition_term(k) = f*s% brunt_B(k)
171 : end do
172 :
173 45 : if (s% brunt_N2_coefficient /= 1d0) then
174 0 : do k=1,nz
175 0 : s% brunt_N2(k) = s% brunt_N2_coefficient*s% brunt_N2(k)
176 : s% brunt_N2_composition_term(k) = &
177 0 : s% brunt_N2_coefficient*s% brunt_N2_composition_term(k)
178 : end do
179 : end if
180 :
181 90 : end subroutine do_brunt_N2
182 :
183 :
184 45 : subroutine do_brunt_B_MHM_form(s, nzlo, nzhi, ierr)
185 : ! Brassard from Mike Montgomery (MHM)
186 45 : use star_utils, only: get_face_values
187 : use interp_1d_def
188 :
189 : type (star_info), pointer :: s
190 : integer, intent(in) :: nzlo, nzhi
191 : integer, intent(out) :: ierr
192 :
193 :
194 45 : real(dp), allocatable, dimension(:) :: T_face, rho_face, chiT_face, chiRho_face
195 : integer :: nz, species, k, op_err
196 : logical, parameter :: dbg = .false.
197 :
198 : include 'formats'
199 :
200 45 : ierr = 0
201 :
202 45 : nz = s% nz
203 45 : species = s% species
204 :
205 45 : allocate(T_face(nz), rho_face(nz), chiT_face(nz), chiRho_face(nz))
206 :
207 45 : call get_face_values(s, s% chiT, chiT_face, ierr)
208 45 : if (ierr /= 0) return
209 :
210 45 : call get_face_values(s, s% chiRho, chiRho_face, ierr)
211 45 : if (ierr /= 0) return
212 :
213 45 : call get_face_values(s, s% T, T_face, ierr)
214 45 : if (ierr /= 0) return
215 :
216 45 : call get_face_values(s, s% rho, rho_face, ierr)
217 45 : if (ierr /= 0) return
218 :
219 45 : !$OMP PARALLEL DO PRIVATE(k,op_err) SCHEDULE(dynamic,2)
220 : do k=nzlo,nzhi
221 : op_err = 0
222 : call get_brunt_B(&
223 : s, species, nz, k, T_face(k), rho_face(k), chiT_face(k), chiRho_face(k), op_err)
224 : if (op_err /= 0) ierr = op_err
225 : end do
226 : !$OMP END PARALLEL DO
227 :
228 225 : end subroutine do_brunt_B_MHM_form
229 :
230 :
231 56284 : subroutine get_brunt_B(s, species, nz, k, T_face, rho_face, chiT_face, chiRho_face, ierr)
232 45 : use eos_def, only: num_eos_basic_results, num_eos_d_dxa_results, i_lnPgas
233 : use eos_support, only: get_eos
234 :
235 : type (star_info), pointer :: s
236 : integer, intent(in) :: species, nz, k
237 : real(dp), intent(in) :: T_face, rho_face, chiT_face, chiRho_face
238 : integer, intent(out) :: ierr
239 :
240 56284 : real(dp) :: lnP1, lnP2, logRho_face, logT_face, Prad_face, &
241 56284 : alfa, Ppoint, dlnP_dm, delta_lnP, delta_lnMbar
242 : real(dp), dimension(num_eos_basic_results) :: &
243 4502720 : res, d_eos_dlnd, d_eos_dlnT
244 1463384 : real(dp) :: d_eos_dxa(num_eos_d_dxa_results,species)
245 :
246 : logical, parameter :: dbg = .false.
247 :
248 : include 'formats'
249 :
250 56284 : ierr = 0
251 56284 : s% brunt_B(k) = 0
252 56284 : if (k <= 1) return
253 :
254 56239 : logT_face = log10(T_face)
255 56239 : logRho_face = log10(rho_face)
256 56239 : Prad_face = crad * T_face*T_face*T_face*T_face / 3
257 :
258 56239 : if (is_bad_num(logT_face) .or. is_bad_num(logRho_face)) then
259 0 : ierr = -1
260 0 : return
261 : write(*,2) 'logT_face', k, logT_face
262 : write(*,2) 'logRho_face', k, logRho_face
263 : write(*,'(A)')
264 : write(*,2) 's% dq(k-1)', k-1, s% dq(k-1)
265 : write(*,2) 's% dq(k)', k, s% dq(k)
266 : write(*,'(A)')
267 : write(*,2) 's% T(k-1)', k-1, s% T(k-1)
268 : write(*,2) 'T_face', k, T_face
269 : write(*,2) 's% T(k)', k, s% T(k)
270 : write(*,'(A)')
271 : write(*,2) 's% rho(k-1)', k-1, s% rho(k-1)
272 : write(*,2) 'rho_face', k, rho_face
273 : write(*,2) 's% rho(k)', k, s% rho(k)
274 : write(*,'(A)')
275 : write(*,2) 's% chiT(k-1)', k-1, s% chiT(k-1)
276 : write(*,2) 'chiT_face', k, chiT_face
277 : write(*,2) 's% chiT(k)', k, s% chiT(k)
278 : write(*,'(A)')
279 : call mesa_error(__FILE__,__LINE__,'get_brunt_B')
280 : end if
281 :
282 : call get_eos( &
283 : s, 0, s% xa(:,k), &
284 : rho_face, logRho_face, T_face, logT_face, &
285 : res, d_eos_dlnd, d_eos_dlnT, &
286 56239 : d_eos_dxa, ierr)
287 56239 : if (ierr /= 0) return
288 56239 : lnP1 = log(Prad_face + exp(res(i_lnPgas)))
289 :
290 : call get_eos( &
291 : s, 0, s% xa(:,k-1), &
292 : rho_face, logRho_face, T_face, logT_face, &
293 : res, d_eos_dlnd, d_eos_dlnT, &
294 56239 : d_eos_dxa, ierr)
295 56239 : if (ierr /= 0) return
296 56239 : lnP2 = log(Prad_face + exp(res(i_lnPgas)))
297 :
298 56239 : delta_lnP = s% lnPeos(k-1) - s% lnPeos(k)
299 56239 : if (delta_lnP > -1d-50) then
300 2 : alfa = s% dq(k-1)/(s% dq(k-1) + s% dq(k))
301 2 : Ppoint = alfa*s% Peos(k) + (1-alfa)*s% Peos(k-1)
302 2 : dlnP_dm = -s% cgrav(k)*s% m(k)/(pi4*pow4(s% r(k))*Ppoint)
303 2 : delta_lnP = dlnP_dm*s% dm_bar(k)
304 : end if
305 :
306 56239 : s% brunt_B(k) = (lnP1 - lnP2)/delta_lnP/chiT_face
307 :
308 : ! add term accounting for the composition-related gradient in gravitational mass
309 56239 : if (s% use_mass_corrections) then
310 0 : delta_lnMbar = log(s% mass_correction(k-1)) - log(s% mass_correction(k))
311 0 : s% brunt_B(k) = s% brunt_B(k) - chiRho_face*delta_lnMbar/delta_lnP/chiT_face
312 : end if
313 :
314 56239 : if (is_bad_num(s% brunt_B(k))) then
315 0 : ierr = -1
316 0 : s% retry_message = 'bad num for brunt_B'
317 0 : if (s% report_ierr) then
318 0 : write(*,2) 's% brunt_B(k)', k, s% brunt_B(k)
319 0 : write(*,2) 'chiT_face', k, chiT_face
320 0 : write(*,2) 'delta_lnP', k, delta_lnP
321 0 : write(*,2) 's% lnPeos(k)', k, s% lnPeos(k)
322 0 : write(*,2) 's% lnPeos(k-1)', k-1, s% lnPeos(k-1)
323 0 : write(*,2) 'lnP1', k, lnP1
324 0 : write(*,2) 'lnP2', k, lnP2
325 0 : write(*,'(A)')
326 : end if
327 0 : if (s% stop_for_bad_nums) then
328 0 : write(*,2) 's% brunt_B(k)', k, s% brunt_B(k)
329 0 : call mesa_error(__FILE__,__LINE__,'do_brunt_B_MHM_form')
330 : end if
331 : end if
332 :
333 56284 : end subroutine get_brunt_B
334 :
335 : end module brunt
|
