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 adjust_mesh_support
21 :
22 : use star_private_def
23 : use const_def, only: dp
24 :
25 : implicit none
26 :
27 : private
28 : public :: get_gval_info, check_validity
29 :
30 : logical, parameter :: dbg = .false.
31 :
32 : contains
33 :
34 10 : subroutine get_gval_info( &
35 : s, delta_gval_max, gvals1, nz, &
36 : num_gvals, gval_names, &
37 : gval_is_xa_function, gval_is_logT_function, ierr)
38 : use chem_def
39 : use num_lib, only: find0
40 : use rates_def
41 : use alloc
42 : use mesh_functions, only: set_mesh_function_data, max_allowed_gvals
43 : type (star_info), pointer :: s
44 : real(dp), dimension(:), pointer :: delta_gval_max
45 : real(dp), dimension(:), pointer :: gvals1
46 : integer, intent(in) :: nz, num_gvals
47 : integer, intent(out) :: ierr
48 : character (len=32), intent(out) :: gval_names(max_allowed_gvals)
49 : logical, intent(out), dimension(max_allowed_gvals) :: &
50 : gval_is_xa_function, gval_is_logT_function
51 :
52 : integer :: j, k
53 : logical, parameter :: dbg = .false.
54 10 : real(dp), allocatable, dimension(:) :: src
55 : real(dp) :: eps_min_for_delta, &
56 : dlog_eps_dlogP_full_off, dlog_eps_dlogP_full_on
57 : real(dp), dimension(:,:), pointer :: gvals
58 :
59 10 : gvals(1:nz,1:num_gvals) => gvals1(1:nz*num_gvals)
60 :
61 : include 'formats'
62 :
63 10 : ierr = 0
64 :
65 10 : eps_min_for_delta = exp10(s% mesh_dlog_eps_min_for_extra)
66 :
67 10 : dlog_eps_dlogP_full_off = s% mesh_dlog_eps_dlogP_full_off
68 10 : dlog_eps_dlogP_full_on = s% mesh_dlog_eps_dlogP_full_on
69 :
70 : call set_mesh_function_data( &
71 : s, num_gvals, gval_names, &
72 10 : gval_is_xa_function, gval_is_logT_function, gvals1, ierr)
73 10 : if (ierr /= 0) return
74 :
75 10 : allocate(src(nz))
76 :
77 10 : call set_delta_gval_max(src, ierr)
78 10 : if (ierr /= 0) return
79 :
80 30 : call smooth_gvals(nz,src,num_gvals,gvals)
81 :
82 :
83 : contains
84 :
85 :
86 10 : subroutine set_delta_gval_max(src, ierr)
87 : real(dp) :: src(:)
88 : integer, intent(out) :: ierr
89 10 : real(dp) :: P_exp, beta
90 :
91 : logical, parameter :: dbg = .false.
92 :
93 : include 'formats'
94 :
95 10 : ierr = 0
96 12105 : delta_gval_max(1:nz) = 1d0
97 :
98 10 : if (s% mesh_Pgas_div_P_exponent /= 0) then
99 0 : P_exp = s% mesh_Pgas_div_P_exponent
100 0 : do k=1,nz
101 0 : beta = s% Pgas(k)/s% Peos(k)
102 0 : delta_gval_max(k) = delta_gval_max(k)*pow(beta,P_exp)
103 : end do
104 : end if
105 :
106 10 : if (s% use_other_mesh_delta_coeff_factor) then
107 0 : do k=1,nz
108 0 : s% mesh_delta_coeff_factor(k) = delta_gval_max(k)
109 : end do
110 0 : call s% other_mesh_delta_coeff_factor(s% id, ierr)
111 0 : if (ierr /= 0) then
112 0 : write(*,*) 'other_mesh_delta_coeff_factor returned ierr', ierr
113 0 : return
114 : end if
115 0 : do k=1,nz
116 0 : delta_gval_max(k) = s% mesh_delta_coeff_factor(k)
117 : end do
118 : end if
119 :
120 100 : do j=1,num_mesh_logX
121 90 : if (s% mesh_dlogX_dlogP_extra(j) > 0 .and. s% mesh_dlogX_dlogP_extra(j) < 1) &
122 10 : call do_mesh_dlogX_dlogP_coef(s,j)
123 : end do
124 :
125 10 : call do1_dlog_eps_dlogP_coef(s% mesh_dlog_pp_dlogP_extra, ipp)
126 10 : call do1_dlog_eps_dlogP_coef(s% mesh_dlog_cno_dlogP_extra, icno)
127 10 : call do1_dlog_eps_dlogP_coef(s% mesh_dlog_3alf_dlogP_extra, i3alf)
128 :
129 10 : call do1_dlog_eps_dlogP_coef(s% mesh_dlog_burn_c_dlogP_extra, i_burn_c)
130 10 : call do1_dlog_eps_dlogP_coef(s% mesh_dlog_burn_n_dlogP_extra, i_burn_n)
131 10 : call do1_dlog_eps_dlogP_coef(s% mesh_dlog_burn_o_dlogP_extra, i_burn_o)
132 10 : call do1_dlog_eps_dlogP_coef(s% mesh_dlog_burn_ne_dlogP_extra, i_burn_ne)
133 10 : call do1_dlog_eps_dlogP_coef(s% mesh_dlog_burn_na_dlogP_extra, i_burn_na)
134 10 : call do1_dlog_eps_dlogP_coef(s% mesh_dlog_burn_mg_dlogP_extra, i_burn_mg)
135 10 : call do1_dlog_eps_dlogP_coef(s% mesh_dlog_burn_si_dlogP_extra, i_burn_si)
136 10 : call do1_dlog_eps_dlogP_coef(s% mesh_dlog_burn_s_dlogP_extra, i_burn_s)
137 10 : call do1_dlog_eps_dlogP_coef(s% mesh_dlog_burn_ar_dlogP_extra, i_burn_ar)
138 10 : call do1_dlog_eps_dlogP_coef(s% mesh_dlog_burn_ca_dlogP_extra, i_burn_ca)
139 10 : call do1_dlog_eps_dlogP_coef(s% mesh_dlog_burn_ti_dlogP_extra, i_burn_ti)
140 10 : call do1_dlog_eps_dlogP_coef(s% mesh_dlog_burn_cr_dlogP_extra, i_burn_cr)
141 10 : call do1_dlog_eps_dlogP_coef(s% mesh_dlog_burn_fe_dlogP_extra, i_burn_fe)
142 :
143 10 : call do1_dlog_eps_dlogP_coef(s% mesh_dlog_cc_dlogP_extra, icc)
144 10 : call do1_dlog_eps_dlogP_coef(s% mesh_dlog_co_dlogP_extra, ico)
145 10 : call do1_dlog_eps_dlogP_coef(s% mesh_dlog_oo_dlogP_extra, ioo)
146 :
147 10 : call do1_dlog_eps_dlogP_coef(s% mesh_dlog_pnhe4_dlogP_extra, ipnhe4)
148 10 : call do1_dlog_eps_dlogP_coef(s% mesh_dlog_photo_dlogP_extra, iphoto)
149 10 : call do1_dlog_eps_dlogP_coef(s% mesh_dlog_other_dlogP_extra, iother)
150 :
151 10 : if (s% mesh_delta_coeff_factor_smooth_iters > 0) then ! smooth delta_gval_max
152 :
153 12105 : do k=1,nz
154 12105 : src(k) = delta_gval_max(k)
155 : end do
156 :
157 30 : do j=1,s% mesh_delta_coeff_factor_smooth_iters
158 30 : delta_gval_max(1) = (2*src(1) + src(2))/3
159 36255 : do k=2,nz-1
160 36255 : delta_gval_max(k) = (src(k-1) + src(k) + src(k+1))/3
161 : end do
162 30 : delta_gval_max(nz) = (2*src(nz) + src(nz-1))/3
163 30 : if (j == 3) exit
164 20 : src(1) = (2*delta_gval_max(1) + delta_gval_max(2))/3
165 24170 : do k=2,nz-1
166 : src(k) = &
167 24170 : (delta_gval_max(k-1) + delta_gval_max(k) + delta_gval_max(k+1))/3
168 : end do
169 30 : src(nz) = (2*delta_gval_max(nz) + delta_gval_max(nz-1))/3
170 : end do
171 :
172 : end if
173 :
174 10 : end subroutine set_delta_gval_max
175 :
176 :
177 0 : subroutine do_mesh_dlogX_dlogP_coef(s, which)
178 : use chem_lib, only: chem_get_iso_id
179 : type (star_info), pointer :: s
180 : integer, intent(in) :: which
181 : real(dp) :: &
182 0 : logX_min_for_extra, dlogX_dlogP_extra, dlogX_dlogP_full_on, dlogX_dlogP_full_off, &
183 0 : X_min_for_extra, dlogX, dlogP, dlogX_dlogP, coef
184 : integer :: k, cid, j
185 : include 'formats'
186 0 : if (len_trim(s% mesh_logX_species(which)) == 0) return
187 0 : cid = chem_get_iso_id(s% mesh_logX_species(which))
188 0 : if (cid <= 0) return
189 0 : j = s% net_iso(cid)
190 0 : if (j == 0) return
191 0 : logX_min_for_extra = s% mesh_logX_min_for_extra(which)
192 0 : dlogX_dlogP_extra = s% mesh_dlogX_dlogP_extra(which)
193 0 : dlogX_dlogP_full_on = s% mesh_dlogX_dlogP_full_on(which)
194 0 : dlogX_dlogP_full_off = s% mesh_dlogX_dlogP_full_off(which)
195 0 : X_min_for_extra = exp10(max(-50d0, logX_min_for_extra))
196 0 : do k=2, nz
197 0 : if (s% xa(j,k) < X_min_for_extra .or. s% xa(j,k-1) < X_min_for_extra) cycle
198 0 : dlogX = abs(log10(s% xa(j,k)/s% xa(j,k-1)))
199 0 : dlogP = max(1d-10, abs(s% lnPeos(k) - s% lnPeos(k-1))*ln10)
200 0 : dlogX_dlogP = dlogX/dlogP
201 0 : if (dlogX_dlogP <= dlogX_dlogP_full_off) cycle
202 0 : if (dlogX_dlogP >= dlogX_dlogP_full_on) then
203 : coef = dlogX_dlogP_extra
204 : else
205 : coef = 1 - (1 - dlogX_dlogP_extra) * &
206 0 : (dlogX_dlogP - dlogX_dlogP_full_off) / (dlogX_dlogP_full_on - dlogX_dlogP_full_off)
207 : end if
208 0 : if (coef < delta_gval_max(k)) delta_gval_max(k) = coef
209 0 : if (coef < delta_gval_max(k-1)) delta_gval_max(k-1) = coef
210 0 : if (k < nz) then
211 0 : if (coef < delta_gval_max(k+1)) delta_gval_max(k+1) = coef
212 : end if
213 : end do
214 0 : end subroutine do_mesh_dlogX_dlogP_coef
215 :
216 :
217 220 : subroutine do1_dlog_eps_dlogP_coef(dlog_eps_dlogP_extra, cat)
218 : real(dp), intent(in) :: dlog_eps_dlogP_extra
219 : integer, intent(in) :: cat
220 : integer :: k
221 220 : real(dp) :: eps, epsm1, dlog_eps, dlogP, dlog_eps_dlogP, &
222 220 : extra, new_max, maxv
223 : include 'formats'
224 :
225 410 : if (dlog_eps_dlogP_extra <= 0 .or. dlog_eps_dlogP_extra >=1) return
226 254415 : maxv = maxval(s% eps_nuc_categories(cat,1:nz))
227 210 : if (maxv < eps_min_for_delta) return
228 :
229 24190 : do k=2, nz
230 :
231 24170 : eps = s% eps_nuc_categories(cat,k)
232 24170 : if (eps < eps_min_for_delta) cycle
233 :
234 8340 : epsm1 = s% eps_nuc_categories(cat,k-1)
235 8340 : if (epsm1 < eps_min_for_delta) cycle
236 :
237 216320 : maxv = maxval(s% eps_nuc_categories(:,k))
238 8320 : if (maxv /= eps) cycle
239 :
240 4688 : dlogP = (s% lnPeos(k) - s% lnPeos(k-1))/ln10
241 4688 : if (abs(dlogP) < 1d-50) cycle
242 :
243 4688 : dlog_eps = abs(log10(eps/epsm1))
244 :
245 4688 : if (is_bad(dlog_eps)) then
246 0 : write(*,3) 'adjust mesh support ' // trim(category_name(cat)), &
247 0 : cat, k, s% eps_nuc_categories(cat,k)
248 0 : stop
249 : end if
250 :
251 4688 : dlog_eps_dlogP = dlog_eps/dlogP
252 4688 : if (dlog_eps_dlogP <= dlog_eps_dlogP_full_off) cycle
253 :
254 4592 : if (dlog_eps_dlogP >= dlog_eps_dlogP_full_on) then
255 : extra = dlog_eps_dlogP_extra
256 : else
257 : extra = 1 - (1 - dlog_eps_dlogP_extra) * &
258 : (dlog_eps_dlogP - dlog_eps_dlogP_full_off) / &
259 2168 : (dlog_eps_dlogP_full_on - dlog_eps_dlogP_full_off)
260 : end if
261 4592 : new_max = extra
262 4592 : if (new_max < delta_gval_max(k)) then
263 1572 : delta_gval_max(k) = new_max
264 : end if
265 4592 : if (new_max < delta_gval_max(k-1)) then
266 1547 : delta_gval_max(k-1) = new_max
267 : end if
268 4612 : if (k < nz) then
269 4582 : if (new_max < delta_gval_max(k+1)) delta_gval_max(k+1) = new_max
270 : end if
271 :
272 : end do
273 0 : end subroutine do1_dlog_eps_dlogP_coef
274 :
275 :
276 : function blend_coef(coef_start, alfa) result(coef)
277 : real(dp), intent(in) :: coef_start, alfa
278 : real(dp) :: coef
279 : real(dp) :: beta
280 :
281 : ! this implements the following piecewise blend
282 : ! 0 < alfa < 0.5 : constant (coeff_start)
283 : ! 0.5 < alfa < 1 : linear (coeff_start -> 1)
284 :
285 : beta = 2d0 * (alfa - 0.5d0)
286 : if (beta < 0d0) then
287 : coef = coef_start
288 : else
289 : coef = coef_start*(1d0 - beta) + beta
290 : end if
291 :
292 : end function blend_coef
293 :
294 :
295 : end subroutine get_gval_info
296 :
297 :
298 : subroutine single_peak(nz,src)
299 : integer, intent(in) :: nz
300 : real(dp), pointer :: src(:)
301 : integer :: k, kmax
302 : real(dp) :: prev, val
303 : kmax = maxloc(src(1:nz), dim=1)
304 : val = src(kmax)
305 : do k=kmax+1,nz
306 : prev = val
307 : val = src(k)
308 : if (val > prev) val = prev
309 : src(k) = val
310 : end do
311 : val = src(kmax)
312 : do k=kmax-1,1,-1
313 : prev = val
314 : val = src(k)
315 : if (val > prev) val = prev
316 : src(k) = val
317 : end do
318 : end subroutine single_peak
319 :
320 :
321 : subroutine increasing_inward(nz,src)
322 : integer, intent(in) :: nz
323 : real(dp), pointer :: src(:)
324 : integer :: k
325 : real(dp) :: prev, val
326 : val = src(1)
327 : do k=2,nz
328 : prev = val
329 : val = src(k)
330 : if (val < prev) val = prev
331 : src(k) = val
332 : end do
333 : end subroutine increasing_inward
334 :
335 :
336 : subroutine decreasing_outward(nz,src)
337 : integer, intent(in) :: nz
338 : real(dp), pointer :: src(:)
339 : integer :: k
340 : real(dp) :: prev, val
341 : val = src(nz)
342 : do k=nz-1,1,-1
343 : prev = val
344 : val = src(k)
345 : if (val > prev) val = prev
346 : src(k) = val
347 : end do
348 : end subroutine decreasing_outward
349 :
350 :
351 : subroutine increasing_outward(nz,src)
352 : integer, intent(in) :: nz
353 : real(dp), pointer :: src(:)
354 : integer :: k
355 : real(dp) :: prev, val
356 : val = src(nz)
357 : do k=nz-1,1,-1
358 : prev = val
359 : val = src(k)
360 : if (val < prev) val = prev
361 : src(k) = val
362 : end do
363 : end subroutine increasing_outward
364 :
365 :
366 10 : subroutine smooth_gvals(nz,src,num_gvals,gvals)
367 : integer, intent(in) :: nz, num_gvals
368 : real(dp) :: src(:), gvals(:,:)
369 : integer :: k, i
370 :
371 40 : do i=1,num_gvals
372 :
373 36315 : do k=1,nz
374 36315 : src(k) = gvals(k,i)
375 : end do
376 :
377 30 : gvals(1,i) = (2*src(1) + src(2))/3
378 36255 : do k=2,nz-1
379 36255 : gvals(k,i) = (src(k-1) + src(k) + src(k+1))/3
380 : end do
381 30 : gvals(nz,i) = (2*src(nz) + src(nz-1))/3
382 :
383 30 : src(1) = (2*gvals(1,i) + gvals(2,i))/3
384 36255 : do k=2,nz-1
385 36255 : src(k) = (gvals(k-1,i) + gvals(k,i) + gvals(k+1,i))/3
386 : end do
387 30 : src(nz) = (2*gvals(nz,i) + gvals(nz-1,i))/3
388 :
389 30 : gvals(1,i) = (2*src(1) + src(2))/3
390 36255 : do k=2,nz-1
391 36255 : gvals(k,i) = (src(k-1) + src(k) + src(k+1))/3
392 : end do
393 40 : gvals(nz,i) = (2*src(nz) + src(nz-1))/3
394 :
395 : end do
396 :
397 10 : end subroutine smooth_gvals
398 :
399 :
400 : subroutine set_boundary_values(s, src, dest, j)
401 : type (star_info), pointer :: s
402 : real(dp), pointer :: src(:), dest(:, :)
403 : integer, intent(in) :: j
404 : integer :: k, nz
405 : nz = s% nz
406 : dest(1,j) = src(1)
407 : do k=2,nz
408 : dest(k,j) = (src(k-1)+src(k))/2
409 : end do
410 : end subroutine set_boundary_values
411 :
412 :
413 10 : subroutine check_validity(s, ierr)
414 : type (star_info), pointer :: s
415 : integer, intent(out) :: ierr
416 :
417 : integer :: k, nz
418 :
419 : include 'formats'
420 :
421 10 : ierr = 0
422 10 : nz = s% nz
423 :
424 12095 : do k=1, nz-1
425 12085 : if (s% xh(s% i_lnR, k) <= s% xh(s% i_lnR, k+1)) then
426 0 : ierr = -1
427 0 : s% retry_message = 'at start of remesh: negative cell volume for cell'
428 0 : if (s% report_ierr) then
429 0 : write(*, *) 'at start of remesh: negative cell volume for cell', k
430 0 : write(*, *)
431 0 : write(*,2) 's% xh(s% i_lnR, k)', k, s% xh(s% i_lnR, k)
432 0 : write(*,2) 's% xh(s% i_lnR, k+1)', k+1, s% xh(s% i_lnR, k+1)
433 0 : write(*, *)
434 0 : write(*, *) 's% model_number', s% model_number
435 0 : write(*, *) 's% nz', s% nz
436 0 : write(*, *) 's% num_retries', s% num_retries
437 0 : write(*, *)
438 : end if
439 0 : return
440 : end if
441 12095 : if (s% dq(k) <= 0) then
442 0 : ierr = -1
443 0 : s% retry_message = 'at start of remesh: non-positive cell mass for cell'
444 0 : if (s% report_ierr) then
445 0 : write(*, *) 'at start of remesh: non-positive cell mass for cell', k
446 0 : write(*, *) 's% model_number', s% model_number
447 0 : write(*, *) 's% nz', s% nz
448 0 : write(*, *) 's% num_retries', s% num_retries
449 0 : write(*, *)
450 : end if
451 0 : return
452 : end if
453 : end do
454 :
455 : end subroutine check_validity
456 :
457 : end module adjust_mesh_support
|
