Line data Source code
1 : ! ***********************************************************************
2 : !
3 : ! Copyright (C) 2012-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 solver_support
21 :
22 : use star_private_def
23 : use utils_lib, only: is_bad
24 : use const_def, only: dp, msun, secyer, ln10, four_thirds_pi
25 : use num_def
26 :
27 : implicit none
28 :
29 : private
30 : public :: force_another_iteration
31 : public :: eval_equations
32 : public :: set_xscale_info
33 : public :: sizequ
34 : public :: sizeB
35 : public :: inspectb
36 : public :: bdomain
37 :
38 : logical, parameter :: dbg = .false.
39 :
40 : contains
41 :
42 11 : subroutine set_xscale_info(s, nvar, ierr)
43 : type (star_info), pointer :: s
44 : integer, intent(in) :: nvar
45 : integer, intent(out) :: ierr
46 :
47 : integer :: i, k, nz, nvar_hydro
48 : real(dp), parameter :: xscale_min = 1
49 :
50 : include 'formats'
51 :
52 11 : ierr = 0
53 :
54 : if (dbg) write(*, *) 'set_xscale'
55 11 : nvar_hydro = s% nvar_hydro
56 11 : nz = s% nz
57 13098 : do k=1,nz
58 170142 : do i=1,nvar
59 170131 : if (i <= nvar_hydro) then ! structure variable
60 52348 : if (i == s% i_j_rot) then
61 0 : s% x_scale(i,k) = 10d0*sqrt(s% cgrav(k)*s% m(k)*s% r_start(k))
62 : else
63 52348 : s% x_scale(i,k) = max(xscale_min, abs(s% xh_start(i,k)))
64 : end if
65 : else ! abundance variable
66 104696 : s% x_scale(i,k) = max(s% xa_scale, s% xa_start(i-nvar_hydro,k))
67 : end if
68 : end do
69 : end do
70 :
71 : contains
72 :
73 : subroutine dump_xscale
74 : integer :: k, j
75 : include 'formats'
76 : !write(*,1) 's% xa_scale', s% xa_scale
77 : do k=1,s% nz
78 : do j=1,nvar
79 : write(*,2) 'xscale ' // trim(s% nameofvar(j)), k, s% x_scale(j,k)
80 : end do
81 : write(*,'(A)')
82 : end do
83 : call mesa_error(__FILE__,__LINE__,'set_xscale')
84 : end subroutine dump_xscale
85 :
86 : end subroutine set_xscale_info
87 :
88 :
89 44 : subroutine eval_equations(s, nvar, ierr)
90 : use hydro_eqns, only: eval_equ
91 : use mix_info, only: set_dxdt_mix
92 : use star_utils, only: update_time, total_times
93 : type (star_info), pointer :: s
94 : integer, intent(in) :: nvar
95 : integer, intent(out) :: ierr
96 :
97 : integer :: cnt, i, j, k, nz
98 : real(dp) :: dt
99 : include 'formats'
100 :
101 44 : ierr = 0
102 44 : nz = s% nz
103 :
104 : if (dbg) write(*, *) 'eval_equations'
105 :
106 44 : dt = s% dt
107 :
108 44 : if (s% solver_iter == 0) then
109 : if (dbg) write(*, *) 'skip set_solver_vars on call before 1st iter'
110 : else
111 : if (dbg) write(*, *) 'call set_solver_vars'
112 33 : call set_solver_vars(s, nvar, dt, ierr)
113 33 : if (ierr /= 0) then
114 0 : if (s% report_ierr) &
115 0 : write(*,2) 'eval_equations: set_solver_vars returned ierr', ierr
116 0 : return
117 : end if
118 : end if
119 :
120 44 : call set_dxdt_mix(s)
121 :
122 44 : if (ierr == 0) then
123 52392 : do k=1,nz
124 680568 : do j=1,nvar
125 628176 : s% equ(j,k) = 0d0
126 628176 : s% residual_weight(j,k) = 1d0
127 680524 : s% correction_weight(j,k) = 1d0
128 : end do
129 : end do
130 : if (dbg) write(*, *) 'call eval_equ'
131 44 : call eval_equ(s, nvar, ierr)
132 44 : if (ierr /= 0) then
133 0 : if (s% report_ierr) &
134 0 : write(*, *) 'eval_equations: eval_equ returned ierr', ierr
135 0 : return
136 : end if
137 : end if
138 :
139 44 : if (ierr /= 0) return
140 :
141 44 : if (.not. s% stop_for_bad_nums) return
142 :
143 44 : cnt = 0
144 52392 : do i=1,nz
145 680568 : do j=1, nvar
146 680524 : if (is_bad_num(s% equ(j, i))) then
147 0 : cnt = cnt + 1
148 0 : s% retry_message = 'eval_equations: equ has a bad num'
149 0 : if (s% report_ierr) then
150 0 : write(*,4) 'eval_equations: equ has a bad num ' // trim(s% nameofequ(j)), &
151 0 : j, i, nvar, s% equ(j, i)
152 0 : write(*,2) 'cell', i
153 0 : write(*,2) 'nz', s% nz
154 : end if
155 0 : if (s% stop_for_bad_nums) call mesa_error(__FILE__,__LINE__,'eval_equations')
156 : end if
157 : end do
158 : end do
159 :
160 44 : if (cnt > 0) then
161 0 : ierr = -1
162 0 : return
163 : end if
164 :
165 :
166 : contains
167 :
168 :
169 : subroutine dump_eval_equ
170 : integer :: k, j
171 : include 'formats'
172 : do k=1,s% nz
173 : do j=1,nvar
174 : write(*,2) 'dx ' // trim(s% nameofvar(j)), k, s% solver_dx(j, k)
175 : end do
176 : write(*,'(A)')
177 : end do
178 : call mesa_error(__FILE__,__LINE__,'dump_eval_equ')
179 44 : end subroutine dump_eval_equ
180 :
181 :
182 : end subroutine eval_equations
183 :
184 :
185 44 : subroutine sizequ(s, nvar, equ_norm, equ_max, k_max, j_max, ierr) ! equ = residuals
186 : type (star_info), pointer :: s
187 : integer, intent(in) :: nvar
188 : real(dp), intent(out) :: equ_norm, equ_max
189 : integer, intent(out) :: k_max, j_max, ierr
190 :
191 : integer :: j, k, num_terms, n, nz, nvar_hydro, nvar_chem, skip_eqn1, skip_eqn2, skip_eqn3
192 44 : real(dp) :: sumequ, absq
193 :
194 : logical :: dbg
195 :
196 : include 'formats'
197 :
198 44 : ierr = 0
199 :
200 44 : equ_norm = 0d0
201 44 : equ_max = 0d0
202 44 : k_max = 0
203 44 : j_max = 0
204 :
205 44 : dbg = s% solver_check_everything
206 :
207 44 : nvar_hydro = min(nvar, s% nvar_hydro)
208 44 : nvar_chem = s% nvar_chem
209 :
210 44 : nz = s% nz
211 44 : n = nz
212 44 : num_terms = 0
213 44 : sumequ = 0
214 44 : skip_eqn1 = 0
215 44 : skip_eqn2 = 0
216 44 : skip_eqn3 = 0
217 44 : if (s% convergence_ignore_equL_residuals) skip_eqn1 = s% i_equL
218 44 : if (s% convergence_ignore_alpha_RTI_residuals) skip_eqn2 = s% i_dalpha_RTI_dt
219 44 : if (s% do_burn .or. s% do_mix) then
220 44 : num_terms = num_terms + nvar*nz
221 44 : if (skip_eqn1 > 0) num_terms = num_terms - nz
222 44 : if (skip_eqn2 > 0) num_terms = num_terms - nz
223 : if (skip_eqn3 > 0) num_terms = num_terms - nz
224 52392 : do k = 1, nz
225 680568 : do j = 1, nvar
226 628176 : if (j == skip_eqn1 .or. j == skip_eqn2 .or. j == skip_eqn3) cycle
227 628176 : if (is_bad(s% equ(j,k)) .or. is_bad(s% residual_weight(j,k))) then
228 0 : ierr = 1
229 0 : return
230 : end if
231 628176 : absq = abs(s% equ(j,k)*s% residual_weight(j,k))
232 628176 : sumequ = sumequ + absq
233 680524 : if (absq > equ_max) then
234 3567 : equ_max = absq
235 3567 : j_max = j
236 3567 : k_max = k
237 : end if
238 : end do
239 : end do
240 : else
241 0 : if (skip_eqn1 == 0 .and. skip_eqn2 == 0) then
242 0 : num_terms = num_terms + nvar_hydro*nz
243 0 : else if (skip_eqn1 > 0 .and. skip_eqn2 > 0) then
244 0 : num_terms = num_terms + (nvar_hydro-2)*nz
245 : else
246 0 : num_terms = num_terms + (nvar_hydro-1)*nz
247 : end if
248 0 : do k = 1, nz
249 0 : do j = 1, nvar_hydro
250 0 : if (j == skip_eqn1 .or. j == skip_eqn2) cycle
251 0 : absq = abs(s% equ(j,k)*s% residual_weight(j,k))
252 0 : sumequ = sumequ + absq
253 0 : if (is_bad(sumequ)) then
254 0 : if (dbg) then
255 0 : write(*,3) trim(s% nameofequ(j)) // ' sumequ', j, k, sumequ
256 0 : call mesa_error(__FILE__,__LINE__,'sizeq 1')
257 : end if
258 0 : ierr = -1
259 0 : if (s% report_ierr) &
260 0 : write(*,3) 'bad equ(j,k)*s% residual_weight(j,k) ' // trim(s% nameofequ(j)), &
261 0 : j, k, s% equ(j,k)*s% residual_weight(j,k)
262 0 : if (s% stop_for_bad_nums) call mesa_error(__FILE__,__LINE__,'sizeq 2')
263 0 : return
264 : end if
265 0 : if (absq > equ_max) then
266 0 : equ_max = absq
267 0 : j_max = j
268 0 : k_max = k
269 : end if
270 : end do
271 : end do
272 : end if
273 44 : if (s% do_burn .or. s% do_mix) then
274 44 : num_terms = num_terms + nvar_chem*nz
275 52392 : do k = 1, nz
276 471176 : do j = nvar_hydro+1, nvar
277 418784 : absq = abs(s% equ(j,k)*s% residual_weight(j,k))
278 418784 : sumequ = sumequ + absq
279 471132 : if (absq > equ_max) then
280 0 : equ_max = absq
281 0 : j_max = j
282 0 : k_max = k
283 : end if
284 : end do
285 : end do
286 : end if
287 :
288 44 : equ_norm = sumequ/num_terms
289 44 : if (dbg) write(*,4) trim(s% nameofequ(j_max)) // ' sizequ equ_max norm', &
290 0 : k_max, s% solver_iter, s% model_number, equ_max, equ_norm
291 :
292 0 : if (dbg) call dump_equ
293 :
294 : return
295 : call dump_equ
296 : call mesa_error(__FILE__,__LINE__,'sizequ')
297 :
298 : contains
299 :
300 0 : subroutine dump_equ
301 : integer :: k, j
302 : include 'formats'
303 0 : do k=1,s% nz
304 0 : do j=1,nvar
305 0 : write(*,3) 'equ ' // trim(s% nameofequ(j)), &
306 0 : k, s% solver_iter, s% equ(j, k)
307 : end do
308 0 : write(*,'(A)')
309 : !if (k == 6) exit
310 : end do
311 0 : end subroutine dump_equ
312 :
313 : end subroutine sizequ
314 :
315 :
316 33 : subroutine sizeB(s, nvar, B, max_correction, correction_norm, max_zone, max_var, ierr)
317 : type (star_info), pointer :: s
318 : integer, intent(in) :: nvar
319 : real(dp), pointer, dimension(:,:) :: B ! (nvar, nz)
320 : real(dp), intent(out) :: correction_norm ! a measure of the average correction
321 : real(dp), intent(out) :: max_correction ! magnitude of the max correction
322 : integer, intent(out) :: max_zone, max_var, ierr
323 :
324 : integer :: k, i, nz, num_terms, j, n, nvar_hydro, jmax, num_xa_terms, &
325 : skip1, skip2, skip3, skip4, skip5
326 33 : real(dp) :: abs_corr, sum_corr, sum_xa_corr, x_limit, &
327 33 : max_abs_correction, max_abs_corr_for_k, max_abs_xa_corr_for_k
328 : logical :: found_NaN, found_bad_num, report
329 : real(dp), parameter :: frac = 0.1d0
330 : logical, parameter :: dbg = .false.
331 : logical, parameter :: check_for_bad_nums = .true.
332 : logical, parameter :: save_max_abs_corr_for_k = .true.
333 :
334 : include 'formats'
335 :
336 : if (dbg) write(*, *) 'enter sizeB'
337 :
338 33 : ierr = 0
339 33 : nz = s% nz
340 33 : n = nz
341 33 : nvar_hydro = min(nvar, s% nvar_hydro)
342 :
343 33 : if (s% include_L_in_correction_limits) then
344 33 : skip1 = 0
345 39294 : do k=1,nz
346 39294 : s% correction_weight(s% i_lum,k) = 1d0/(frac*s% L_start(1) + abs(s% L(k)))
347 : end do
348 : else
349 0 : skip1 = s% i_lum
350 : end if
351 :
352 33 : if (s% u_flag .and. s% include_u_in_correction_limits) then
353 0 : skip2 = 0
354 0 : do k=1,nz
355 0 : s% correction_weight(s% i_u,k) = 1d0/(frac*s% csound_start(k) + abs(s% u(k)))
356 : end do
357 33 : else if (s% v_flag .and. s% include_v_in_correction_limits) then
358 0 : skip2 = 0
359 0 : do k=1,nz
360 0 : s% correction_weight(s% i_v,k) = 1d0/(frac*s% csound_start(k) + abs(s% v(k)))
361 : end do
362 33 : else if (s% u_flag) then
363 0 : skip2 = s% i_u
364 : else
365 33 : skip2 = s% i_v
366 : end if
367 :
368 33 : if (s% RSP2_flag .and. s% include_w_in_correction_limits) then
369 0 : skip3 = 0
370 0 : do k=1,nz
371 0 : if (abs(s% w(k)) < 1d0) then
372 0 : s% correction_weight(s% i_w,k) = 0d0
373 : else
374 0 : s% correction_weight(s% i_w,k) = 1d0/(1d3 + abs(s% w(k))) ! don't sweat the small stuff
375 : end if
376 : end do
377 : else
378 33 : skip3 = s% i_w
379 : end if
380 33 : skip4 = s% i_Hp
381 :
382 33 : skip5 = 0
383 :
384 33 : max_zone = 0
385 33 : max_var = 0
386 33 : num_terms = 0
387 33 : num_xa_terms = 0
388 33 : sum_corr = 0
389 33 : sum_xa_corr = 0
390 33 : max_correction = 0
391 33 : max_abs_correction = 0
392 33 : x_limit = s% correction_xa_limit
393 33 : found_NaN = .false.
394 33 : found_bad_num = .false.
395 33 : report = s% report_ierr
396 39294 : cell_loop: do k = 1, nz
397 39261 : max_abs_corr_for_k = 0
398 39261 : max_abs_xa_corr_for_k = 0
399 :
400 39261 : if (s% do_burn .or. s% do_mix) then
401 : jmax = nvar
402 : else
403 39261 : jmax = nvar_hydro
404 : end if
405 510393 : var_loop: do j = 1, jmax
406 : if (j == skip1 .or. &
407 : j == skip2 .or. &
408 : j == skip3 .or. &
409 : j == skip4 .or. &
410 471132 : j == skip5 .or. &
411 : j == s% i_alpha_RTI) cycle var_loop
412 : if (check_for_bad_nums) then
413 471132 : if (is_bad_num(B(j,k)*s% correction_weight(j,k))) then
414 0 : found_bad_num = .true.
415 0 : if (report) write(*,2) 'sizeB: bad num for correction ' // &
416 0 : s% nameofvar(j), k, B(j,k)*s% correction_weight(j,k)
417 0 : if (s% stop_for_bad_nums) then
418 0 : found_NaN = .true.
419 0 : write(*,3) s% nameofvar(j) // ' B(j,k)*s% correction_weight(j,k)', &
420 0 : j, k, B(j,k)*s% correction_weight(j,k)
421 0 : call mesa_error(__FILE__,__LINE__,'sizeB')
422 : end if
423 :
424 0 : max_zone = k
425 0 : max_var = j
426 0 : exit cell_loop
427 :
428 : cycle var_loop
429 : end if
430 : end if
431 471132 : if (j > nvar_hydro) then
432 314088 : if (s% xa_start(j-nvar_hydro,k) < x_limit) cycle var_loop
433 : end if
434 :
435 280146 : abs_corr = abs(B(j,k)*s% correction_weight(j,k))
436 280146 : if (is_bad_num(abs_corr)) then
437 0 : found_bad_num = .true.
438 0 : if (report) write(*,3) 'B(j,k)*s% correction_weight(j,k)', &
439 0 : j, k, B(j,k)*s% correction_weight(j,k)
440 0 : if (s% stop_for_bad_nums) found_NaN = .true.
441 : end if
442 : if (abs_corr > max_abs_corr_for_k &
443 : .and. .not. (j > nvar_hydro .and. s% ignore_species_in_max_correction)) &
444 : max_abs_corr_for_k = abs_corr
445 : if (abs_corr > max_abs_correction &
446 280146 : .and. .not. (j > nvar_hydro .and. s% ignore_species_in_max_correction)) then
447 11963 : max_correction = B(j,k)*s% correction_weight(j,k)
448 11963 : max_abs_correction = abs_corr
449 11963 : max_zone = k
450 11963 : max_var = j
451 : end if
452 319407 : if (j > nvar_hydro) then
453 123102 : num_xa_terms = num_xa_terms + 1
454 123102 : sum_xa_corr = sum_xa_corr + abs_corr
455 123102 : if (abs_corr > max_abs_xa_corr_for_k) &
456 52001 : max_abs_xa_corr_for_k = abs_corr
457 : else
458 157044 : num_terms = num_terms + 1
459 157044 : sum_corr = sum_corr + abs_corr
460 : end if
461 : end do var_loop
462 39261 : if (num_xa_terms > 0) then
463 39261 : num_terms = num_terms + 1
464 39261 : sum_corr = sum_corr + sum_xa_corr/num_xa_terms
465 : end if
466 39261 : if (s% do_burn .or. s% do_mix) then
467 353349 : species_loop: do j = nvar_hydro+1, nvar
468 314088 : i = j - s% nvar_hydro
469 : if (check_for_bad_nums) then
470 314088 : if (is_bad_num(B(j,k)*s% correction_weight(j,k))) then
471 0 : found_bad_num = .true.
472 0 : if (report) write(*,3) 'chem B(j,k)*s% correction_weight(j,k)', j, k, &
473 0 : B(j,k)*s% correction_weight(j,k)
474 0 : if (s% stop_for_bad_nums) then
475 0 : found_NaN = .true.
476 0 : write(*,3) 'chem B(j,k)*s% correction_weight(j,k)', &
477 0 : j, k, B(j,k)*s% correction_weight(j,k)
478 0 : call mesa_error(__FILE__,__LINE__,'sizeB')
479 0 : max_zone = k
480 0 : max_var = j
481 0 : exit cell_loop
482 : end if
483 : end if
484 : end if
485 : ! recall that correction dx = B*xscale, so B is a relative correction
486 353349 : if (s% xa_start(i,k) >= x_limit) then
487 123102 : abs_corr = abs(B(j,k)*s% correction_weight(j,k))
488 123102 : if (.not. s% ignore_species_in_max_correction) then
489 : if (abs_corr > max_abs_corr_for_k) max_abs_corr_for_k = abs_corr
490 123102 : if (abs_corr > max_abs_correction) then
491 0 : max_abs_correction = abs_corr
492 0 : max_correction = B(j,k)*s% correction_weight(j,k)
493 0 : max_zone = k
494 0 : max_var = j
495 : end if
496 : end if
497 123102 : sum_corr = sum_corr + abs_corr
498 123102 : num_terms = num_terms + 1
499 : end if
500 : end do species_loop
501 : end if
502 39294 : s% max_abs_xa_corr(k) = max_abs_xa_corr_for_k
503 : end do cell_loop
504 :
505 33 : if (found_bad_num) then
506 0 : ierr = -1
507 0 : if (found_NaN .and. s% stop_for_bad_nums) then
508 0 : write(*,*) 'found bad num'
509 0 : call mesa_error(__FILE__,__LINE__,'sizeB')
510 : end if
511 22 : if (.not. dbg) return
512 : end if
513 :
514 33 : if (is_bad_num(sum_corr)) then
515 0 : ierr = -1
516 0 : if (s% stop_for_bad_nums) then
517 0 : if (report) write(*,*) 'sum_corr', sum_corr
518 0 : call mesa_error(__FILE__,__LINE__,'sizeB')
519 : end if
520 0 : if (.not. dbg) return
521 : write(*,*) 'sum_corr', sum_corr
522 : call mesa_error(__FILE__,__LINE__,'sizeB')
523 : end if
524 :
525 33 : correction_norm = sum_corr/num_terms !sqrt(sum_corr/num_terms)
526 : if (dbg) then
527 : write(*,2) 'sizeB: iter, correction_norm, max_correction', &
528 : s% solver_iter, correction_norm, max_correction
529 : if (max_correction > 1d50 .or. is_bad_num(correction_norm)) then
530 : call show_stuff
531 : call mesa_error(__FILE__,__LINE__,'sizeB')
532 : end if
533 : end if
534 :
535 33 : if (s% solver_show_correction_info) call show_stuff
536 :
537 33 : abs_corr = max_abs_correction
538 :
539 33 : s% abs_max_corr2 = s% abs_max_corr1; s% abs_max_corr1 = abs_corr
540 33 : s% max_var2 = s% max_var1; s% max_var1 = max_var
541 33 : s% max_zone2 = s% max_zone1; s% max_zone1 = max_zone
542 :
543 33 : if (ierr /= 0) call mesa_error(__FILE__,__LINE__,'ierr in sizeB')
544 :
545 33 : if (is_bad_num(max_correction)) then
546 0 : ierr = -1
547 0 : if (s% stop_for_bad_nums) then
548 0 : if (report) write(*,*) 'max_correction', max_correction
549 0 : call mesa_error(__FILE__,__LINE__,'sizeB')
550 : end if
551 0 : if (.not. dbg) return
552 : write(*,*) 'max_correction', max_correction
553 : call mesa_error(__FILE__,__LINE__,'sizeB')
554 : end if
555 :
556 33 : if (s% solver_iter < 3) return
557 : ! check for flailing
558 : if ( &
559 : abs_corr > s% tol_max_correction .and. &
560 : abs_corr > s% abs_max_corr1 .and. s% abs_max_corr1 > s% abs_max_corr2 .and. &
561 : max_zone == s% max_zone1 .and. s% max_zone1 == s% max_zone2 .and. &
562 11 : max_var == s% max_var1 .and. s% max_var1 == s% max_var2) then
563 0 : if (s% solver_show_correction_info) then
564 0 : write(*,*) 'give up because diverging'
565 : end if
566 0 : max_correction = 1d99
567 : end if
568 :
569 :
570 : contains
571 :
572 :
573 0 : subroutine show_stuff
574 : integer :: j, k
575 0 : real(dp) :: dx, prev, new
576 : include 'formats'
577 0 : if (s% solver_iter == 1) then
578 0 : write(*,'(A)')
579 : write(*,'(4a7,12a16,99a13)') &
580 0 : 'model', 'iter', 'var', 'zone', &
581 0 : 'corr norm', 'max corr', 'xscale', &
582 0 : 'dx', 'new-prev', 'new', 'prev', &
583 0 : 'mass loc', 'log dt/yr', 'lgE', 'lgT', 'lgRho'
584 : end if
585 0 : k = max_zone
586 0 : j = max_var
587 0 : if (j > nvar_hydro) then
588 0 : prev = s% xa_start(j - nvar_hydro,k)
589 : else
590 0 : prev = s% xh_start(j,k)
591 : end if
592 0 : dx = B(j,k)*s% correction_weight(j,k)*s% x_scale(j,k)
593 0 : new = prev + dx
594 : write(*,'(2i7,a7,i7,12e16.8,99f13.8)') &
595 0 : s% model_number, s% solver_iter, trim(s% nameofvar(max_var)), k, &
596 0 : correction_norm, B(j,k)*s% correction_weight(j,k), s% x_scale(j,k), &
597 0 : dx, new - prev, new, prev, &
598 0 : s% m(k)/Msun, log10(s% dt/secyer), &
599 0 : s% lnE(k)/ln10, s% lnT(k)/ln10, &
600 0 : s% lnd(k)/ln10
601 0 : end subroutine show_stuff
602 :
603 :
604 : subroutine dump_B
605 : integer :: k, j
606 : include 'formats'
607 : do k=1,s% nz
608 : do j=1,nvar
609 : write(*,2) 'B ' // trim(s% nameofequ(j)), k, B(j, k)
610 : end do
611 : write(*,'(A)')
612 : end do
613 : call mesa_error(__FILE__,__LINE__,'dump_equ')
614 : end subroutine dump_B
615 :
616 :
617 : end subroutine sizeB
618 :
619 :
620 : ! the proposed change to dx is B*xscale*correction_factor
621 : ! edit correction_factor and/or B as necessary so that the new dx will be valid.
622 : ! set ierr nonzero if things are beyond repair.
623 33 : subroutine Bdomain(s, nvar, B, correction_factor, ierr)
624 : use const_def, only: dp
625 : use chem_def, only: chem_isos
626 : use star_utils, only: current_min_xa_hard_limit, rand
627 : type (star_info), pointer :: s
628 : integer, intent(in) :: nvar
629 : real(dp), pointer, dimension(:,:) :: B ! (nvar, nz)
630 : real(dp), intent(inout) :: correction_factor
631 : integer, intent(out) :: ierr
632 : integer :: i, j, k, nz, species, bad_j, bad_k
633 33 : real(dp) :: alpha, min_alpha, new_xa, old_xa, dxa, eps, min_xa_hard_limit, &
634 33 : dlum_surf, old_lum_surf, new_lum_surf
635 : include 'formats'
636 33 : ierr = 0
637 33 : min_alpha = 1d0
638 33 : nz = s% nz
639 :
640 33 : if (s% RSP2_flag) & ! clip change in w to maintain non-negativity.
641 0 : call clip_so_non_negative(s% i_w, 0d0)
642 :
643 33 : if (s% RTI_flag) & ! clip change in alpha_RTI to maintain non-negativity.
644 0 : call clip_so_non_negative(s% i_alpha_RTI, 0d0)
645 :
646 33 : if (s% i_lum>=0 .and. s% scale_max_correction_for_negative_surf_lum) then
647 : !ensure surface luminosity does not become negative
648 0 : dlum_surf = B(s% i_lum,1)*s% x_scale(s% i_lum,1)
649 0 : old_lum_surf = s% xh_start(s% i_lum,1) + s% solver_dx(s% i_lum,1)
650 0 : new_lum_surf = old_lum_surf + dlum_surf
651 0 : if (new_lum_surf < 0d0 .and. old_lum_surf > 0d0) then
652 : correction_factor = min(correction_factor, &
653 0 : -s% max_frac_for_negative_surf_lum*old_lum_surf/dlum_surf)
654 : end if
655 : end if
656 :
657 : !if (s% w_div_wc_flag) then
658 : ! do k=1,nz
659 : ! old_xa = s% xh_start(s% i_w_div_wc,k) + dx(s% i_w_div_wc,k)
660 : ! dxa = B(s% i_w_div_wc,k)*xscale(s% i_w_div_wc,k)*correction_factor
661 : ! if (dxa > 0.9d0*(1-old_xa)) then
662 : ! if(k==91)write(*,*) "problems", old_xa, dxa, old_xa+dxa
663 : ! correction_factor = min(correction_factor,(0.9d0*(1-old_xa))/B(s% i_w_div_wc,k)*xscale(s%i_w_div_wc,k))
664 : ! dxa = B(s% i_w_div_wc,k)*xscale(s% i_w_div_wc,k)*correction_factor
665 : ! if(k==91)write(*,*) "need to reduce correction", k, old_xa+dxa
666 : ! end if
667 : ! end do
668 : !end if
669 :
670 33 : if ((.not. s% do_solver_damping_for_neg_xa) .or. &
671 : (.not. (s% do_mix .or. s% do_burn))) then
672 : if (min_alpha < 1d0) then
673 : min_alpha = max(min_alpha, s% corr_coeff_limit)
674 : correction_factor = min_alpha*correction_factor
675 : end if
676 : return
677 : end if
678 :
679 33 : bad_j = 0
680 33 : bad_k = 0
681 33 : if (nvar > s% nvar_hydro) then
682 33 : species = s% species
683 33 : min_xa_hard_limit = current_min_xa_hard_limit(s)
684 33 : eps = -0.5d0*min_xa_hard_limit ! allow xa to be this much below 0d0
685 39294 : do k=1,nz
686 353382 : do j=1,species
687 314088 : i = j + s% nvar_hydro
688 314088 : old_xa = s% xa_start(j,k) + s% solver_dx(i,k)
689 314088 : if (old_xa <= 1d-90) cycle
690 314088 : dxa = B(i,k)*s% x_scale(i,k)*correction_factor
691 314088 : new_xa = old_xa + dxa
692 314088 : if (new_xa >= 0d0) cycle
693 0 : alpha = -(old_xa + eps)/dxa
694 : ! so dxa*alpha = -old_xa - eps,
695 : ! and therefore old_xa + alpha*dxa = -eps = 0.5*min_xa_hard_limit
696 39261 : if (alpha < min_alpha) then
697 0 : min_alpha = alpha
698 0 : bad_j = j
699 0 : bad_k = k
700 : end if
701 : end do
702 : end do
703 : end if
704 :
705 33 : min_alpha = max(min_alpha, s% corr_coeff_limit)
706 33 : correction_factor = min_alpha*correction_factor
707 33 : if (s% trace_solver_damping .and. min_alpha < 1d0 .and. bad_j > 0) then
708 : write(*,4) 'solver damping to avoid negative mass fractions: ' // &
709 0 : trim(chem_isos% name(s% chem_id(bad_j))), bad_k, &
710 0 : s% model_number, s% solver_iter, min_alpha
711 : end if
712 :
713 : contains
714 :
715 0 : subroutine clip_so_non_negative(i,minval)
716 : integer, intent(in) :: i
717 : real(dp), intent(in) :: minval
718 0 : real(dp) :: dval, old_val, new_val
719 0 : do k = 1, s% nz
720 0 : dval = B(i,k)*s% x_scale(i,k)*correction_factor
721 0 : old_val = s% xh_start(i,k) + s% solver_dx(i,k)
722 0 : new_val = old_val + dval
723 0 : if (dval >= 0) cycle
724 0 : if (new_val >= 0d0) cycle
725 0 : dval = minval - old_val
726 0 : B(i,k) = dval/(s% x_scale(i,k)*correction_factor)
727 : end do
728 33 : end subroutine clip_so_non_negative
729 :
730 : end subroutine Bdomain
731 :
732 :
733 33 : subroutine inspectB(s, nvar, B, ierr)
734 : type (star_info), pointer :: s
735 : integer, intent(in) :: nvar
736 : real(dp), pointer, dimension(:,:) :: B ! (nvar, nz)
737 : integer, intent(out) :: ierr
738 :
739 : integer, parameter :: inspectB_iter_stop = -1
740 : include 'formats'
741 :
742 : if (dbg) write(*, *) 'inspectB', s% solver_iter
743 33 : ierr = 0
744 33 : if (s% solver_iter == inspectB_iter_stop) then
745 0 : call dumpB
746 0 : call mesa_error(__FILE__,__LINE__,'debug: inspectB')
747 : end if
748 :
749 : contains
750 :
751 0 : subroutine dumpB
752 : integer :: k, j
753 : include 'formats'
754 0 : do k=1,s% nz
755 0 : do j=1,nvar
756 0 : write(*,2) 'B ' // trim(s% nameofvar(j)), k, B(j, k)
757 0 : write(*,2) 'xscale ' // trim(s% nameofvar(j)), k, s% x_scale(j, k)
758 0 : write(*,2) 'dx ' // trim(s% nameofvar(j)), k, s% solver_dx(j, k)
759 : end do
760 0 : write(*,'(A)')
761 : end do
762 0 : call mesa_error(__FILE__,__LINE__,'dumpB')
763 0 : end subroutine dumpB
764 :
765 : end subroutine inspectB
766 :
767 :
768 : ! about to declare victory... but may want to do another iteration
769 : ! 1 means force another iteration
770 : ! 0 means don't need to force another
771 : ! -1 means failure. solver returns with non-convergence.
772 11 : integer function force_another_iteration(s, iter, itermin)
773 : type (star_info), pointer :: s
774 : integer, intent(in) :: iter ! have finished this many iterations and have converged
775 : integer, intent(in) :: itermin ! this is the requested minimum. iter may be < itermin.
776 :
777 : include 'formats'
778 :
779 11 : if (iter < itermin) then
780 0 : force_another_iteration = 1
781 : return
782 : end if
783 11 : force_another_iteration = 0
784 :
785 : if (s% k_below_just_added > 1 .and. &
786 11 : s% num_surf_revisions < s% max_num_surf_revisions .and. &
787 : abs(s% lnS(1) - s% surf_lnS) > &
788 : s% max_abs_rel_change_surf_lnS*max(s% lnS(1),s% surf_lnS)) then
789 0 : s% surf_lnS = s% lnS(1)
790 0 : s% num_surf_revisions = s% num_surf_revisions + 1
791 0 : force_another_iteration = 1
792 0 : s% used_extra_iter_in_solver_for_accretion = .true.
793 0 : return
794 : end if
795 :
796 : end function force_another_iteration
797 :
798 :
799 33 : subroutine set_solver_vars(s, nvar, dt, ierr)
800 : type (star_info), pointer :: s
801 : integer, intent(in) :: nvar
802 : real(dp), intent(in) :: dt
803 : integer, intent(out) :: ierr
804 33 : s% num_solver_setvars = s% num_solver_setvars + 1
805 33 : call set_vars_for_solver(s, nvar, 1, s% nz, dt, ierr)
806 33 : end subroutine set_solver_vars
807 :
808 :
809 33 : subroutine set_vars_for_solver(s, nvar, nzlo, nzhi, dt, ierr)
810 : use const_def, only: secyer, Msun, Lsun, Rsun
811 : use star_utils, only: set_rmid, set_dm_bar, set_m_and_dm, set_rv_info
812 : use star_utils, only: current_min_xa_hard_limit, current_sum_xa_hard_limit, &
813 : lookup_nameofvar
814 : use hydro_vars, only: set_hydro_vars
815 : use hydro_rotation, only: set_i_rot, set_omega
816 : use mix_info, only: get_convection_sigmas
817 : use chem_def
818 : type (star_info), pointer :: s
819 : integer, intent(in) :: nvar, nzlo, nzhi
820 : real(dp), intent(in) :: dt
821 : integer, intent(out) :: ierr
822 :
823 : logical, parameter :: &
824 : skip_basic_vars = .true., &
825 : skip_micro_vars = .false., &
826 : skip_kap = .false., &
827 : skip_neu = .false., &
828 : skip_net = .false., &
829 : skip_eos = .false., &
830 : skip_mlt = .false., &
831 : skip_grads = .true., &
832 : skip_rotation = .true., &
833 : skip_m_grav_and_grav = .true., &
834 : skip_brunt = .true., &
835 : skip_mixing_info = .true., &
836 : skip_set_cz_bdy_mass = .true., &
837 : skip_other_cgrav = .true.
838 : logical :: do_chem, try_again, do_edit_lnR, report_dx
839 : integer :: j, k, kk, klo, khi, i_var, &
840 : i_lnd, i_lnT, i_lnR, i_lum, i_w, i_Hp, i_v, &
841 : i_u, i_alpha_RTI, i_w_div_wc, i_j_rot, &
842 : fe56, nvar_chem, species, nz, nvar_hydro
843 33 : real(dp), dimension(:, :), pointer :: xh_start, xa_start
844 : integer :: op_err, kbad, &
845 : cnt, max_fixes, loc(2), k_lo, k_hi
846 33 : real(dp) :: xavg, dx_for_i_var, x_for_i_var, &
847 33 : dq_sum, d_dxdt_dx, min_xa_hard_limit, sum_xa_hard_limit
848 : logical :: do_lnd, do_lnT, do_lnR, do_lum, do_w, &
849 : do_u, do_v, do_alpha_RTI, do_w_div_wc, do_j_rot
850 :
851 : include 'formats'
852 :
853 33 : ierr = 0
854 :
855 33 : nz = s% nz
856 33 : nvar_chem = s% nvar_chem
857 33 : species = s% species
858 33 : nvar_hydro = s% nvar_hydro
859 33 : d_dxdt_dx = 1d0/s% dt
860 :
861 33 : xh_start => s% xh_start
862 33 : xa_start => s% xa_start
863 :
864 : report_dx = &
865 : s% solver_test_partials_dx_0 > 0d0 .and. &
866 : s% solver_test_partials_k > 0 .and. &
867 : s% solver_call_number == s% solver_test_partials_call_number .and. &
868 : s% solver_test_partials_iter_number == s% solver_iter .and. &
869 0 : len_trim(s% solver_test_partials_show_dx_var_name) > 0
870 :
871 : if (report_dx) then
872 0 : k = s% solver_test_partials_k
873 0 : i_var = lookup_nameofvar(s, s% solver_test_partials_show_dx_var_name)
874 0 : if (i_var > 0) then
875 0 : if (i_var > nvar_hydro) then
876 0 : dx_for_i_var = s% solver_dx(i_var,k)
877 0 : x_for_i_var = xa_start(i_var-nvar_hydro,k) + s% solver_dx(i_var,k)
878 : else
879 0 : dx_for_i_var = s% solver_dx(i_var,k)
880 0 : x_for_i_var = xh_start(i_var,k) + s% solver_dx(i_var,k)
881 : end if
882 : write(*,3) 'dx, x for var name ' // &
883 0 : trim(s% solver_test_partials_show_dx_var_name), &
884 0 : k, s% solver_iter, dx_for_i_var, x_for_i_var
885 : end if
886 : end if
887 :
888 33 : do_chem = (s% do_burn .or. s% do_mix)
889 :
890 : if (dbg .and. .not. skip_mixing_info) write(*,2) 'redo mix info iter', s% solver_iter
891 :
892 33 : min_xa_hard_limit = current_min_xa_hard_limit(s)
893 33 : sum_xa_hard_limit = current_sum_xa_hard_limit(s)
894 :
895 33 : i_lnd = s% i_lnd
896 33 : i_lnT = s% i_lnT
897 33 : i_lnR = s% i_lnR
898 33 : i_lum = s% i_lum
899 33 : i_w = s% i_w
900 33 : i_Hp = s% i_Hp
901 33 : i_v = s% i_v
902 33 : i_u = s% i_u
903 33 : i_alpha_RTI = s% i_alpha_RTI
904 33 : i_w_div_wc = s% i_w_div_wc
905 33 : i_j_rot = s% i_j_rot
906 :
907 33 : do_lnd = i_lnd > 0 .and. i_lnd <= nvar
908 33 : do_lnT = i_lnT > 0 .and. i_lnT <= nvar
909 33 : do_lnR = i_lnR > 0 .and. i_lnR <= nvar
910 33 : do_lum = i_lum > 0 .and. i_lum <= nvar
911 33 : do_w = i_w > 0 .and. i_w <= nvar
912 33 : do_v = i_v > 0 .and. i_v <= nvar
913 33 : do_u = i_u > 0 .and. i_u <= nvar
914 33 : do_alpha_RTI = i_alpha_RTI > 0 .and. i_alpha_RTI <= nvar
915 33 : do_w_div_wc = i_w_div_wc > 0 .and. i_w_div_wc <= nvar
916 33 : do_j_rot = i_j_rot > 0 .and. i_j_rot <= nvar
917 :
918 33 : fe56 = s% net_iso(ife56)
919 : if (fe56 /= 0) fe56 = nvar_hydro+fe56
920 :
921 33 : if (nvar > nvar_hydro) then
922 39294 : do k=1,nz
923 353382 : do j=1,species
924 314088 : s% xa_sub_xa_start(j,k) = s% solver_dx(j+nvar_hydro,k)
925 353349 : s% xa(j,k) = xa_start(j,k) + s% solver_dx(j+nvar_hydro,k)
926 : end do
927 : end do
928 : max_fixes = 0 !5
929 : do cnt=1,max_fixes
930 : loc = minloc(s% xa(1:species,1:nz))
931 : j = loc(1)
932 : k = loc(2)
933 : if (s% xa(j,k) >= 1d-3*min_xa_hard_limit) exit ! too good to fix
934 : if (s% xa(j,k) < min_xa_hard_limit) then
935 : if (s% report_ierr) then
936 : khi = nz
937 : do kk=k+1,nz
938 : if (s% xa(j,kk) < min_xa_hard_limit) cycle
939 : khi = kk-1; exit
940 : end do
941 : klo = 1
942 : do kk=k-1,1,-1
943 : if (s% xa(j,kk) < min_xa_hard_limit) cycle
944 : klo = kk+1; exit
945 : end do
946 : do k=klo,khi
947 : write(*,2) &
948 : 'negative ' // trim(chem_isos% name(s% chem_id(j))), &
949 : k, s% xa(j,k), xa_start(j,k), s% solver_dx(nvar_hydro+j,k), s% m(k)/Msun
950 : end do
951 : end if
952 : s% retry_message = 'some abundance < min_xa_hard_limit'
953 : ierr = -1
954 : return
955 : exit ! too bad to fix
956 : end if
957 : if (k == 1) then
958 : k_lo = 1; k_hi = 2
959 : else if (k == nz) then
960 : k_lo = nz-1; k_hi = nz
961 : else if (s% sig(k) > s% sig(k+1)) then
962 : k_lo = k-1; k_hi = k
963 : else if (s% sig(k+1) > 0) then
964 : k_lo = k; k_hi = k+1
965 : else
966 : exit
967 : end if
968 : try_again = .true.
969 : do while (try_again .and. sum(s% xa(j,k_lo:k_hi)*s% dq(k_lo:k_hi)) < 0)
970 : try_again = .false.
971 : if (k_lo > 1) then
972 : if (s% sig(k_lo) > 0) then
973 : k_lo = k_lo-1
974 : try_again = .true.
975 : end if
976 : end if
977 : if (k_hi < nz) then
978 : if (s% sig(k_hi+1) > 0) then
979 : k_hi = k_hi+1
980 : try_again = .true.
981 : end if
982 : end if
983 : end do
984 : !write(*,3) 'no extend', k_lo, k_hi
985 : if (.not. try_again) exit
986 : dq_sum = sum(s% dq(k_lo:k_hi))
987 : if (s% report_ierr) then
988 : write(*,5) 'fix xa(j,k_lo:k_hi)', j, k_lo, k, k_hi, s% xa(j,k), &
989 : sum(s% xa(j,k_lo:k_hi)*s% dq(k_lo:k_hi))/dq_sum, dq_sum
990 : !stop
991 : end if
992 : do j=1,species
993 : xavg = sum(s% xa(j,k_lo:k_hi)*s% dq(k_lo:k_hi))/dq_sum
994 : do kk=k_lo,k_hi
995 : s% xa(j,kk) = xavg
996 : end do
997 : end do
998 : end do
999 : end if
1000 :
1001 33 : if (s% solver_test_partials_k > 0 .and. &
1002 : s% solver_test_partials_k <= nz) then
1003 0 : k = s% solver_test_partials_k
1004 : end if
1005 :
1006 33 : kbad = 0
1007 33 : !$OMP PARALLEL DO PRIVATE(k,op_err) SCHEDULE(dynamic,2)
1008 : do k=1,nz
1009 : op_err = 0
1010 : call set1(k,.false.,op_err)
1011 : if (op_err /= 0) then
1012 : kbad = k; ierr = op_err
1013 : end if
1014 : end do
1015 : !$OMP END PARALLEL DO
1016 :
1017 33 : if (ierr /= 0) then
1018 0 : if (s% report_ierr) then
1019 0 : do k=1,nz ! report the errors sequentially
1020 0 : call set1(k,.true.,op_err)
1021 : end do
1022 0 : write(*,3) 'set_vars_for_solver failed: model, nz', &
1023 0 : s% model_number, nz ! kbad depends on num threads
1024 : end if
1025 0 : return
1026 : end if
1027 :
1028 33 : if (s% solver_test_partials_k > 0 .and. &
1029 : s% solver_test_partials_k <= nz) then
1030 0 : k = s% solver_test_partials_k
1031 : end if
1032 :
1033 33 : if (do_lum) then
1034 39294 : do k=1,nz
1035 39294 : if (is_bad_num(s% L(k))) then
1036 0 : if (s% report_ierr) write(*,2) 'set_vars_for_solver L', k, s% L(k), &
1037 0 : xh_start(i_lum,k) + s% solver_dx(i_lum,k), &
1038 0 : xh_start(i_lum,k), s% solver_dx(i_lum,k)
1039 0 : s% retry_message = 'bad num for some L'
1040 0 : ierr = -1
1041 0 : if (s% stop_for_bad_nums) then
1042 0 : write(*,2) 'set_vars_for_solver L', k, s% L(k)
1043 0 : call mesa_error(__FILE__,__LINE__,'set_vars_for_solver')
1044 : end if
1045 0 : return
1046 : stop
1047 : end if
1048 : end do
1049 : end if
1050 :
1051 33 : if (ierr /= 0) then
1052 0 : if (s% report_ierr) then
1053 :
1054 0 : do k=1,nz
1055 0 : if (abs(1d0 - sum(s% xa(:,k))) > 1d-3) then
1056 0 : write(*,2) 'set_vars_for_solver: bad xa sum', k, &
1057 0 : sum(s% xa(:,k)), sum(xa_start(:,k)), sum(s% solver_dx(nvar_hydro+1:nvar,k))
1058 0 : write(*,'(51x,a)') 'xa, xa_start+dx, xa_start, dx'
1059 0 : do j=1,species
1060 0 : write(*,2) trim(chem_isos% name(s% chem_id(j))), k, &
1061 0 : s% xa(j,k), xa_start(j,k) + s% solver_dx(nvar_hydro+j,k), &
1062 0 : xa_start(j,k), s% solver_dx(nvar_hydro+j,k)
1063 : end do
1064 :
1065 : exit
1066 :
1067 : end if
1068 : end do
1069 0 : write(*,'(A)')
1070 :
1071 : end if
1072 0 : return
1073 : end if
1074 :
1075 33 : do_edit_lnR = do_lnR .and. (.not. s% doing_check_partials)
1076 33 : if (do_edit_lnR) call edit_lnR(s, xh_start, s% solver_dx)
1077 39294 : do k=1,nz
1078 39261 : if (do_edit_lnR) s% r(k) = exp(s% lnR(k))
1079 39261 : call set_rv_info(s,k)
1080 : ! note: m_grav is held constant during solver iterations
1081 39294 : s% grav(k) = s% cgrav(k)*s% m_grav(k)/(s% r(k)*s% r(k))
1082 : end do
1083 :
1084 33 : if (do_lnR) then
1085 33 : call set_rmid(s, 1, nz, ierr)
1086 33 : if (ierr /= 0) then
1087 0 : if (s% report_ierr) write(*, *) 'set_rmid returned ierr', ierr
1088 0 : return
1089 : end if
1090 : end if
1091 :
1092 : if (dbg) write(*, *) 'call set_hydro_vars'
1093 : call set_hydro_vars( &
1094 : s, 1, nz, skip_basic_vars, skip_micro_vars, &
1095 : skip_m_grav_and_grav, skip_eos, skip_net, skip_neu, skip_kap, &
1096 : skip_grads, skip_rotation, skip_brunt, skip_other_cgrav, &
1097 33 : skip_mixing_info, skip_set_cz_bdy_mass, skip_mlt, ierr)
1098 33 : if (ierr /= 0) then
1099 0 : if (s% report_ierr) write(*,2) 'set_hydro_vars returned ierr', ierr
1100 0 : return
1101 : end if
1102 :
1103 66 : if (s% rotation_flag) then
1104 : ! Moments of inertia are kept fixed as those at the start of the hydro solver
1105 : ! neither omege nor irot enter the equations directly
1106 : ! call set_i_rot(s)
1107 0 : call set_omega(s, 'hydro_mtx')
1108 : end if
1109 :
1110 :
1111 : contains
1112 :
1113 :
1114 39261 : subroutine set1(k,report,ierr)
1115 : !use chem_def, only: chem_isos
1116 : integer, intent(in) :: k
1117 : logical, intent(in) :: report
1118 : integer, intent(out) :: ierr
1119 :
1120 510393 : real(dp) :: x(nvar)
1121 : ! setting x = xh_start + dx is necessary because of numerical issues.
1122 : ! we want to ensure that we've calculated the variables using exactly the
1123 : ! same values for x as will be returned as the final result.
1124 : integer :: j, k_below_just_added
1125 39261 : real(dp) :: v
1126 :
1127 : include 'formats'
1128 39261 : ierr = 0
1129 39261 : v = 0
1130 :
1131 39261 : k_below_just_added = 1
1132 :
1133 196305 : do j=1,min(nvar, nvar_hydro)
1134 196305 : x(j) = xh_start(j,k) + s% solver_dx(j,k)
1135 : !write(*,2) 'new ' // s% nameofvar(j), k, x(j)
1136 : end do
1137 :
1138 39261 : if (do_lnT) then
1139 :
1140 39261 : s% lnT(k) = x(i_lnT)
1141 39261 : s% T(k) = exp(s% lnT(k))
1142 39261 : s% dxh_lnT(k) = s% solver_dx(i_lnT,k)
1143 : if (abs(s% lnT(k) - s% lnT_start(k)) > &
1144 39261 : ln10*s% hydro_mtx_max_allowed_abs_dlogT .and. &
1145 : s% min_logT_for_hydro_mtx_max_allowed < &
1146 : ln10*min(s% lnT(k),s% lnT_start(k))) then
1147 0 : if (report) &
1148 0 : write(*,4) 'hydro_mtx: change too large, dlogT, logT, logT_start', &
1149 0 : s% model_number, k, s% solver_iter, &
1150 0 : (s% lnT(k) - s% lnT_start(k))/ln10, &
1151 0 : s% lnT(k)/ln10, s% lnT_start(k)/ln10
1152 0 : write(s% retry_message, *) 'abs(dlogT) > hydro_mtx_max_allowed_abs_dlogT', k
1153 0 : ierr = -1
1154 0 : return
1155 : end if
1156 39261 : if (s% lnT(k) > ln10*s% hydro_mtx_max_allowed_logT .and. &
1157 : s% min_logT_for_hydro_mtx_max_allowed < &
1158 : ln10*min(s% lnT(k),s% lnT_start(k))) then
1159 0 : if (report) &
1160 0 : write(*,4) 'hydro_mtx: logT too large', &
1161 0 : s% model_number, k, s% solver_iter, &
1162 0 : s% lnT(k)/ln10, s% lnT_start(k)/ln10
1163 0 : write(s% retry_message, *) 'logT > hydro_mtx_max_allowed_logT', k
1164 0 : ierr = -1
1165 0 : return
1166 : end if
1167 39261 : if (s% lnT(k) < ln10*s% hydro_mtx_min_allowed_logT) then
1168 0 : if (report) &
1169 0 : write(*,4) 'hydro_mtx: logT too small', &
1170 0 : s% model_number, k, s% solver_iter, &
1171 0 : s% lnT(k)/ln10, s% lnT_start(k)/ln10
1172 0 : write(s% retry_message, *) 'logT < hydro_mtx_min_allowed_logT', k
1173 0 : ierr = -1
1174 0 : return
1175 : end if
1176 39261 : s% T(k) = exp(s% lnT(k))
1177 39261 : if (is_bad_num(s% T(k))) then
1178 0 : s% retry_message = 'bad num for T'
1179 0 : if (s% stop_for_bad_nums) then
1180 0 : write(*,2) 'set_vars_for_solver T', k, s% T(k)
1181 0 : call mesa_error(__FILE__,__LINE__,'set_vars_for_solver')
1182 : end if
1183 0 : if (report) write(*,2) 'bad num T', k, s% T(k)
1184 0 : ierr = -1
1185 : end if
1186 :
1187 : end if
1188 :
1189 39261 : if (do_lum) then
1190 39261 : s% L(k) = x(i_lum)
1191 39261 : if (is_bad_num(s% L(k))) then
1192 0 : s% retry_message = 'bad num for L'
1193 0 : ierr = -1
1194 0 : if (s% stop_for_bad_nums) then
1195 0 : write(*,2) 'set_vars_for_solver L', k, s% L(k)
1196 0 : call mesa_error(__FILE__,__LINE__,'set_vars_for_solver')
1197 : end if
1198 0 : if (report) write(*,2) 'bad num L', k, s% L(k)
1199 : end if
1200 : end if
1201 :
1202 39261 : if (do_w) then
1203 0 : s% w(k) = x(i_w)
1204 0 : if (s% w(k) < 0d0) then
1205 : !write(*,4) 'set_vars_for_solver: fix w < 0', k, &
1206 : ! s% solver_iter, s% model_number, s% w(k)
1207 0 : s% w(k) = s% RSP2_w_fix_if_neg
1208 : end if
1209 0 : if (is_bad_num(s% w(k))) then
1210 0 : s% retry_message = 'bad num for w'
1211 0 : ierr = -1
1212 0 : if (s% stop_for_bad_nums) then
1213 0 : write(*,2) 'set_vars_for_solver w', k, s% w(k)
1214 0 : call mesa_error(__FILE__,__LINE__,'set_vars_for_solver')
1215 : end if
1216 0 : if (report) write(*,2) 'bad num w', k, s% w(k)
1217 : end if
1218 0 : s% Hp_face(k) = x(i_Hp)
1219 0 : if (is_bad_num(s% Hp_face(k))) then
1220 0 : s% retry_message = 'bad num for Hp_face'
1221 0 : ierr = -1
1222 0 : if (s% stop_for_bad_nums) then
1223 0 : write(*,2) 'set_vars_for_solver Hp_face', k, s% Hp_face(k)
1224 0 : call mesa_error(__FILE__,__LINE__,'set_vars_for_solver')
1225 : end if
1226 0 : if (report) write(*,2) 'bad num Hp_face', k, s% Hp_face(k)
1227 : end if
1228 : end if
1229 :
1230 39261 : if (do_v) then
1231 0 : s% v(k) = x(i_v)
1232 0 : s% dxh_v(k) = s% solver_dx(i_v,k)
1233 0 : if (is_bad_num(s% v(k))) then
1234 0 : s% retry_message = 'bad num for v'
1235 0 : if (s% stop_for_bad_nums) then
1236 0 : write(*,2) 'set_vars_for_solver v', k, s% v(k)
1237 0 : call mesa_error(__FILE__,__LINE__,'set_vars_for_solver')
1238 : end if
1239 0 : if (report) write(*,2) 'bad num v', k, s% v(k)
1240 0 : ierr = -1
1241 : end if
1242 : end if
1243 :
1244 39261 : if (do_u) then
1245 0 : s% u(k) = x(i_u)
1246 0 : s% dxh_u(k) = s% solver_dx(i_u,k)
1247 0 : if (is_bad_num(s% u(k))) then
1248 0 : s% retry_message = 'bad num for u'
1249 0 : if (s% stop_for_bad_nums) then
1250 0 : write(*,2) 'set_vars_for_solver u', k, s% u(k)
1251 0 : call mesa_error(__FILE__,__LINE__,'set_vars_for_solver')
1252 : end if
1253 0 : if (report) write(*,2) 'bad num u', k, s% u(k)
1254 0 : ierr = -1
1255 : end if
1256 : end if
1257 :
1258 39261 : if (do_alpha_RTI) then
1259 0 : s% alpha_RTI(k) = max(0d0, x(i_alpha_RTI))
1260 0 : s% dxh_alpha_RTI(k) = s% solver_dx(i_alpha_RTI,k)
1261 0 : if (is_bad_num(s% alpha_RTI(k))) then
1262 0 : s% retry_message = 'bad num for alpha_RTI'
1263 0 : if (s% stop_for_bad_nums) then
1264 0 : write(*,2) 'set_vars_for_solver alpha_RTI', k, s% alpha_RTI(k)
1265 0 : call mesa_error(__FILE__,__LINE__,'set_vars_for_solver')
1266 : end if
1267 0 : if (report) write(*,2) 'bad num alpha_RTI', k, s% alpha_RTI(k)
1268 0 : ierr = -1
1269 : end if
1270 0 : if (s% alpha_RTI(k) < 0d0) then
1271 0 : s% alpha_RTI(k) = 0d0
1272 0 : x(i_alpha_RTI) = 0d0
1273 : end if
1274 : end if
1275 :
1276 39261 : if (do_w_div_wc) then
1277 0 : s% w_div_w_crit_roche(k) = x(i_w_div_wc)
1278 0 : if (s% w_div_w_crit_roche(k) > 0.99d0) then
1279 0 : s% w_div_w_crit_roche(k) = 0.99d0
1280 : end if
1281 0 : if (s% w_div_w_crit_roche(k) < -0.99d0) then
1282 0 : s% w_div_w_crit_roche(k) = -0.99d0
1283 : end if
1284 0 : if (is_bad_num(s% w_div_w_crit_roche(k))) then
1285 0 : s% retry_message = 'bad num for w_div_w_crit_roche'
1286 0 : if (s% stop_for_bad_nums) then
1287 0 : write(*,2) 'set_vars_for_solver w_div_w_crit_roche', k, s% w_div_w_crit_roche(k)
1288 0 : call mesa_error(__FILE__,__LINE__,'set_vars_for_solver')
1289 : end if
1290 0 : if (report) write(*,2) 'bad num w_div_w_crit_roche', k, s% w_div_w_crit_roche(k)
1291 0 : ierr = -1
1292 : end if
1293 : end if
1294 :
1295 39261 : if (do_j_rot) then
1296 0 : s% j_rot(k) = x(i_j_rot)
1297 0 : if (is_bad_num(s% j_rot(k))) then
1298 0 : s% retry_message = 'bad num for j_rot'
1299 0 : if (s% stop_for_bad_nums) then
1300 0 : write(*,2) 'set_vars_for_solver j_rot', k, s% j_rot(k)
1301 0 : call mesa_error(__FILE__,__LINE__,'set_vars_for_solver')
1302 : end if
1303 0 : if (report) write(*,2) 'bad num j_rot', k, s% j_rot(k)
1304 0 : ierr = -1
1305 : end if
1306 : end if
1307 :
1308 39261 : if (do_lnR) then
1309 39261 : s% lnR(k) = x(i_lnR)
1310 39261 : s% r(k) = exp(s% lnR(k))
1311 39261 : s% dxh_lnR(k) = s% solver_dx(i_lnR,k)
1312 39261 : if (is_bad_num(s% lnR(k))) then
1313 0 : s% retry_message = 'bad num for lnR'
1314 0 : if (s% stop_for_bad_nums) then
1315 0 : write(*,2) 'set_vars_for_solver r lnR solver_dx', &
1316 0 : k, s% r(k), s% lnR(k), s% solver_dx(i_lnR,k)
1317 0 : call mesa_error(__FILE__,__LINE__,'set_vars_for_solver')
1318 : end if
1319 0 : if (report) write(*,2) 'bad num r lnR solver_dx', &
1320 0 : k, s% r(k), s% lnR(k), s% solver_dx(i_lnR,k)
1321 0 : ierr = -1
1322 : end if
1323 : end if
1324 :
1325 39261 : if (do_lnd) then
1326 39261 : s% lnd(k) = x(i_lnd)
1327 39261 : s% rho(k) = exp(s% lnd(k))
1328 39261 : s% dxh_lnd(k) = s% solver_dx(i_lnd,k)
1329 39261 : if (s% lnd(k) < ln10*s% hydro_mtx_min_allowed_logRho) then
1330 0 : write(s% retry_message, *) 'logRho < hydro_mtx_min_allowed_logRho', k
1331 0 : if (report) &
1332 0 : write(*,4) 'hydro_mtx: logRho too small', &
1333 0 : s% model_number, k, s% solver_iter, &
1334 0 : s% lnd(k)/ln10, s% lnd_start(k)/ln10
1335 0 : ierr = -1
1336 0 : return
1337 : end if
1338 39261 : if (s% lnd(k) > ln10*s% hydro_mtx_max_allowed_logRho .and. &
1339 : s% min_logT_for_hydro_mtx_max_allowed < &
1340 : ln10*min(s% lnT(k),s% lnT_start(k))) then
1341 0 : write(s% retry_message, *) 'logRho > hydro_mtx_max_allowed_logRho', k
1342 0 : if (s% report_ierr .or. report) &
1343 0 : write(*,4) 'hydro_mtx: logRho too large', &
1344 0 : s% model_number, k, s% solver_iter, &
1345 0 : s% lnd(k)/ln10, s% lnd_start(k)/ln10
1346 0 : ierr = -1
1347 0 : return
1348 : end if
1349 : if (abs(s% lnd(k) - s% lnd_start(k)) > &
1350 39261 : ln10*s% hydro_mtx_max_allowed_abs_dlogRho .and. &
1351 : s% min_logT_for_hydro_mtx_max_allowed < &
1352 : ln10*min(s% lnT(k),s% lnT_start(k))) then
1353 0 : write(s% retry_message, *) 'abs(dlogRho) > hydro_mtx_max_allowed_abs_dlogRho', k
1354 0 : if (s% report_ierr .or. report) &
1355 : write(*,4) &
1356 0 : 'hydro_mtx: dlogRho, logRho, logRho_start', &
1357 0 : s% model_number, k, s% solver_iter, &
1358 0 : (s% lnd(k) - s% lnd_start(k))/ln10, &
1359 0 : s% lnd(k)/ln10, s% lnd_start(k)/ln10
1360 0 : ierr = -1
1361 0 : return
1362 : end if
1363 39261 : if (is_bad_num(s% rho(k))) then
1364 0 : write(s% retry_message, *) 'bad num for rho', k
1365 0 : if (s% stop_for_bad_nums) then
1366 0 : write(*,2) 'set_vars_for_solver rho', k, s% rho(k)
1367 0 : call mesa_error(__FILE__,__LINE__,'set_vars_for_solver')
1368 : end if
1369 0 : if (report) write(*,2) 'bad num rho', k, s% rho(k)
1370 0 : ierr = -1
1371 : end if
1372 : end if
1373 :
1374 39261 : if (do_chem) &
1375 : call check1_chem( &
1376 39261 : s, k, min_xa_hard_limit, sum_xa_hard_limit, report, ierr)
1377 :
1378 39294 : end subroutine set1
1379 :
1380 :
1381 : end subroutine set_vars_for_solver
1382 :
1383 :
1384 39261 : subroutine check1_chem( &
1385 : s, k, min_xa_hard_limit, sum_xa_hard_limit, report, ierr)
1386 : use chem_def, only: chem_isos
1387 : type (star_info), pointer :: s
1388 : integer, intent(in) :: k
1389 : real(dp), intent(in) :: min_xa_hard_limit, sum_xa_hard_limit
1390 : logical, intent(in) :: report
1391 : integer, intent(out) :: ierr
1392 :
1393 : integer :: j, species
1394 39261 : real(dp) :: sum_xa
1395 : logical :: okay
1396 :
1397 : include 'formats'
1398 :
1399 39261 : ierr = 0
1400 39261 : species = s% species
1401 39261 : okay = .true.
1402 39261 : if (min_xa_hard_limit > -1d50) then
1403 353349 : do j=1,species
1404 353349 : if (s% xa(j,k) < min_xa_hard_limit) then
1405 0 : s% retry_message = 'some xa < min_xa_hard_limit'
1406 0 : if (report .or. s% report_bad_negative_xa) then
1407 : write(*,3) &
1408 : 'bad negative xa, min_xa_hard_limit, sig, logT ' // &
1409 0 : trim(chem_isos% name(s% chem_id(j))), j, k, &
1410 0 : s% xa(j,k), min_xa_hard_limit, s% sig(k), s% lnT(k)/ln10
1411 0 : okay = .false.
1412 : end if
1413 0 : ierr = -1
1414 0 : !$omp critical (hydro_mtx_crit1)
1415 0 : s% why_Tlim = Tlim_neg_X
1416 : !$omp end critical (hydro_mtx_crit1)
1417 0 : if (.not. report) return
1418 : end if
1419 : end do
1420 : end if
1421 :
1422 353349 : do j=1,species
1423 353349 : s% xa(j,k) = max(0d0, min(1d0, s% xa(j,k)))
1424 : end do
1425 :
1426 353349 : sum_xa = sum(s% xa(1:species,k))
1427 39261 : if (is_bad_num(sum_xa)) then
1428 0 : ierr = -1
1429 0 : s% retry_message = 'bad num for mass fractions'
1430 0 : if (s% stop_for_bad_nums) then
1431 0 : write(*,2) 'set_vars_for_solver sum_xa', k, sum_xa
1432 0 : call mesa_error(__FILE__,__LINE__,'set_vars_for_solver')
1433 : end if
1434 0 : if (report) then
1435 0 : write(*,'(a60,i8,99f20.10)') 'bad num sum X', k, s% m(k)/Msun, sum_xa
1436 : end if
1437 0 : !$omp critical (hydro_mtx_crit2)
1438 0 : if (s% why_Tlim /= Tlim_neg_X) s% why_Tlim = Tlim_bad_Xsum
1439 : !$omp end critical (hydro_mtx_crit2)
1440 0 : if (.not. report) return
1441 : end if
1442 39261 : if (abs(sum_xa - 1d0) > sum_xa_hard_limit) then
1443 0 : s% retry_message = 'mass fractions > sum_xa_hard_limit'
1444 0 : if (report) then
1445 : write(*,2) &
1446 0 : 'bad sumX', k, &
1447 0 : sum_xa - 1d0, sum_xa_hard_limit, &
1448 0 : sum(s% xa_start(1:species,k)), &
1449 0 : sum_xa - sum(s% xa_start(1:species,k))
1450 : end if
1451 0 : ierr = -1
1452 0 : !$omp critical (hydro_mtx_crit3)
1453 0 : if (s% why_Tlim /= Tlim_neg_X) s% why_Tlim = Tlim_bad_Xsum
1454 : !$omp end critical (hydro_mtx_crit3)
1455 0 : okay = .false.
1456 0 : if (.not. report) return
1457 : end if
1458 :
1459 39261 : if (abs(sum_xa - 1d0) > 1d-12) then
1460 : !jmax = maxloc(s% xa(1:species,k),dim=1)
1461 : !xsum = sum(s% xa(1:species,k)) - s% xa(jmax,k)
1462 : !if (1d0 > xsum) then
1463 : ! s% xa(jmax,k) = 1d0 - xsum
1464 : !else
1465 5598 : do j=1,species
1466 5598 : s% xa(j,k) = s% xa(j,k)/sum_xa
1467 : end do
1468 : !end if
1469 : end if
1470 :
1471 39261 : if (s% xa_clip_limit > 0) then
1472 353349 : do j=1,species
1473 353349 : if (s% xa(j,k) < s% xa_clip_limit) s% xa(j,k) = 0d0
1474 : end do
1475 : end if
1476 :
1477 39261 : end subroutine check1_chem
1478 :
1479 :
1480 : subroutine dump_struct(s)
1481 : type (star_info), pointer :: s
1482 : integer :: k, j
1483 :
1484 : include 'formats'
1485 :
1486 : do k=1,s% nz
1487 : write(*,2) 'dq', k, s% dq(k)
1488 : write(*,2) 'm', k, s% m(k)
1489 : write(*,2) 'T', k, s% T(k)
1490 : write(*,2) 'rho', k, s% rho(k)
1491 : write(*,2) 'Pgas', k, s% Pgas(k)
1492 : write(*,2) 'L', k, s% L(k)
1493 : write(*,2) 'r', k, s% r(k)
1494 : write(*,2) 'grada', k, s% grada(k)
1495 : write(*,2) 'opacity', k, s% opacity(k)
1496 : write(*,2) 'd_opacity_dlnd', k, s% d_opacity_dlnd(k)
1497 : write(*,2) 'd_opacity_dlnT', k, s% d_opacity_dlnT(k)
1498 : write(*,2) 'eps_nuc', k, s% eps_nuc(k)
1499 : write(*,2) 'd_epsnuc_dlnd', k, s% d_epsnuc_dlnd(k)
1500 : write(*,2) 'd_epsnuc_dlnT', k, s% d_epsnuc_dlnT(k)
1501 : write(*,2) 'non_nuc_neu', k, s% non_nuc_neu(k)
1502 : write(*,2) 'd_nonnucneu_dlnd', k, s% d_nonnucneu_dlnd(k)
1503 : write(*,2) 'd_nonnucneu_dlnT', k, s% d_nonnucneu_dlnT(k)
1504 : write(*,2) 'eps_grav', k, s% eps_grav_ad(k)% val
1505 : write(*,2) 'gradT', k, s% gradT(k)
1506 : !write(*,2) '', k, s% (k)
1507 : do j=1,s% species
1508 : write(*,3) 'xa(j,k)', j, k, s% xa(j,k)
1509 : end do
1510 : end do
1511 :
1512 :
1513 39261 : end subroutine dump_struct
1514 :
1515 :
1516 33 : subroutine edit_lnR(s, xh_start, dx)
1517 : ! uses mass and density to set radius
1518 : type (star_info), pointer :: s
1519 : real(dp), dimension(:, :) :: xh_start
1520 : real(dp), dimension(:, :) :: dx
1521 33 : real(dp) :: vol00, volp1, cell_vol
1522 : integer :: k, nz
1523 : include 'formats'
1524 33 : vol00 = four_thirds_pi*s% R_center*s% R_center*s% R_center
1525 33 : nz = s% nz
1526 39294 : do k=nz, 1, -1
1527 39261 : volp1 = vol00
1528 39261 : cell_vol = s% dm(k)/s% rho(k)
1529 39261 : vol00 = volp1 + cell_vol
1530 39261 : s% lnR(k) = log(vol00/four_thirds_pi)/3
1531 39294 : dx(s% i_lnR,k) = s% lnR(k) - xh_start(s% i_lnR,k)
1532 : end do
1533 33 : call edit_dlnR_dt_above_k_below_just_added(s, xh_start)
1534 33 : end subroutine edit_lnR
1535 :
1536 :
1537 33 : subroutine edit_dlnR_dt_above_k_below_just_added(s, xh_start)
1538 : type (star_info), pointer :: s
1539 : real(dp), dimension(:, :) :: xh_start
1540 : integer :: k_below_just_added
1541 : real(dp) :: lnR_start
1542 33 : k_below_just_added = s% k_below_just_added
1543 : if (k_below_just_added == 1) return
1544 : lnR_start = xh_start(s% i_lnR,1)
1545 : end subroutine edit_dlnR_dt_above_k_below_just_added
1546 :
1547 : end module solver_support
|
