Line data Source code
1 : ! ***********************************************************************
2 : !
3 : ! Copyright (C) 2013-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 :
21 : module star_solver
22 :
23 : use star_private_def
24 : use const_def, only: dp, i8
25 : use num_def
26 : use mtx_def
27 : use mtx_lib, only: block_multiply_xa
28 : use utils_lib, only: is_bad
29 : use solver_support
30 :
31 : implicit none
32 :
33 : private
34 : public :: solver
35 :
36 : contains
37 :
38 11 : subroutine solver( &
39 : s, nvar, skip_global_corr_coeff_limit, &
40 : gold_tolerances_level, tol_max_correction, tol_correction_norm, &
41 : convergence_failure, ierr)
42 : use alloc, only: non_crit_get_quad_array, non_crit_return_quad_array
43 : use utils_lib, only: realloc_if_needed_1, quad_realloc_if_needed_1, fill_with_NaNs
44 :
45 : type (star_info), pointer :: s
46 : ! the primary variables
47 : integer, intent(in) :: nvar ! number of variables per zone
48 : logical, intent(in) :: skip_global_corr_coeff_limit
49 :
50 : ! convergence criteria
51 : integer, intent(in) :: gold_tolerances_level ! 0, 1, or 2
52 : real(dp), intent(in) :: tol_max_correction, tol_correction_norm
53 : ! a trial solution is considered to have converged if
54 : ! max_correction <= tol_max_correction and
55 : !
56 : ! either
57 : ! (correction_norm <= tol_correction_norm)
58 : ! .and. (residual_norm <= tol_residual_norm)
59 : ! or
60 : ! (correction_norm*residual_norm <= tol_corr_resid_product)
61 : ! .and. (abs(slope) <= tol_abs_slope_min)
62 : !
63 : ! where "slope" is slope of the line for line search in the solver,
64 : ! and is analogous to the slope of df/ddx in a 1D solver root finder.
65 :
66 : ! output
67 : logical, intent(out) :: convergence_failure
68 : integer, intent(out) :: ierr ! 0 means okay.
69 :
70 : integer :: ldAF, neqns
71 :
72 : include 'formats'
73 :
74 : ierr = 0
75 :
76 11 : neqns = nvar*s% nz
77 11 : ldAF = 3*nvar
78 :
79 11 : call realloc_if_needed_1(s% AF1,ldAF*neqns,(ldAF+2)*200,ierr)
80 11 : if (ierr /= 0) return
81 :
82 11 : if (s% fill_arrays_with_NaNs) call fill_with_NaNs(s% AF1)
83 :
84 : call do_solver( &
85 : s, nvar, s% AF1, ldAF, neqns, skip_global_corr_coeff_limit, &
86 : gold_tolerances_level, tol_max_correction, tol_correction_norm, &
87 11 : convergence_failure, ierr)
88 :
89 11 : end subroutine solver
90 :
91 :
92 11 : subroutine do_solver( &
93 : s, nvar, AF1, ldAF, neq, skip_global_corr_coeff_limit, &
94 : gold_tolerances_level, tol_max_correction, tol_correction_norm, &
95 : convergence_failure, ierr)
96 :
97 : type (star_info), pointer :: s
98 :
99 : integer, intent(in) :: nvar, ldAF, neq
100 : logical, intent(in) :: skip_global_corr_coeff_limit
101 :
102 : real(dp), pointer, dimension(:) :: AF1 ! =(ldAF, neq)
103 :
104 : ! controls
105 : integer, intent(in) :: gold_tolerances_level
106 : real(dp), intent(in) :: tol_max_correction, tol_correction_norm
107 :
108 : ! output
109 : logical, intent(out) :: convergence_failure
110 : integer, intent(out) :: ierr
111 :
112 : ! info saved in work arrays
113 :
114 : real(dp), dimension(:,:), pointer :: dxsave=>null(), ddxsave=>null(), B=>null(), grad_f=>null(), soln=>null()
115 : real(dp), dimension(:), pointer :: dxsave1=>null(), ddxsave1=>null(), B1=>null(), grad_f1=>null(), &
116 : row_scale_factors1=>null(), col_scale_factors1=>null(), soln1=>null(), &
117 : save_ublk1=>null(), save_dblk1=>null(), save_lblk1=>null()
118 : real(dp), dimension(:,:), pointer :: rhs=>null()
119 : integer, dimension(:), pointer :: ipiv1=>null()
120 : real(dp), dimension(:,:), pointer :: ddx=>null(), xder=>null()
121 :
122 : integer, dimension(:), pointer :: ipiv_blk1=>null()
123 11 : character (len=s%nz) :: equed1
124 :
125 : real(dp), dimension(:), pointer :: lblk1=>null(), dblk1=>null(), ublk1=>null()
126 : real(dp), dimension(:), pointer :: lblkF1=>null(), dblkF1=>null(), ublkF1=>null()
127 :
128 : ! locals
129 : real(dp) :: &
130 11 : coeff, f, slope, residual_norm, max_residual, &
131 11 : corr_norm_min, resid_norm_min, correction_factor, temp_correction_factor, &
132 11 : correction_norm, max_correction, &
133 11 : tol_residual_norm, tol_max_residual, &
134 11 : tol_residual_norm2, tol_max_residual2, &
135 11 : tol_residual_norm3, tol_max_residual3, &
136 11 : tol_abs_slope_min, tol_corr_resid_product, &
137 11 : min_corr_coeff, max_corr_min, max_resid_min, max_abs_correction
138 : integer :: nz, iter, max_tries, tiny_corr_cnt, i, &
139 : force_iter_value, iter_for_resid_tol2, iter_for_resid_tol3, &
140 : max_corr_k, max_corr_j, max_resid_k, max_resid_j
141 : integer(i8) :: time0
142 : character (len=strlen) :: err_msg
143 : logical :: first_try, dbg_msg, passed_tol_tests, &
144 : doing_extra, disabled_resid_tests, pass_resid_tests, &
145 : pass_corr_tests_without_coeff, pass_corr_tests_with_coeff
146 :
147 : integer, parameter :: num_tol_msgs = 15
148 : character (len=32) :: tol_msg(num_tol_msgs)
149 : character (len=64) :: message
150 :
151 : real(dp), pointer, dimension(:) :: equ1=>null()
152 : real(dp), pointer, dimension(:,:) :: equ=>null() ! (nvar,nz)
153 : real(dp), pointer, dimension(:,:) :: AF=>null() ! (ldAF,neq)
154 : real(dp), pointer, dimension(:,:,:) :: ublk=>null(), dblk=>null(), lblk=>null() ! (nvar,nvar,nz)
155 : real(dp), dimension(:,:,:), pointer :: lblkF=>null(), dblkF=>null(), ublkF=>null() ! (nvar,nvar,nz)
156 :
157 11 : call do_solver_work()
158 : ! Split it this way so we can guarantee cleanup() gets called once do_solver_work finishes
159 : ! Otherwise all the pointers will leak memory.
160 22 : call cleanup()
161 :
162 : contains
163 :
164 55 : subroutine do_solver_work
165 : include 'formats'
166 :
167 11 : err_msg = ''
168 :
169 11 : nz = s% nz
170 :
171 11 : AF(1:ldAF,1:neq) => AF1(1:ldAF*neq)
172 :
173 11 : tol_msg(1) = 'avg corr'
174 11 : tol_msg(2) = 'max corr '
175 11 : tol_msg(3) = 'avg+max corr'
176 11 : tol_msg(4) = 'avg resid'
177 11 : tol_msg(5) = 'avg corr+resid'
178 11 : tol_msg(6) = 'max corr, avg resid'
179 11 : tol_msg(7) = 'avg+max corr, avg resid'
180 11 : tol_msg(8) = 'max resid'
181 11 : tol_msg(9) = 'avg corr, max resid'
182 11 : tol_msg(10) = 'max corr+resid'
183 11 : tol_msg(11) = 'avg+max corr, max resid'
184 11 : tol_msg(12) = 'avg+max resid'
185 11 : tol_msg(13) = 'avg corr, avg+max resid'
186 11 : tol_msg(14) = 'max corr, avg+max resid'
187 11 : tol_msg(15) = 'avg+max corr+resid'
188 :
189 11 : ierr = 0
190 11 : iter = 0
191 11 : s% solver_iter = iter
192 :
193 11 : call set_param_defaults
194 11 : dbg_msg = s% report_solver_progress
195 :
196 11 : if (gold_tolerances_level == 2) then
197 0 : tol_residual_norm = s% gold2_tol_residual_norm1
198 0 : tol_max_residual = s% gold2_tol_max_residual1
199 0 : tol_residual_norm2 = s% gold2_tol_residual_norm2
200 0 : tol_max_residual2 = s% gold2_tol_max_residual2
201 0 : tol_residual_norm3 = s% gold2_tol_residual_norm3
202 0 : tol_max_residual3 = s% gold2_tol_max_residual3
203 11 : else if (gold_tolerances_level == 1) then
204 11 : tol_residual_norm = s% gold_tol_residual_norm1
205 11 : tol_max_residual = s% gold_tol_max_residual1
206 11 : tol_residual_norm2 = s% gold_tol_residual_norm2
207 11 : tol_max_residual2 = s% gold_tol_max_residual2
208 11 : tol_residual_norm3 = s% gold_tol_residual_norm3
209 11 : tol_max_residual3 = s% gold_tol_max_residual3
210 : else
211 0 : tol_residual_norm = s% tol_residual_norm1
212 0 : tol_max_residual = s% tol_max_residual1
213 0 : tol_residual_norm2 = s% tol_residual_norm2
214 0 : tol_max_residual2 = s% tol_max_residual2
215 0 : tol_residual_norm3 = s% tol_residual_norm3
216 0 : tol_max_residual3 = s% tol_max_residual3
217 : end if
218 :
219 11 : tol_abs_slope_min = -1 ! unused
220 11 : tol_corr_resid_product = -1 ! unused
221 11 : if (skip_global_corr_coeff_limit) then
222 11 : min_corr_coeff = 1
223 : else
224 0 : min_corr_coeff = s% corr_coeff_limit
225 : end if
226 :
227 11 : if (gold_tolerances_level == 2) then
228 0 : iter_for_resid_tol2 = s% gold2_iter_for_resid_tol2
229 0 : iter_for_resid_tol3 = s% gold2_iter_for_resid_tol3
230 11 : else if (gold_tolerances_level == 1) then
231 11 : iter_for_resid_tol2 = s% gold_iter_for_resid_tol2
232 11 : iter_for_resid_tol3 = s% gold_iter_for_resid_tol3
233 : else
234 0 : iter_for_resid_tol2 = s% iter_for_resid_tol2
235 0 : iter_for_resid_tol3 = s% iter_for_resid_tol3
236 : end if
237 :
238 11 : call pointers(ierr)
239 11 : if (ierr /= 0) then
240 : return
241 : end if
242 :
243 11 : doing_extra = .false.
244 11 : passed_tol_tests = .false. ! goes true when pass the tests
245 11 : convergence_failure = .false. ! goes true when time to give up
246 11 : coeff = 1.d0
247 :
248 11 : residual_norm=0
249 11 : max_residual=0
250 11 : corr_norm_min=1d99
251 11 : max_corr_min=1d99
252 11 : max_resid_min=1d99
253 11 : resid_norm_min=1d99
254 11 : correction_factor=0
255 11 : f=0d0
256 11 : slope=0d0
257 :
258 11 : call set_xscale_info(s, nvar, ierr)
259 11 : if (ierr /= 0) then
260 0 : if (dbg_msg) &
261 0 : write(*, *) 'solver failure: set_xscale_info returned ierr', ierr
262 0 : convergence_failure = .true.
263 0 : return
264 : end if
265 :
266 11 : call do_equations(ierr)
267 11 : if (ierr /= 0) then
268 0 : if (dbg_msg) &
269 0 : write(*, *) 'solver failure: eval_equations returned ierr', ierr
270 0 : convergence_failure = .true.
271 0 : return
272 : end if
273 :
274 11 : call sizequ(s, nvar, residual_norm, max_residual, max_resid_k, max_resid_j, ierr)
275 11 : if (ierr /= 0) then
276 0 : if (dbg_msg) &
277 0 : write(*, *) 'solver failure: sizequ returned ierr', ierr
278 0 : convergence_failure = .true.
279 0 : return
280 : end if
281 :
282 11 : first_try = .true.
283 11 : iter = 1
284 11 : s% solver_iter = iter
285 11 : if (s% doing_first_model_of_run) then
286 1 : max_tries = s% max_tries1
287 10 : else if (s% retry_cnt > 20) then
288 0 : max_tries = s% max_tries_after_20_retries
289 10 : else if (s% retry_cnt > 10) then
290 0 : max_tries = s% max_tries_after_10_retries
291 10 : else if (s% retry_cnt > 5) then
292 0 : max_tries = s% max_tries_after_5_retries
293 10 : else if (s% retry_cnt > 0) then
294 0 : max_tries = s% max_tries_for_retry
295 : else
296 10 : max_tries = s% solver_max_tries_before_reject
297 : end if
298 11 : tiny_corr_cnt = 0
299 :
300 11 : s% num_solver_iterations = 0
301 :
302 44 : iter_loop: do while (.not. passed_tol_tests)
303 :
304 33 : if (dbg_msg .and. first_try) write(*, *)
305 :
306 33 : max_resid_j = -1
307 33 : max_corr_j = -1
308 :
309 33 : if (iter >= iter_for_resid_tol2) then
310 0 : if (iter < iter_for_resid_tol3) then
311 0 : tol_residual_norm = tol_residual_norm2
312 0 : tol_max_residual = tol_max_residual2
313 0 : if (dbg_msg .and. iter == iter_for_resid_tol2) &
314 0 : write(*,1) 'tol2 residual tolerances: norm, max', &
315 0 : tol_residual_norm, tol_max_residual
316 : else
317 0 : tol_residual_norm = tol_residual_norm3
318 0 : tol_max_residual = tol_max_residual3
319 0 : if (dbg_msg .and. iter == iter_for_resid_tol3) &
320 0 : write(*,1) 'tol3 residual tolerances: norm, max', &
321 0 : tol_residual_norm, tol_max_residual
322 : end if
323 33 : else if (dbg_msg .and. iter == 1) then
324 0 : write(*,2) 'solver_call_number', s% solver_call_number
325 0 : write(*,2) 'gold tolerances level', gold_tolerances_level
326 0 : write(*,1) 'correction tolerances: norm, max', &
327 0 : tol_correction_norm, tol_max_correction
328 0 : write(*,1) 'tol1 residual tolerances: norm, max', &
329 0 : tol_residual_norm, tol_max_residual
330 : end if
331 :
332 33 : call solver_test_partials(nvar, xder, ierr)
333 33 : if (ierr /= 0) then
334 0 : call write_msg('solver_test_partials returned ierr /= 0')
335 0 : convergence_failure = .true.
336 0 : exit iter_loop
337 : end if
338 :
339 33 : s% num_solver_iterations = s% num_solver_iterations + 1
340 : if (s% model_number == 1 .and. &
341 33 : s% num_solver_iterations > 60 .and. &
342 : mod(s% num_solver_iterations,10) == 0) &
343 0 : write(*,*) 'first model is slow to converge: num tries', &
344 0 : s% num_solver_iterations
345 :
346 33 : if (.not. solve_equ()) then ! either singular or horribly ill-conditioned
347 0 : write(err_msg, '(a, i5, 3x, a)') 'info', ierr, 'bad_matrix'
348 0 : call oops(err_msg)
349 0 : exit iter_loop
350 : end if
351 :
352 33 : call inspectB(s, nvar, soln, ierr)
353 33 : if (ierr /= 0) then
354 0 : call oops('inspectB returned ierr')
355 0 : exit iter_loop
356 : end if
357 :
358 : ! compute size of scaled correction B
359 : call sizeB(s, nvar, soln, &
360 33 : max_correction, correction_norm, max_corr_k, max_corr_j, ierr)
361 33 : if (ierr /= 0) then
362 0 : call oops('correction rejected by sizeB')
363 0 : exit iter_loop
364 : end if
365 :
366 33 : correction_norm = abs(correction_norm)
367 33 : max_abs_correction = abs(max_correction)
368 33 : corr_norm_min = min(correction_norm, corr_norm_min)
369 33 : max_corr_min = min(max_abs_correction, max_corr_min)
370 :
371 33 : if (is_bad_num(correction_norm) .or. is_bad_num(max_abs_correction)) then
372 : ! bad news -- bogus correction
373 0 : call oops('bad result from sizeB -- correction info either NaN or Inf')
374 0 : if (s% stop_for_bad_nums) then
375 0 : write(*,1) 'correction_norm', correction_norm
376 0 : write(*,1) 'max_correction', max_correction
377 0 : call mesa_error(__FILE__,__LINE__,'solver')
378 : end if
379 : exit iter_loop
380 : end if
381 :
382 33 : if (.not. s% ignore_too_large_correction) then
383 33 : if ((correction_norm > s% corr_param_factor*s% scale_correction_norm) .and. &
384 : .not. s% doing_first_model_of_run) then
385 0 : call oops('avg corr too large')
386 0 : exit iter_loop
387 : end if
388 : end if
389 :
390 : ! shrink the correction if it is too large
391 33 : correction_factor = 1d0
392 33 : temp_correction_factor = 1d0
393 :
394 33 : if (correction_norm*correction_factor > s% scale_correction_norm) then
395 0 : correction_factor = min(correction_factor,s% scale_correction_norm/correction_norm)
396 : end if
397 :
398 33 : if (max_abs_correction*correction_factor > s% scale_max_correction) then
399 0 : temp_correction_factor = s% scale_max_correction/max_abs_correction
400 : end if
401 :
402 33 : if (iter > s% solver_itermin_until_reduce_min_corr_coeff) then
403 0 : if (min_corr_coeff == 1d0 .and. &
404 : s% solver_reduced_min_corr_coeff < 1d0) then
405 0 : min_corr_coeff = s% solver_reduced_min_corr_coeff
406 : end if
407 : end if
408 :
409 33 : correction_factor = max(min_corr_coeff, correction_factor)
410 33 : if (.not. s% ignore_min_corr_coeff_for_scale_max_correction) then
411 33 : temp_correction_factor = max(min_corr_coeff, temp_correction_factor)
412 : end if
413 33 : correction_factor = min(correction_factor, temp_correction_factor)
414 :
415 : ! fix B if out of definition domain
416 33 : call Bdomain(s, nvar, soln, correction_factor, ierr)
417 33 : if (ierr /= 0) then ! correction cannot be fixed
418 0 : call oops('correction rejected by Bdomain')
419 0 : exit iter_loop
420 : end if
421 :
422 33 : if (min_corr_coeff < 1d0) then
423 : ! compute gradient of f = equ<dot>jacobian
424 : ! NOTE: NOT jacobian<dot>equ
425 0 : call block_multiply_xa(nvar, nz, lblk1, dblk1, ublk1, equ1, grad_f1)
426 :
427 0 : slope = eval_slope(nvar, nz, grad_f, soln)
428 0 : if (is_bad_num(slope) .or. slope > 0d0) then ! a very bad sign
429 0 : if (is_bad_num(slope) .and. s% stop_for_bad_nums) then
430 0 : write(*,1) 'slope', slope
431 0 : call mesa_error(__FILE__,__LINE__,'solver')
432 : end if
433 0 : slope = 0d0
434 0 : min_corr_coeff = 1d0
435 : end if
436 :
437 : else
438 :
439 33 : slope = 0d0
440 :
441 : end if
442 :
443 33 : f = 0d0
444 : call adjust_correction( &
445 33 : min_corr_coeff, correction_factor, grad_f1, f, slope, coeff, err_msg, ierr)
446 33 : if (ierr /= 0) then
447 0 : call oops(err_msg)
448 0 : exit iter_loop
449 : end if
450 33 : s% solver_adjust_iter = 0
451 :
452 : ! coeff is factor by which adjust_correction rescaled the correction vector
453 33 : if (coeff > s% tiny_corr_factor*min_corr_coeff .or. min_corr_coeff >= 1d0) then
454 33 : tiny_corr_cnt = 0
455 : else
456 0 : tiny_corr_cnt = tiny_corr_cnt + 1
457 : end if
458 :
459 : ! check the residuals for the equations
460 :
461 33 : call sizequ(s, nvar, residual_norm, max_residual, max_resid_k, max_resid_j, ierr)
462 33 : if (ierr /= 0) then
463 0 : call oops('sizequ returned ierr')
464 0 : exit iter_loop
465 : end if
466 :
467 33 : if (is_bad_num(residual_norm)) then
468 0 : call oops('residual_norm is a a bad number (NaN or Infinity)')
469 0 : if (s% stop_for_bad_nums) then
470 0 : write(*,1) 'residual_norm', residual_norm
471 0 : call mesa_error(__FILE__,__LINE__,'solver')
472 : end if
473 : exit iter_loop
474 : end if
475 :
476 33 : if (is_bad_num(max_residual)) then
477 0 : call oops('max_residual is a a bad number (NaN or Infinity)')
478 0 : if (s% stop_for_bad_nums) then
479 0 : write(*,1) 'max_residual', max_residual
480 0 : call mesa_error(__FILE__,__LINE__,'solver')
481 : end if
482 : exit iter_loop
483 : end if
484 :
485 33 : residual_norm = abs(residual_norm)
486 33 : max_residual = abs(max_residual)
487 33 : s% residual_norm = residual_norm
488 33 : s% max_residual = max_residual
489 33 : resid_norm_min = min(residual_norm, resid_norm_min)
490 33 : max_resid_min = min(max_residual, max_resid_min)
491 :
492 : disabled_resid_tests = &
493 33 : tol_max_residual > 1d2 .and. tol_residual_norm > 1d2
494 : pass_resid_tests = &
495 : .not. disabled_resid_tests .and. &
496 : max_residual <= tol_max_residual .and. &
497 33 : residual_norm <= tol_residual_norm
498 : pass_corr_tests_without_coeff = &
499 : max_abs_correction <= tol_max_correction .and. &
500 33 : correction_norm <= tol_correction_norm
501 : pass_corr_tests_with_coeff = &
502 : max_abs_correction <= tol_max_correction*coeff .and. &
503 33 : correction_norm <= tol_correction_norm*coeff
504 :
505 : passed_tol_tests = &
506 : (pass_resid_tests .and. pass_corr_tests_with_coeff) .or. &
507 33 : (disabled_resid_tests .and. pass_corr_tests_without_coeff)
508 :
509 33 : if (.not. passed_tol_tests) then
510 :
511 22 : if (iter >= max_tries) then
512 0 : call get_message
513 0 : message = trim(message) // ' -- give up'
514 0 : if (len_trim(s% retry_message) == 0) &
515 0 : s% retry_message = trim(message) // ' in solver'
516 0 : if (dbg_msg) call write_msg(message)
517 0 : if (.not. pass_resid_tests .and. .not. disabled_resid_tests) then
518 0 : if (residual_norm > tol_residual_norm) &
519 0 : write(*,2) 'residual_norm > tol_residual_norm', &
520 0 : s% model_number, residual_norm, tol_residual_norm
521 0 : if (max_residual > tol_max_residual) &
522 0 : write(*,2) 'max_residual > tol_max_residual', &
523 0 : s% model_number, max_residual, tol_max_residual
524 : end if
525 0 : if (disabled_resid_tests) then ! no coeff for corrections
526 0 : if (correction_norm > tol_correction_norm) &
527 0 : write(*,2) 'correction_norm > tol_correction_norm', &
528 0 : s% model_number, correction_norm, tol_correction_norm, coeff
529 0 : if (max_abs_correction > tol_max_correction) &
530 0 : write(*,2) 'max_abs_correction > tol_max_correction', &
531 0 : s% model_number, max_abs_correction, tol_max_correction, coeff
532 : else ! include coeff for corrections
533 0 : if (correction_norm > tol_correction_norm*coeff) &
534 0 : write(*,2) 'correction_norm > tol_correction_norm*coeff', &
535 0 : s% model_number, correction_norm, tol_correction_norm*coeff, coeff
536 0 : if (max_abs_correction > tol_max_correction*coeff) &
537 0 : write(*,2) 'max_abs_correction > tol_max_correction*coeff', &
538 0 : s% model_number, max_abs_correction, tol_max_correction*coeff, coeff
539 : end if
540 0 : convergence_failure = .true.
541 0 : exit iter_loop
542 22 : else if (.not. first_try .and. .not. s% doing_first_model_of_run) then
543 10 : if (correction_norm > s% corr_norm_jump_limit*corr_norm_min) then
544 0 : call oops('avg correction jumped')
545 0 : exit iter_loop
546 10 : else if (residual_norm > s% resid_norm_jump_limit*resid_norm_min) then
547 0 : call oops('avg residual jumped')
548 0 : exit iter_loop
549 10 : else if (max_abs_correction > s% max_corr_jump_limit*max_corr_min) then
550 0 : call oops('max correction jumped')
551 0 : exit iter_loop
552 10 : else if (max_residual > s% max_resid_jump_limit*max_resid_min) then
553 0 : call oops('max residual jumped')
554 0 : exit iter_loop
555 : else if (tiny_corr_cnt >= s% tiny_corr_coeff_limit &
556 10 : .and. min_corr_coeff < 1) then
557 0 : call oops('tiny corrections')
558 0 : exit iter_loop
559 : end if
560 12 : else if (.not. s% doing_first_model_of_run) then
561 10 : if (coeff < min(min_corr_coeff,correction_factor)) then
562 0 : call oops('coeff too small')
563 0 : exit iter_loop
564 : end if
565 : end if
566 : end if
567 :
568 33 : if (dbg_msg) then
569 0 : if (.not. passed_tol_tests) then
570 0 : call get_message
571 : end if
572 0 : if (.not. passed_tol_tests) then
573 0 : call write_msg(message)
574 0 : else if (iter < s% solver_itermin) then
575 0 : call write_msg('iter < itermin')
576 : else
577 0 : call write_msg('okay!')
578 : end if
579 : end if
580 :
581 33 : if (passed_tol_tests .and. (iter+1 < max_tries)) then
582 : ! about to declare victory... but may want to do another iteration
583 11 : force_iter_value = force_another_iteration(s, iter, s% solver_itermin)
584 11 : if (force_iter_value > 0) then
585 0 : passed_tol_tests = .false. ! force another
586 0 : tiny_corr_cnt = 0 ! reset the counter
587 0 : corr_norm_min = 1d99
588 0 : resid_norm_min = 1d99
589 0 : max_corr_min = 1d99
590 0 : max_resid_min = 1d99
591 11 : else if (force_iter_value < 0) then ! failure
592 0 : call oops('force iter')
593 0 : exit iter_loop
594 : end if
595 : end if
596 :
597 33 : if (s% use_other_solver_monitor .and. &
598 : associated(s% other_solver_monitor)) then
599 : call s% other_solver_monitor( &
600 : s% id, iter, passed_tol_tests, &
601 : correction_norm, max_correction, &
602 0 : residual_norm, max_residual, ierr)
603 0 : if (ierr /= 0) then
604 0 : call oops('other_solver_monitor')
605 0 : exit iter_loop
606 : end if
607 : end if
608 :
609 33 : iter=iter+1
610 33 : s% solver_iter = iter
611 33 : first_try = .false.
612 :
613 : end do iter_loop
614 :
615 11 : if (max_residual > s% warning_limit_for_max_residual .and. .not. convergence_failure) &
616 0 : write(*,2) 'WARNING: max_residual > warning_limit_for_max_residual', &
617 0 : s% model_number, max_residual, s% warning_limit_for_max_residual
618 :
619 : end subroutine do_solver_work
620 :
621 :
622 33 : subroutine solver_test_partials(nvar, xder, ierr)
623 : ! create jacobian by using numerical differences for partial derivatives
624 : integer, intent(in) :: nvar
625 : real(dp), pointer, dimension(:,:) :: xder ! (nvar, nz)
626 : integer, intent(out) :: ierr
627 :
628 : integer :: j, k, k_lo, k_hi
629 33 : real(dp), dimension(:,:), pointer :: save_equ, save_dx
630 : logical :: testing_partial
631 :
632 : include 'formats'
633 :
634 33 : ierr = 0
635 : testing_partial = & ! check inlist parameters
636 : s% solver_test_partials_dx_0 > 0d0 .and. &
637 : s% solver_test_partials_k > 0 .and. &
638 : s% solver_call_number == s% solver_test_partials_call_number .and. &
639 0 : s% solver_test_partials_iter_number == iter
640 33 : if (.not. testing_partial) return
641 :
642 0 : call do_equations(ierr)
643 0 : if (ierr /= 0) return
644 :
645 0 : allocate(save_dx(nvar,nz), save_equ(nvar,nz))
646 :
647 0 : do k=1,nz
648 0 : do j=1,nvar
649 0 : save_dx(j,k) = s% solver_dx(j,k)
650 0 : save_equ(j,k) = equ(j,k)
651 : end do
652 : end do
653 :
654 0 : s% doing_check_partials = .true. ! let set_vars_for_solver know
655 0 : k_lo = s% solver_test_partials_k_low
656 0 : if (k_lo > 0 .and. k_lo <= s% nz) then
657 0 : k_hi = s% solver_test_partials_k_high
658 0 : if (k_hi <= 0) then
659 : k_hi = s% nz
660 : else
661 0 : k_hi = min(k_hi,s% nz)
662 : end if
663 0 : do k = k_lo, k_hi
664 0 : call test_cell_partials(s, k, save_dx, save_equ, ierr)
665 0 : if (ierr /= 0) call mesa_error(__FILE__,__LINE__,'failed solver_test_partials')
666 : end do
667 : else
668 0 : k = s% solver_test_partials_k
669 0 : call test_cell_partials(s, k, save_dx, save_equ, ierr)
670 0 : if (ierr /= 0) call mesa_error(__FILE__,__LINE__,'failed solver_test_partials')
671 : end if
672 0 : deallocate(save_dx, save_equ)
673 0 : call mesa_error(__FILE__,__LINE__,'done solver_test_partials')
674 :
675 33 : end subroutine solver_test_partials
676 :
677 :
678 0 : subroutine get_message
679 : include 'formats'
680 0 : i = 0
681 0 : if (correction_norm > tol_correction_norm*coeff) i = i+1
682 0 : if (max_abs_correction > tol_max_correction*coeff) i = i+2
683 0 : if (residual_norm > tol_residual_norm*coeff) i = i+4
684 0 : if (max_residual > tol_max_residual*coeff) i = i+8
685 0 : if (i == 0) then
686 0 : message = 'out of tries'
687 : else
688 0 : message = tol_msg(i)
689 : end if
690 0 : end subroutine get_message
691 :
692 :
693 11 : subroutine set_param_defaults
694 11 : if (s% corr_param_factor == 0) s% corr_param_factor = 10d0
695 11 : if (s% scale_max_correction == 0) s% scale_max_correction = 1d99
696 11 : if (s% corr_norm_jump_limit == 0) s% corr_norm_jump_limit = 1d99
697 11 : if (s% max_corr_jump_limit == 0) s% max_corr_jump_limit = 1d99
698 11 : end subroutine set_param_defaults
699 :
700 :
701 0 : subroutine oops(msg)
702 : character (len=*), intent(in) :: msg
703 : character (len=strlen) :: full_msg
704 : include 'formats'
705 0 : full_msg = trim(msg) // ' -- give up'
706 0 : if (len_trim(s% retry_message) == 0) s% retry_message = trim(full_msg) // ' in solver'
707 0 : call write_msg(full_msg)
708 0 : convergence_failure = .true.
709 0 : end subroutine oops
710 :
711 :
712 132 : subroutine do_equations(ierr)
713 : integer, intent(out) :: ierr
714 44 : call prepare_solver_matrix(nvar, xder, ierr)
715 44 : if (ierr /= 0) return
716 44 : call eval_equations(s, nvar, ierr)
717 44 : if (ierr /= 0) return
718 44 : call s% other_after_solver_setmatrix(s% id, ierr)
719 : end subroutine do_equations
720 :
721 :
722 44 : subroutine prepare_solver_matrix(nvar, xder, ierr)
723 : integer, intent(in) :: nvar
724 : real(dp), pointer, dimension(:,:) :: xder ! (nvar, nz)
725 : integer, intent(out) :: ierr
726 : integer :: nz, neqns
727 : include 'formats'
728 44 : ierr = 0
729 44 : nz = s% nz
730 44 : neqns = nvar*nz
731 44 : s% ublk(1:nvar,1:nvar,1:nz) => ublk1
732 44 : s% dblk(1:nvar,1:nvar,1:nz) => dblk1
733 44 : s% lblk(1:nvar,1:nvar,1:nz) => lblk1
734 44 : end subroutine prepare_solver_matrix
735 :
736 :
737 33 : subroutine adjust_correction( &
738 : min_corr_coeff_in, max_corr_coeff, grad_f, f, slope, coeff, &
739 : err_msg, ierr)
740 : real(dp), intent(in) :: min_corr_coeff_in
741 : real(dp), intent(in) :: max_corr_coeff
742 : real(dp), intent(in) :: grad_f(:) ! (neq) ! gradient df/ddx at xold
743 : real(dp), intent(out) :: f ! 1/2 fvec^2. minimize this.
744 : real(dp), intent(in) :: slope
745 : real(dp), intent(out) :: coeff
746 :
747 : ! the new correction is coeff*xscale*soln
748 : ! with min_corr_coeff <= coeff <= max_corr_coeff
749 : ! if all goes well, the new x will give an improvement in f
750 :
751 : character (len=*), intent(out) :: err_msg
752 : integer, intent(out) :: ierr
753 :
754 : integer :: i, k, iter
755 : logical :: first_time
756 33 : real(dp) :: a1, alam, alam2, a2, disc, f2, &
757 33 : rhs1, rhs2, tmplam, fold, min_corr_coeff
758 33 : real(dp) :: f_target
759 : logical :: skip_eval_f, dbg_adjust
760 :
761 : real(dp), parameter :: alf = 1d-2 ! ensures sufficient decrease in f
762 :
763 : real(dp), parameter :: alam_factor = 0.2d0
764 :
765 : include 'formats'
766 :
767 33 : ierr = 0
768 33 : coeff = 0
769 33 : dbg_adjust = .false.
770 33 : err_msg = ''
771 :
772 33 : skip_eval_f = (min_corr_coeff_in == 1)
773 33 : if (skip_eval_f) then
774 33 : f = 0
775 : else
776 0 : do k=1,nz
777 0 : do i=1,nvar
778 0 : dxsave(i,k) = s% solver_dx(i,k)
779 0 : ddxsave(i,k) = ddx(i,k)
780 : end do
781 : end do
782 0 : f = eval_f(nvar,nz,equ)
783 0 : if (is_bad_num(f)) then
784 0 : ierr = -1
785 0 : write(err_msg,*) 'adjust_correction failed in eval_f'
786 0 : if (dbg_msg) write(*,*) &
787 0 : 'adjust_correction: eval_f(nvar,nz,equ)', eval_f(nvar,nz,equ)
788 0 : if (s% stop_for_bad_nums) then
789 0 : write(*,1) 'f', f
790 0 : call mesa_error(__FILE__,__LINE__,'solver adjust_correction')
791 : end if
792 33 : return
793 : end if
794 : end if
795 33 : fold = f
796 33 : min_corr_coeff = min(min_corr_coeff_in, max_corr_coeff) ! make sure min <= max
797 33 : alam = max_corr_coeff
798 33 : first_time = .true.
799 33 : f2 = 0
800 33 : alam2 = 0
801 : if (dbg_adjust) then
802 : write(*,4) 'max_corr_coeff', k, s% solver_iter, &
803 : s% model_number, max_corr_coeff
804 : write(*,4) 'slope', k, s% solver_iter, &
805 : s% model_number, slope
806 : write(*,4) 'f', k, s% solver_iter, &
807 : s% model_number, f
808 : end if
809 :
810 33 : search_loop: do iter = 1, 1000
811 :
812 33 : coeff = max(min_corr_coeff, alam)
813 33 : s% solver_adjust_iter = iter
814 :
815 33 : call apply_coeff(nvar, nz, dxsave, soln, coeff, skip_eval_f)
816 :
817 33 : call do_equations(ierr)
818 33 : if (ierr /= 0) then
819 0 : if (alam > min_corr_coeff .and. s% model_number == 1) then
820 : ! try again with smaller correction vector.
821 : ! need this to rescue create pre-main-sequence model in some nasty cases.
822 0 : alam = max(alam/10, min_corr_coeff)
823 0 : ierr = 0
824 0 : cycle search_loop
825 : end if
826 0 : write(err_msg,*) 'adjust_correction failed in eval_equations'
827 0 : if (dbg_msg .or. dbg_adjust) &
828 0 : write(*,2) 'adjust_correction: eval_equations returned ierr', &
829 0 : ierr, min_corr_coeff, max_corr_coeff
830 : exit search_loop
831 : end if
832 :
833 33 : if (min_corr_coeff == 1) return
834 :
835 : if (dbg_adjust) then
836 : do k=1,nz
837 : do i=1,nvar
838 : write(*,5) trim(s% nameofequ(i)), k, iter, s% solver_iter, &
839 : s% model_number, equ(i,k)
840 : end do
841 : end do
842 : end if
843 :
844 0 : f = eval_f(nvar,nz,equ)
845 0 : if (is_bad_num(f)) then
846 0 : if (s% stop_for_bad_nums) then
847 0 : write(*,1) 'f', f
848 0 : call mesa_error(__FILE__,__LINE__,'solver adjust_correction eval_f')
849 : end if
850 0 : if (alam > min_corr_coeff) then
851 0 : alam = max(alam/10, min_corr_coeff)
852 0 : ierr = 0
853 0 : cycle search_loop
854 : end if
855 0 : err_msg = 'equ norm is NaN or other bad num'
856 0 : ierr = -1
857 0 : exit search_loop
858 : end if
859 :
860 0 : f_target = max(fold/2, fold + alf*coeff*slope)
861 0 : if (f <= f_target) then
862 : return ! sufficient decrease in f
863 : end if
864 :
865 0 : if (alam <= min_corr_coeff) then
866 : return ! time to give up
867 : end if
868 :
869 : ! reduce alam and try again
870 0 : if (first_time) then
871 0 : tmplam = -slope/(2*(f-fold-slope))
872 0 : first_time = .false.
873 : if (dbg_adjust) then
874 : write(*,5) 'slope', k, iter, s% solver_iter, &
875 : s% model_number, slope
876 : write(*,5) 'f', k, iter, s% solver_iter, &
877 : s% model_number, f
878 : write(*,5) 'fold', k, iter, s% solver_iter, &
879 : s% model_number, fold
880 : write(*,5) '2*(f-fold-slope)', k, iter, s% solver_iter, &
881 : s% model_number, 2*(f-fold-slope)
882 : end if
883 : else ! have two prior f values to work with
884 0 : rhs1 = f - fold - alam*slope
885 0 : rhs2 = f2 - fold - alam2*slope
886 0 : a1 = (rhs1/(alam*alam) - rhs2/(alam2*alam2))/(alam - alam2)
887 0 : a2 = (-alam2*rhs1/(alam*alam) + alam*rhs2/(alam2*alam2))/(alam - alam2)
888 : if (dbg_adjust) then
889 : write(*,5) 'slope', k, iter, s% solver_iter, &
890 : s% model_number, slope
891 : write(*,5) 'f', k, iter, s% solver_iter, &
892 : s% model_number, f
893 : write(*,5) 'f2', k, iter, s% solver_iter, &
894 : s% model_number, f2
895 : write(*,5) 'fold', k, iter, s% solver_iter, &
896 : s% model_number, fold
897 : write(*,5) 'alam', k, iter, s% solver_iter, &
898 : s% model_number, alam
899 : write(*,5) 'alam2', k, iter, s% solver_iter, &
900 : s% model_number, alam2
901 : write(*,5) 'rhs1', k, iter, s% solver_iter, &
902 : s% model_number, rhs1
903 : write(*,5) 'rhs2', k, iter, s% solver_iter, &
904 : s% model_number, rhs2
905 : write(*,5) 'a1', k, iter, s% solver_iter, &
906 : s% model_number, a1
907 : write(*,5) 'a2', k, iter, s% solver_iter, &
908 : s% model_number, a2
909 : end if
910 0 : if (a1 == 0) then
911 0 : tmplam = -slope/(2*a2)
912 : else
913 0 : disc = a2*a2-3*a1*slope
914 0 : if (disc < 0) then
915 0 : tmplam = alam*alam_factor
916 0 : else if (a2 <= 0) then
917 0 : tmplam = (-a2+sqrt(disc))/(3*a1)
918 : else
919 0 : tmplam = -slope/(a2+sqrt(disc))
920 : end if
921 : if (dbg_adjust) then
922 : write(*,5) 'disc', k, iter, s% solver_iter, &
923 : s% model_number, disc
924 : end if
925 : end if
926 0 : if (tmplam > alam*alam_factor) tmplam = alam*alam_factor
927 : end if
928 :
929 0 : alam2 = alam
930 0 : f2 = f
931 0 : alam = max(tmplam, alam*alam_factor, min_corr_coeff)
932 :
933 33 : if (dbg_adjust) then
934 : write(*,5) 'tmplam', k, iter, s% solver_iter, &
935 : s% model_number, tmplam
936 : write(*,5) 'min_corr_coeff', k, iter, s% solver_iter, &
937 : s% model_number, min_corr_coeff
938 : write(*,5) 'alam_factor', k, iter, s% solver_iter, &
939 : s% model_number, alam_factor
940 : end if
941 :
942 : end do search_loop
943 :
944 0 : do k=1,nz
945 0 : do i=1,nvar
946 0 : s% solver_dx(i,k) = dxsave(i,k)
947 0 : ddx(i,k) = ddxsave(i,k)
948 : end do
949 : end do
950 :
951 33 : end subroutine adjust_correction
952 :
953 :
954 33 : subroutine apply_coeff(nvar, nz, dxsave, soln, coeff, just_use_dx)
955 : integer, intent(in) :: nvar, nz
956 : real(dp), intent(in), dimension(:,:) :: dxsave, soln
957 : real(dp), intent(in) :: coeff
958 : logical, intent(in) :: just_use_dx
959 : integer :: i, k
960 : include 'formats'
961 :
962 33 : if (just_use_dx) then
963 33 : if (coeff == 1d0) then
964 39294 : do k=1,nz
965 510426 : do i=1,nvar
966 510393 : s% solver_dx(i,k) = s% solver_dx(i,k) + s% x_scale(i,k)*soln(i,k)
967 : end do
968 : end do
969 : else
970 0 : do k=1,nz
971 0 : do i=1,nvar
972 0 : s% solver_dx(i,k) = s% solver_dx(i,k) + coeff*s% x_scale(i,k)*soln(i,k)
973 : end do
974 : end do
975 : end if
976 33 : return
977 : end if
978 : ! else use dxsave instead of dx
979 0 : if (coeff == 1d0) then
980 0 : do k=1,nz
981 0 : do i=1,nvar
982 0 : s% solver_dx(i,k) = dxsave(i,k) + s% x_scale(i,k)*soln(i,k)
983 : end do
984 : end do
985 : return
986 : end if
987 0 : do k=1,nz
988 0 : do i=1,nvar
989 0 : s% solver_dx(i,k) = dxsave(i,k) + coeff*s% x_scale(i,k)*soln(i,k)
990 : end do
991 : end do
992 : end subroutine apply_coeff
993 :
994 :
995 33 : logical function solve_equ()
996 : use star_utils, only: start_time, update_time
997 : use rsp_def, only: NV, MAX_NZN
998 : integer :: i, ierr
999 33 : real(dp) :: total_time
1000 :
1001 : include 'formats'
1002 33 : ierr = 0
1003 33 : solve_equ = .true.
1004 :
1005 33 : if (s% doing_timing) then
1006 0 : call start_time(s, time0, total_time)
1007 : end if
1008 :
1009 33 : !$OMP PARALLEL DO SIMD
1010 : do i=1,neq
1011 : b1(i) = -equ1(i)
1012 : end do
1013 : !$OMP END PARALLEL DO SIMD
1014 :
1015 33 : if (s% use_DGESVX_in_bcyclic) then
1016 0 : !$OMP PARALLEL DO SIMD
1017 : do i = 1, nvar*nvar*nz
1018 : save_ublk1(i) = ublk1(i)
1019 : save_dblk1(i) = dblk1(i)
1020 : save_lblk1(i) = lblk1(i)
1021 : end do
1022 : !$OMP END PARALLEL DO SIMD
1023 : end if
1024 :
1025 33 : call factor_mtx(ierr)
1026 33 : if (ierr == 0) call solve_mtx(ierr)
1027 :
1028 33 : if (s% use_DGESVX_in_bcyclic) then
1029 0 : !$OMP PARALLEL DO SIMD
1030 : do i = 1, nvar*nvar*nz
1031 : ublk1(i) = save_ublk1(i)
1032 : dblk1(i) = save_dblk1(i)
1033 : lblk1(i) = save_lblk1(i)
1034 : end do
1035 : !$OMP END PARALLEL DO SIMD
1036 : end if
1037 :
1038 33 : if (s% doing_timing) then
1039 0 : call update_time(s, time0, total_time, s% time_solver_matrix)
1040 : end if
1041 33 : if (ierr /= 0) then
1042 0 : solve_equ = .false.
1043 0 : b(1:nvar,1:nz) = 0d0
1044 : end if
1045 :
1046 33 : end function solve_equ
1047 :
1048 :
1049 33 : subroutine factor_mtx(ierr)
1050 33 : use star_bcyclic, only: bcyclic_factor
1051 : integer, intent(out) :: ierr
1052 : call bcyclic_factor( &
1053 : s, nvar, nz, lblk1, dblk1, ublk1, lblkF1, dblkF1, ublkF1, ipiv_blk1, &
1054 : B1, row_scale_factors1, col_scale_factors1, &
1055 33 : equed1, iter, ierr)
1056 33 : end subroutine factor_mtx
1057 :
1058 :
1059 33 : subroutine solve_mtx(ierr)
1060 33 : use star_bcyclic, only: bcyclic_solve
1061 : integer, intent(out) :: ierr
1062 : call bcyclic_solve( &
1063 : s, nvar, nz, lblk1, dblk1, ublk1, lblkF1, dblkF1, ublkF1, ipiv_blk1, &
1064 : B1, soln1, row_scale_factors1, col_scale_factors1, equed1, &
1065 33 : iter, ierr)
1066 33 : end subroutine solve_mtx
1067 :
1068 :
1069 0 : subroutine test_cell_partials(s, k, save_dx, save_equ, ierr)
1070 33 : use star_utils, only: lookup_nameofvar, lookup_nameofequ
1071 : type (star_info), pointer :: s
1072 : integer, intent(in) :: k
1073 : real(dp), pointer, dimension(:,:) :: save_dx, save_equ
1074 : integer, intent(out) :: ierr
1075 : integer :: i_equ, i_var, i_var_sink, i_var_xa_index, i_var_sink_xa_index
1076 : include 'formats'
1077 0 : ierr = 0
1078 0 : write(*,'(A)')
1079 0 : i_equ = lookup_nameofequ(s, s% solver_test_partials_equ_name)
1080 0 : if (i_equ == 0 .and. len_trim(s% solver_test_partials_equ_name) > 0) then
1081 0 : if (s% solver_test_partials_equ_name == 'lnE') then ! testing eos
1082 0 : i_equ = -1
1083 0 : else if (s% solver_test_partials_equ_name == 'eps_nuc') then
1084 0 : i_equ = -2
1085 0 : else if (s% solver_test_partials_equ_name == 'opacity') then
1086 0 : i_equ = -3
1087 0 : else if (s% solver_test_partials_equ_name == 'lnP') then
1088 0 : i_equ = -4
1089 0 : else if (s% solver_test_partials_equ_name == 'non_nuc_neu') then
1090 0 : i_equ = -5
1091 0 : else if (s% solver_test_partials_equ_name == 'gradT') then
1092 0 : i_equ = -6
1093 0 : else if (s% solver_test_partials_equ_name == 'mlt_vc') then
1094 0 : i_equ = -7
1095 0 : else if (s% solver_test_partials_equ_name == 'grad_ad') then
1096 0 : i_equ = -8
1097 : end if
1098 0 : else if (i_equ /= 0) then
1099 0 : write(*,1) 'equ name ' // trim(s% solver_test_partials_equ_name)
1100 : end if
1101 0 : i_var = lookup_nameofvar(s, s% solver_test_partials_var_name)
1102 0 : if (i_var /= 0) write(*,1) 'var name ' // trim(s% solver_test_partials_var_name)
1103 0 : if (i_var > s% nvar_hydro) then ! get index in xa
1104 0 : i_var_xa_index = i_var - s% nvar_hydro
1105 : else
1106 0 : i_var_xa_index = 0
1107 : end if
1108 0 : i_var_sink = lookup_nameofvar(s, s% solver_test_partials_sink_name)
1109 0 : i_var_sink_xa_index = 0
1110 0 : if (i_var_sink > 0 .and. i_var > s% nvar_hydro) then
1111 0 : write(*,1) 'sink name ' // trim(s% solver_test_partials_sink_name)
1112 0 : if (i_var_sink > s% nvar_hydro) then ! get index in xa
1113 0 : i_var_sink_xa_index = i_var_sink - s% nvar_hydro
1114 : else
1115 0 : write(*,*) 'ERROR: sink name must be a chem name for the current net'
1116 0 : ierr = -1
1117 : return
1118 : end if
1119 : end if
1120 0 : if (s% solver_test_partials_equ_name == 'all') then
1121 0 : do i_equ = 1, s% nvar_hydro
1122 : call test_equ_partials(s, &
1123 : i_equ, i_var, i_var_sink, i_var_xa_index, i_var_sink_xa_index, &
1124 0 : k, save_dx, save_equ, ierr)
1125 0 : if (ierr /= 0) call mesa_error(__FILE__,__LINE__,'failed solver_test_partials')
1126 : end do
1127 : else
1128 : call test_equ_partials(s, &
1129 : i_equ, i_var, i_var_sink, i_var_xa_index, i_var_sink_xa_index, &
1130 0 : k, save_dx, save_equ, ierr)
1131 0 : if (ierr /= 0) call mesa_error(__FILE__,__LINE__,'failed solver_test_partials')
1132 : end if
1133 0 : end subroutine test_cell_partials
1134 :
1135 :
1136 0 : subroutine test_equ_partials(s, &
1137 : i_equ, i_var, i_var_sink, i_var_xa_index, i_var_sink_xa_index, &
1138 : k, save_dx, save_equ, ierr)
1139 : type (star_info), pointer :: s
1140 : integer, intent(in) :: &
1141 : i_equ, i_var, i_var_sink, i_var_xa_index, i_var_sink_xa_index, k
1142 : real(dp), pointer, dimension(:,:) :: save_dx, save_equ
1143 : integer, intent(out) :: ierr
1144 : integer :: i, j_var_xa_index, j_var_sink_xa_index
1145 : include 'formats'
1146 0 : if (i_equ /= 0) then
1147 0 : if (s% solver_test_partials_var_name == 'all') then
1148 0 : do i = 1, s% nvar_hydro
1149 : call test3_partials(s, &
1150 : i_equ, i, i_var_sink, i_var_xa_index, i_var_sink_xa_index, &
1151 0 : k, save_dx, save_equ, ierr)
1152 0 : if (ierr /= 0) call mesa_error(__FILE__,__LINE__,'failed solver_test_partials')
1153 0 : write(*,'(A)')
1154 : end do
1155 0 : else if (i_var == 0) then
1156 0 : write(*,*) 'failed to recognize variable name'
1157 : else
1158 : call test3_partials(s, &
1159 : i_equ, i_var, i_var_sink, i_var_xa_index, i_var_sink_xa_index, &
1160 0 : k, save_dx, save_equ, ierr)
1161 0 : if (ierr /= 0) call mesa_error(__FILE__,__LINE__,'failed solver_test_partials')
1162 : end if
1163 : else ! i_equ == 0
1164 0 : if (i_var /= 0) then
1165 : call test1_partial(s, &
1166 : i_equ, i_var, i_var_sink, i_var_xa_index, i_var_sink_xa_index, &
1167 0 : k, 0, s% solver_test_partials_dval_dx, save_dx, save_equ, ierr)
1168 : else ! i_var == 0
1169 0 : if (s% solver_test_partials_var <= 0) then
1170 0 : write(*,2) 'need to set solver_test_partials_var', s% solver_test_partials_var
1171 0 : write(*,2) 'for solver_test_partials_k', s% solver_test_partials_k
1172 0 : call mesa_error(__FILE__,__LINE__,'failed solver_test_partials')
1173 : end if
1174 0 : if (s% solver_test_partials_var > s% nvar_hydro) then
1175 0 : j_var_xa_index = s% solver_test_partials_var - s% nvar_hydro
1176 0 : if (s% solver_test_partials_dx_sink > s% nvar_hydro) then
1177 0 : j_var_sink_xa_index = s% solver_test_partials_dx_sink - s% nvar_hydro
1178 : else
1179 0 : write(*,*) 'set solver_test_partials_dx_sink to variable index, not to xa index', &
1180 0 : s% solver_test_partials_dx_sink
1181 0 : call mesa_error(__FILE__,__LINE__,'failed solver_test_partials')
1182 : end if
1183 : else
1184 0 : j_var_xa_index = 0
1185 0 : j_var_sink_xa_index = 0
1186 : end if
1187 : call test1_partial(s, &
1188 : i_equ, s% solver_test_partials_var, s% solver_test_partials_dx_sink, &
1189 : j_var_xa_index, j_var_sink_xa_index, &
1190 0 : k, 0, s% solver_test_partials_dval_dx, save_dx, save_equ, ierr)
1191 : end if
1192 0 : if (ierr /= 0) call mesa_error(__FILE__,__LINE__,'failed solver_test_partials')
1193 : end if
1194 0 : write(*,'(A)')
1195 0 : end subroutine test_equ_partials
1196 :
1197 :
1198 0 : subroutine get_lnE_partials(s, &
1199 : k, i_var, i_var_sink, i_var_xa_index, i_var_sink_xa_index, &
1200 : dvardx0_m1, dvardx0_00, dvardx0_p1)
1201 : use eos_def, only: i_lnE
1202 : type (star_info), pointer :: s
1203 : integer, intent(in) :: k, i_var, i_var_sink, i_var_xa_index, i_var_sink_xa_index
1204 : real(dp), intent(out) :: dvardx0_m1, dvardx0_00, dvardx0_p1
1205 0 : dvardx0_m1 = 0d0; dvardx0_00 = 0d0; dvardx0_p1 = 0d0
1206 0 : if (i_var_xa_index > 0) then
1207 : dvardx0_00 = s% d_eos_dxa(i_lnE,i_var_xa_index,k) - &
1208 0 : s% d_eos_dxa(i_lnE,i_var_sink_xa_index,k)
1209 0 : else if (i_var == s% i_lnd) then
1210 0 : dvardx0_00 = s% dE_dRho_for_partials(k)*s% rho(k)/s% energy(k)
1211 0 : else if (i_var == s% i_lnT) then
1212 0 : dvardx0_00 = s% Cv_for_partials(k)*s% T(k)/s% energy(k)
1213 : end if
1214 0 : end subroutine get_lnE_partials
1215 :
1216 :
1217 0 : subroutine get_lnP_partials(s, &
1218 : k, i_var, i_var_sink, i_var_xa_index, i_var_sink_xa_index, &
1219 : dvardx0_m1, dvardx0_00, dvardx0_p1)
1220 0 : use eos_def, only: i_lnPgas
1221 : type (star_info), pointer :: s
1222 : integer, intent(in) :: k, i_var, i_var_sink, i_var_xa_index, i_var_sink_xa_index
1223 : real(dp), intent(out) :: dvardx0_m1, dvardx0_00, dvardx0_p1
1224 0 : dvardx0_m1 = 0d0; dvardx0_00 = 0d0; dvardx0_p1 = 0d0
1225 0 : if (i_var_xa_index > 0) then
1226 : dvardx0_00 = s% Pgas(k)/s% Peos(k) * &
1227 0 : (s% d_eos_dxa(i_lnPgas,i_var_xa_index,k) - s% d_eos_dxa(i_lnPgas,i_var_sink_xa_index,k))
1228 0 : else if (i_var == s% i_lnd) then
1229 0 : dvardx0_00 = s% chiRho_for_partials(k)
1230 0 : else if (i_var == s% i_lnT) then
1231 0 : dvardx0_00 = s% chiT_for_partials(k)
1232 : end if
1233 0 : end subroutine get_lnP_partials
1234 :
1235 :
1236 0 : subroutine get_grad_ad_partials(s, &
1237 : k, i_var, i_var_sink, i_var_xa_index, i_var_sink_xa_index, &
1238 : dvardx0_m1, dvardx0_00, dvardx0_p1)
1239 0 : use eos_def, only: i_grad_ad
1240 : type (star_info), pointer :: s
1241 : integer, intent(in) :: k, i_var, i_var_sink, i_var_xa_index, i_var_sink_xa_index
1242 : real(dp), intent(out) :: dvardx0_m1, dvardx0_00, dvardx0_p1
1243 0 : dvardx0_m1 = 0d0; dvardx0_00 = 0d0; dvardx0_p1 = 0d0
1244 0 : if (i_var_xa_index > 0) then
1245 : dvardx0_00 = &
1246 0 : (s% d_eos_dxa(i_grad_ad,i_var_xa_index,k) - s% d_eos_dxa(i_grad_ad,i_var_sink_xa_index,k))
1247 0 : else if (i_var == s% i_lnd) then
1248 0 : dvardx0_00 = s% d_eos_dlnd(i_grad_ad,k)
1249 0 : else if (i_var == s% i_lnT) then
1250 0 : dvardx0_00 = s% d_eos_dlnT(i_grad_ad,k)
1251 : end if
1252 0 : end subroutine get_grad_ad_partials
1253 :
1254 :
1255 0 : subroutine get_eps_nuc_partials(s, &
1256 : k, i_var, i_var_sink, i_var_xa_index, i_var_sink_xa_index, &
1257 : dvardx0_m1, dvardx0_00, dvardx0_p1)
1258 : type (star_info), pointer :: s
1259 : integer, intent(in) :: k, i_var, i_var_sink, i_var_xa_index, i_var_sink_xa_index
1260 : real(dp), intent(out) :: dvardx0_m1, dvardx0_00, dvardx0_p1
1261 0 : dvardx0_m1 = 0d0; dvardx0_00 = 0d0; dvardx0_p1 = 0d0
1262 0 : if (i_var > s% nvar_hydro) then
1263 0 : dvardx0_00 = s% d_epsnuc_dx(i_var_xa_index,k) - s% d_epsnuc_dx(i_var_sink_xa_index,k)
1264 0 : else if (i_var == s% i_lnd) then
1265 0 : dvardx0_00 = s% d_epsnuc_dlnd(k)
1266 0 : else if (i_var == s% i_lnT) then
1267 0 : dvardx0_00 = s% d_epsnuc_dlnT(k)
1268 : end if
1269 0 : end subroutine get_eps_nuc_partials
1270 :
1271 :
1272 0 : subroutine get_non_nuc_neu_partials(s, &
1273 : k, i_var, i_var_sink, i_var_xa_index, i_var_sink_xa_index, &
1274 : dvardx0_m1, dvardx0_00, dvardx0_p1)
1275 : type (star_info), pointer :: s
1276 : integer, intent(in) :: k, i_var, i_var_sink, i_var_xa_index, i_var_sink_xa_index
1277 : real(dp), intent(out) :: dvardx0_m1, dvardx0_00, dvardx0_p1
1278 0 : dvardx0_m1 = 0d0; dvardx0_00 = 0d0; dvardx0_p1 = 0d0
1279 0 : if (i_var > s% nvar_hydro) then
1280 : dvardx0_00 = 0d0
1281 0 : else if (i_var == s% i_lnd) then
1282 0 : dvardx0_00 = s% d_nonnucneu_dlnd(k)
1283 0 : else if (i_var == s% i_lnT) then
1284 0 : dvardx0_00 = s% d_nonnucneu_dlnT(k)
1285 : end if
1286 0 : end subroutine get_non_nuc_neu_partials
1287 :
1288 :
1289 0 : subroutine get_gradT_partials(s, &
1290 : k, i_var, i_var_sink, i_var_xa_index, i_var_sink_xa_index, &
1291 : dvardx0_m1, dvardx0_00, dvardx0_p1)
1292 : type (star_info), pointer :: s
1293 : integer, intent(in) :: k, i_var, i_var_sink, i_var_xa_index, i_var_sink_xa_index
1294 : real(dp), intent(out) :: dvardx0_m1, dvardx0_00, dvardx0_p1
1295 0 : dvardx0_m1 = 0d0; dvardx0_00 = 0d0; dvardx0_p1 = 0d0
1296 0 : if (i_var > s% nvar_hydro) then
1297 : dvardx0_00 = 0d0
1298 0 : else if (i_var == s% i_lnd) then
1299 0 : dvardx0_m1 = s% gradT_ad(k)%d1Array(i_lnd_m1)
1300 0 : dvardx0_00 = s% gradT_ad(k)%d1Array(i_lnd_00)
1301 0 : else if (i_var == s% i_lnT) then
1302 0 : dvardx0_m1 = s% gradT_ad(k)%d1Array(i_lnT_m1)
1303 0 : dvardx0_00 = s% gradT_ad(k)%d1Array(i_lnT_00)
1304 0 : else if (i_var == s% i_lnR) then
1305 0 : dvardx0_00 = s% gradT_ad(k)%d1Array(i_lnR_00)
1306 0 : else if (i_var == s% i_lum) then
1307 0 : dvardx0_00 = s% gradT_ad(k)%d1Array(i_L_00)
1308 0 : else if (i_var == s% i_w_div_wc) then
1309 0 : dvardx0_00 = s% gradT_ad(k)%d1Array(i_w_00)
1310 : end if
1311 0 : end subroutine get_gradT_partials
1312 :
1313 :
1314 0 : subroutine get_mlt_vc_partials(s, &
1315 : k, i_var, i_var_sink, i_var_xa_index, i_var_sink_xa_index, &
1316 : dvardx0_m1, dvardx0_00, dvardx0_p1)
1317 : type (star_info), pointer :: s
1318 : integer, intent(in) :: k, i_var, i_var_sink, i_var_xa_index, i_var_sink_xa_index
1319 : real(dp), intent(out) :: dvardx0_m1, dvardx0_00, dvardx0_p1
1320 0 : dvardx0_m1 = 0d0; dvardx0_00 = 0d0; dvardx0_p1 = 0d0
1321 0 : if (i_var > s% nvar_hydro) then
1322 : dvardx0_00 = 0d0
1323 0 : else if (i_var == s% i_lnd) then
1324 0 : dvardx0_m1 = s% mlt_vc_ad(k)%d1Array(i_lnd_m1)
1325 0 : dvardx0_00 = s% mlt_vc_ad(k)%d1Array(i_lnd_00)
1326 0 : else if (i_var == s% i_lnT) then
1327 0 : dvardx0_m1 = s% mlt_vc_ad(k)%d1Array(i_lnT_m1)
1328 0 : dvardx0_00 = s% mlt_vc_ad(k)%d1Array(i_lnT_00)
1329 0 : else if (i_var == s% i_lnR) then
1330 0 : dvardx0_00 = s% mlt_vc_ad(k)%d1Array(i_lnR_00)
1331 0 : else if (i_var == s% i_lum) then
1332 0 : dvardx0_00 = s% mlt_vc_ad(k)%d1Array(i_L_00)
1333 : end if
1334 0 : end subroutine get_mlt_vc_partials
1335 :
1336 :
1337 0 : subroutine get_opacity_partials(s, &
1338 : k, i_var, i_var_sink, i_var_xa_index, i_var_sink_xa_index, &
1339 : dvardx0_m1, dvardx0_00, dvardx0_p1)
1340 : type (star_info), pointer :: s
1341 : integer, intent(in) :: k, i_var, i_var_sink, i_var_xa_index, i_var_sink_xa_index
1342 : real(dp), intent(out) :: dvardx0_m1, dvardx0_00, dvardx0_p1
1343 0 : dvardx0_m1 = 0d0; dvardx0_00 = 0d0; dvardx0_p1 = 0d0
1344 0 : if (i_var > s% nvar_hydro) then
1345 : dvardx0_00 = 0d0 ! s% d_opacity_dx(i_var_xa_index,k) - s% d_opacity_dx(i_var_sink_xa_index,k)
1346 0 : else if (i_var == s% i_lnd) then
1347 0 : dvardx0_00 = s% d_opacity_dlnd(k)
1348 0 : else if (i_var == s% i_lnT) then
1349 0 : dvardx0_00 = s% d_opacity_dlnT(k)
1350 : end if
1351 0 : end subroutine get_opacity_partials
1352 :
1353 :
1354 0 : subroutine test3_partials(s, &
1355 : i_equ, i_var, i_var_sink, i_var_xa_index, i_var_sink_xa_index, &
1356 : k, save_dx, save_equ, ierr)
1357 : type (star_info), pointer :: s
1358 : integer, intent(in) :: &
1359 : i_equ, i_var, i_var_sink, i_var_xa_index, i_var_sink_xa_index, k
1360 : real(dp), pointer, dimension(:,:) :: save_dx, save_equ
1361 : integer, intent(out) :: ierr
1362 : real(dp) :: dvardx0_m1, dvardx0_00, dvardx0_p1
1363 0 : dvardx0_m1 = 0d0
1364 0 : dvardx0_00 = 0d0
1365 0 : dvardx0_p1 = 0d0
1366 0 : if (i_equ > 0) then
1367 0 : if (i_var > s% nvar_hydro) then ! testing abundance
1368 0 : if (k > 1) dvardx0_m1 = &
1369 0 : s% lblk(i_equ,i_var,k)/s% x_scale(i_var,k-1) - s% lblk(i_equ,i_var_sink,k)/s% x_scale(i_var_sink,k-1)
1370 : dvardx0_00 = &
1371 0 : s% dblk(i_equ,i_var,k)/s% x_scale(i_var,k) - s% dblk(i_equ,i_var_sink,k)/s% x_scale(i_var_sink,k)
1372 0 : if (k < s% nz) dvardx0_p1 = &
1373 0 : s% ublk(i_equ,i_var,k)/s% x_scale(i_var,k+1) - s% ublk(i_equ,i_var_sink,k)/s% x_scale(i_var_sink,k+1)
1374 : else
1375 0 : if (k > 1) dvardx0_m1 = s% lblk(i_equ,i_var,k)/s% x_scale(i_var,k-1)
1376 0 : dvardx0_00 = s% dblk(i_equ,i_var,k)/s% x_scale(i_var,k)
1377 0 : if (k < s% nz) dvardx0_p1 = s% ublk(i_equ,i_var,k)/s% x_scale(i_var,k+1)
1378 : end if
1379 0 : else if (i_equ == -1) then ! 'lnE'
1380 : call get_lnE_partials(s, &
1381 : k, i_var, i_var_sink, i_var_xa_index, i_var_sink_xa_index, &
1382 0 : dvardx0_m1, dvardx0_00, dvardx0_p1)
1383 0 : elseif (i_equ == -2) then ! 'eps_nuc'
1384 : call get_eps_nuc_partials(s, &
1385 : k, i_var, i_var_sink, i_var_xa_index, i_var_sink_xa_index, &
1386 0 : dvardx0_m1, dvardx0_00, dvardx0_p1)
1387 0 : else if (i_equ == -3) then ! 'opacity'
1388 : call get_opacity_partials(s, &
1389 : k, i_var, i_var_sink, i_var_xa_index, i_var_sink_xa_index, &
1390 0 : dvardx0_m1, dvardx0_00, dvardx0_p1)
1391 0 : else if (i_equ == -4) then ! 'lnP'
1392 : call get_lnP_partials(s, &
1393 : k, i_var, i_var_sink, i_var_xa_index, i_var_sink_xa_index, &
1394 0 : dvardx0_m1, dvardx0_00, dvardx0_p1)
1395 0 : else if (i_equ == -5) then ! 'non_nuc_neu'
1396 : call get_non_nuc_neu_partials(s, &
1397 : k, i_var, i_var_sink, i_var_xa_index, i_var_sink_xa_index, &
1398 0 : dvardx0_m1, dvardx0_00, dvardx0_p1)
1399 0 : else if (i_equ == -6) then ! 'gradT'
1400 : call get_gradT_partials(s, &
1401 : k, i_var, i_var_sink, i_var_xa_index, i_var_sink_xa_index, &
1402 0 : dvardx0_m1, dvardx0_00, dvardx0_p1)
1403 0 : else if (i_equ == -7) then ! 'mlt_vc'
1404 : call get_mlt_vc_partials(s, &
1405 : k, i_var, i_var_sink, i_var_xa_index, i_var_sink_xa_index, &
1406 0 : dvardx0_m1, dvardx0_00, dvardx0_p1)
1407 0 : else if (i_equ == -8) then ! 'grad_ad'
1408 : call get_grad_ad_partials(s, &
1409 : k, i_var, i_var_sink, i_var_xa_index, i_var_sink_xa_index, &
1410 0 : dvardx0_m1, dvardx0_00, dvardx0_p1)
1411 : end if
1412 0 : if (k > 1) then
1413 : call test1_partial(s, &
1414 : i_equ, i_var, i_var_sink, i_var_xa_index, i_var_sink_xa_index, &
1415 0 : k, -1, dvardx0_m1, save_dx, save_equ, ierr)
1416 0 : if (ierr /= 0) call mesa_error(__FILE__,__LINE__,'test3_partials')
1417 : end if
1418 : call test1_partial(s, &
1419 : i_equ, i_var, i_var_sink, i_var_xa_index, i_var_sink_xa_index, &
1420 0 : k, 0, dvardx0_00, save_dx, save_equ, ierr)
1421 0 : if (ierr /= 0) call mesa_error(__FILE__,__LINE__,'test3_partials')
1422 0 : if (k < s% nz) then
1423 : call test1_partial(s, &
1424 : i_equ, i_var, i_var_sink, i_var_xa_index, i_var_sink_xa_index, &
1425 0 : k, 1, dvardx0_p1, save_dx, save_equ, ierr)
1426 0 : if (ierr /= 0) call mesa_error(__FILE__,__LINE__,'test3_partials')
1427 : end if
1428 0 : end subroutine test3_partials
1429 :
1430 :
1431 0 : subroutine test1_partial(s, &
1432 : i_equ, i_var, i_var_sink, i_var_xa_index, i_var_sink_xa_index, &
1433 : k, k_off, dvardx_0, save_dx, save_equ, ierr)
1434 : use chem_def, only: chem_isos
1435 : type (star_info), pointer :: s
1436 : integer, intent(in) :: &
1437 : i_equ, i_var, i_var_sink, i_var_xa_index, i_var_sink_xa_index, k, k_off
1438 : real(dp), intent(in) :: dvardx_0
1439 : real(dp), pointer, dimension(:,:) :: save_dx, save_equ
1440 : character (len=3) :: k_off_str
1441 : integer, intent(out) :: ierr
1442 : character (len = 32) :: equ_str
1443 0 : real(dp) :: dx_0, err, dvardx, xdum, uncertainty
1444 : include 'formats'
1445 0 : ierr = 0
1446 :
1447 0 : if (i_var > s% nvar_hydro) then ! testing abundance
1448 : dx_0 = s% solver_test_partials_dx_0 * &
1449 : max(abs(s% xa_start(i_var_xa_index,k) + s% solver_dx(i_var,k)), &
1450 : abs(s% xa_start(i_var_xa_index,k)), &
1451 0 : 1d-99)
1452 0 : write(*,1) 'var name ' // chem_isos% name(s% chem_id(i_var_xa_index))
1453 0 : write(*,1) 'sink name ' // chem_isos% name(s% chem_id(i_var_sink_xa_index))
1454 : else
1455 : dx_0 = s% solver_test_partials_dx_0 * &
1456 : max(abs(s% xh_start(i_var,k) + s% solver_dx(i_var,k)), &
1457 0 : abs(s% xh_start(i_var,k)))
1458 0 : if (dx_0 == 0d0) dx_0 = s% solver_test_partials_dx_0
1459 : end if
1460 : dvardx = dfridr(s, &
1461 : i_equ, i_var, i_var_sink, i_var_xa_index, i_var_sink_xa_index, &
1462 0 : k, k_off, dx_0, save_dx, err)
1463 0 : if (dvardx == 0d0 .and. abs(dvardx_0) < 1d-14) then
1464 0 : xdum = 0d0
1465 0 : else if (dvardx_0 == 0d0 .and. abs(dvardx) < 1d-14) then
1466 0 : xdum = 0d0
1467 0 : else if (dvardx == 0d0 .or. dvardx_0 == 0d0) then
1468 0 : xdum = 1d0
1469 : else
1470 0 : xdum = abs(dvardx - dvardx_0)/min(abs(dvardx),abs(dvardx_0))
1471 : end if
1472 0 : if (ierr /= 0) then
1473 0 : write(*,*) 'test1_partial failed'
1474 0 : call mesa_error(__FILE__,__LINE__,'setmatrix')
1475 : end if
1476 0 : if (i_equ /= 0) then
1477 0 : if (k_off == 0) then
1478 0 : k_off_str = ') '
1479 0 : else if (k_off == -1) then
1480 0 : k_off_str = '-1)'
1481 0 : else if (k_off == 1) then
1482 0 : k_off_str = '+1)'
1483 : end if
1484 0 : if (dvardx /= 0d0) then
1485 0 : uncertainty = abs(err/dvardx)
1486 : else
1487 0 : uncertainty = 0d0
1488 : end if
1489 0 : if (xdum > 1d-5 .and. uncertainty < 1d-6) then
1490 0 : write(*, '(a5,1x)', advance='no') '*****'
1491 0 : else if (uncertainty > 1d-7) then
1492 0 : write(*, '(a5,1x)', advance='no') '?????'
1493 : else
1494 0 : write(*, '(6x)', advance='no')
1495 : end if
1496 0 : if (i_equ > 0) then
1497 0 : equ_str = s% nameofequ(i_equ)
1498 : else if (i_equ == -1) then
1499 0 : equ_str = 'lnE'
1500 : else if (i_equ == -2) then
1501 0 : equ_str = 'eps_nuc'
1502 : else if (i_equ == -3) then
1503 0 : equ_str = 'opacity'
1504 : else if (i_equ == -4) then
1505 0 : equ_str = 'lnP'
1506 : else if (i_equ == -5) then
1507 0 : equ_str = 'non_nuc_neu'
1508 : else if (i_equ == -6) then
1509 0 : equ_str = 'gradT'
1510 : else if (i_equ == -7) then
1511 0 : equ_str = 'mlt_vc'
1512 : else if (i_equ == -8) then
1513 0 : equ_str = 'grad_ad'
1514 : else
1515 0 : equ_str = 'unknown'
1516 : end if
1517 : write(*,'(a25,3x,i5,4x,a,f12.5,4x,a,f12.5,99(4x,a,1pe22.12))') &
1518 : 'd_' // trim(equ_str) // '(k)/d_' // trim(s% nameofvar(i_var)) // &
1519 0 : '(k' // trim(k_off_str), &
1520 0 : k, 'lg rel diff', safe_log10(xdum), 'lg uncertainty', safe_log10(uncertainty), &
1521 0 : 'Analytic', dvardx_0, 'Numeric', dvardx
1522 : ! write(*,'(a70,2x,i5,f10.3,3x,a,f10.3,99(3x,a,1pe26.16))') &
1523 : ! 'log dfridr rel_diff partials wrt ' // trim(s% nameofvar(i_var)) // &
1524 : ! '(k' // k_off_str // ' of ' // trim(equ_str) // '(k)', &
1525 : ! k, safe_log10(xdum), 'log uncertainty', safe_log10(uncertainty), &
1526 : ! 'analytic', dvardx_0, 'numeric', dvardx, &
1527 : ! 'analytic/numeric', abs(dvardx_0)/max(1d-99,abs(dvardx))
1528 :
1529 : else
1530 0 : write(*,'(A)')
1531 0 : write(*,1) 'analytic and numeric partials wrt ' // trim(s% nameofvar(i_var)), &
1532 0 : dvardx_0, dvardx
1533 0 : write(*,1) 'log dfridr relative uncertainty for numeric partial', &
1534 0 : safe_log10(err/max(1d-99,abs(dvardx)))
1535 0 : if (dvardx_0 /= 0d0) write(*,1) 'rel_diff', xdum
1536 : end if
1537 0 : end subroutine test1_partial
1538 :
1539 :
1540 0 : real(dp) function dfridr_func(s, &
1541 : i_equ, i_var, i_var_sink, i_var_xa_index, i_var_sink_xa_index, &
1542 0 : k, k_off, delta_x, save_dx) result(val)
1543 : type (star_info), pointer :: s
1544 : integer, intent(in) :: &
1545 : i_equ, i_var, i_var_sink, i_var_xa_index, i_var_sink_xa_index, k, k_off
1546 : real(dp), intent(in) :: delta_x
1547 : real(dp), pointer, dimension(:,:) :: save_dx
1548 : include 'formats'
1549 0 : s% solver_dx(i_var,k+k_off) = save_dx(i_var,k+k_off) + delta_x
1550 0 : if (i_var_xa_index > 0) then ! changing abundance
1551 : !write(*,2) 'new dx, x for abundance', i_var, &
1552 : ! dx(i_var,k+k_off), s% xa(i_var - s% nvar_hydro,k+k_off)
1553 0 : if (i_var_sink_xa_index <= 0 .or. i_var_sink_xa_index > s% species) then
1554 0 : write(*,2) 'bad i_var_sink_xa_index', i_var_sink_xa_index
1555 0 : call mesa_error(__FILE__,__LINE__,'star_solver dfridr_func')
1556 : end if
1557 0 : s% solver_dx(i_var_sink,k+k_off) = save_dx(i_var_sink,k+k_off) - delta_x
1558 : end if
1559 0 : call do_equations(ierr)
1560 0 : if (ierr /= 0) then
1561 : !exit
1562 0 : write(*,3) 'call eval_equations failed in dfridr_func'
1563 0 : call mesa_error(__FILE__,__LINE__,'setmatrix')
1564 : end if
1565 0 : if (i_equ > 0) then
1566 0 : val = equ(i_equ,k) ! testing partial of residual for cell k equation
1567 : else if (i_equ == 0) then
1568 0 : val = s% solver_test_partials_val
1569 : else if (i_equ == -1) then
1570 0 : val = s% lnE(k)
1571 : else if (i_equ == -2) then
1572 0 : val = s% eps_nuc(k)
1573 : else if (i_equ == -3) then
1574 0 : val = s% opacity(k)
1575 : else if (i_equ == -4) then
1576 0 : val = s% lnPeos(k)
1577 : else if (i_equ == -5) then
1578 0 : val = s% non_nuc_neu(k)
1579 : else if (i_equ == -6) then
1580 0 : val = s% gradT(k)
1581 : else if (i_equ == -7) then
1582 0 : val = s% mlt_vc(k)
1583 : else if (i_equ == -8) then
1584 0 : val = s% grada(k)
1585 : else
1586 : val = 0d0
1587 : end if
1588 0 : s% solver_dx(i_var,k+k_off) = save_dx(i_var,k+k_off)
1589 0 : if (i_var_sink > 0) & ! restore sink abundance
1590 0 : s% solver_dx(i_var_sink,k+k_off) = save_dx(i_var_sink,k+k_off)
1591 0 : end function dfridr_func
1592 :
1593 :
1594 0 : real(dp) function dfridr(s, &
1595 : i_equ, i_var, i_var_sink, i_var_xa_index, i_var_sink_xa_index, &
1596 : k, k_off, hx, save_dx, err)
1597 : type (star_info), pointer :: s
1598 : integer, intent(in) :: &
1599 : i_equ, i_var, i_var_sink, i_var_xa_index, i_var_sink_xa_index, k, k_off
1600 : real(dp), intent(in) :: hx
1601 : real(dp), pointer, dimension(:,:) :: save_dx
1602 : real(dp), intent(out) :: err
1603 : ! this routine returns the first derivative of a function func(x)
1604 : ! at the point x, by ridders method of polynomial extrapolation.
1605 : ! value hx is the initial step size;
1606 : ! it should be an increment for which func changes substantially.
1607 : ! an estimate of the error in the first derivative is returned in err.
1608 : integer, parameter :: ntab = 20
1609 : integer :: i,j
1610 0 : real(dp) :: errt,fac,hh,a(ntab,ntab),f1,f2
1611 : real(dp), parameter :: con2=2d0, con=sqrt(con2), big=1d50, safe=2d0
1612 : include 'formats'
1613 0 : dfridr = 0d0
1614 0 : hh = hx
1615 : ! 2nd order central difference
1616 : f1 = dfridr_func(s, &
1617 : i_equ, i_var, i_var_sink, i_var_xa_index, i_var_sink_xa_index, &
1618 0 : k, k_off, hh, save_dx)
1619 : !write(*,2) 'f1', 1, f1, save_dx(i_var,k) + hh
1620 : f2 = dfridr_func(s, &
1621 : i_equ, i_var, i_var_sink, i_var_xa_index, i_var_sink_xa_index, &
1622 0 : k, k_off, -hh, save_dx)
1623 : !write(*,2) 'f2', 1, f2, save_dx(i_var,k) - hh
1624 0 : a(1,1) = (f1 - f2)/(2d0*hh)
1625 : !write(*,2) 'dfdx', 1, a(1,1), &
1626 : ! hh, (save_dx(s% solver_test_partials_var,s% solver_test_partials_k) + hh)/ln10, &
1627 : ! save_dx(s% solver_test_partials_var,s% solver_test_partials_k)/ln10
1628 0 : err = big
1629 : ! successive columns in the neville tableu will go to smaller stepsizes
1630 : ! and higher orders of extrapolation
1631 0 : do i=2,ntab
1632 0 : hh = hh/con
1633 : f1 = dfridr_func(s, &
1634 : i_equ, i_var, i_var_sink, i_var_xa_index, i_var_sink_xa_index, &
1635 0 : k, k_off, hh, save_dx)
1636 : !write(*,2) 'f1', i, f1, save_dx(i_var,k) + hh
1637 : f2 = dfridr_func(s, &
1638 : i_equ, i_var, i_var_sink, i_var_xa_index, i_var_sink_xa_index, &
1639 0 : k, k_off, -hh, save_dx)
1640 : !write(*,2) 'f2', i, f2, save_dx(i_var,k) - hh
1641 0 : a(1,i) = (f1 - f2)/(2d0*hh)
1642 : !write(*,2) 'dfdx', i, a(1,i), &
1643 : ! hh, (save_dx(s% solver_test_partials_var,s% solver_test_partials_k) + hh)/ln10, &
1644 : ! save_dx(s% solver_test_partials_var,s% solver_test_partials_k)/ln10
1645 : ! compute extrapolations of various orders; the error strategy is to compare
1646 : ! each new extrapolation to one order lower but both at the same stepsize
1647 : ! and at the previous stepsize
1648 0 : fac = con2
1649 0 : do j=2,i
1650 0 : a(j,i) = (a(j-1,i)*fac - a(j-1,i-1))/(fac-1d0)
1651 0 : fac = con2*fac
1652 0 : errt = max(abs(a(j,i)-a(j-1,i)),abs(a(j,i)-a(j-1,i-1)))
1653 0 : if (errt <= err) then
1654 0 : err = errt
1655 0 : dfridr = a(j,i)
1656 : !write(*,1) 'dfridr', err/dfridr
1657 : end if
1658 : end do
1659 : ! if higher order is worse by a significant factor safe, then bail
1660 0 : if (abs(a(i,i) - a(i-1,i-1)) >= safe*err) then
1661 : !write(*,1) 'higher order is worse'
1662 0 : return
1663 : end if
1664 : end do
1665 : end function dfridr
1666 :
1667 :
1668 : subroutine set_xtras(x,num_xtra)
1669 : real(dp) :: x(:,:)
1670 : integer, intent(in) :: num_xtra
1671 : integer :: k
1672 : include 'formats'
1673 : if (.not. s% u_flag) then
1674 : x(1,1:nz) = 0
1675 : x(2,1:nz) = 0
1676 : return
1677 : end if
1678 : do k=1,nz
1679 : x(1,k) = s% u_face_ad(k)%val
1680 : if (is_bad_num(x(1,k))) then
1681 : write(*,2) 'exit_setmatrix x(1,k)', k, x(1,k)
1682 : stop
1683 : end if
1684 : end do
1685 : do k=1,nz
1686 : x(2,k) = s% P_face_ad(k)%val
1687 : if (is_bad_num(x(2,k))) then
1688 : write(*,2) 'exit_setmatrix x(2,k)', k, x(2,k)
1689 : stop
1690 : end if
1691 : end do
1692 : end subroutine set_xtras
1693 :
1694 :
1695 0 : subroutine store_mix_type_str(str, integer_string, i, k)
1696 : character (len=5) :: str
1697 : character (len=10) :: integer_string
1698 : integer, intent(in) :: i, k
1699 : integer :: mix_type, j
1700 0 : if (k < 1 .or. k > s% nz) then
1701 0 : str(i:i) = 'x'
1702 0 : return
1703 : end if
1704 0 : mix_type = s% mixing_type(k)
1705 0 : if (mix_type < 10) then
1706 0 : j = mix_type+1
1707 0 : str(i:i) = integer_string(j:j)
1708 : else
1709 0 : str(i:i) = '?'
1710 : end if
1711 : end subroutine store_mix_type_str
1712 :
1713 :
1714 0 : subroutine write_msg(msg)
1715 : use const_def, only: secyer
1716 : character(*) :: msg
1717 :
1718 : integer :: k
1719 : character (len=64) :: max_resid_str, max_corr_str
1720 : character (len=5) :: max_resid_mix_type_str, max_corr_mix_type_str
1721 : character (len=10) :: integer_string
1722 : include 'formats'
1723 :
1724 0 : if (.not. dbg_msg) return
1725 :
1726 0 : if (max_resid_j < 0) then
1727 0 : call sizequ(s, nvar, residual_norm, max_residual, max_resid_k, max_resid_j, ierr)
1728 : end if
1729 :
1730 0 : if (max_resid_j > 0) then
1731 0 : write(max_resid_str,*) 'max resid ' // trim(s% nameofequ(max_resid_j))
1732 : else
1733 0 : max_resid_str = ''
1734 : end if
1735 :
1736 0 : if (max_corr_j < 0) then
1737 : call sizeB(s, nvar, B, &
1738 0 : max_correction, correction_norm, max_corr_k, max_corr_j, ierr)
1739 : end if
1740 :
1741 0 : if (max_corr_j > 0) then
1742 0 : write(max_corr_str,*) 'max corr ' // trim(s% nameofvar(max_corr_j))
1743 : else
1744 0 : max_corr_str = ''
1745 : end if
1746 :
1747 0 : integer_string = '0123456789'
1748 0 : k = max_corr_k
1749 0 : call store_mix_type_str(max_corr_mix_type_str, integer_string, 1, k-2)
1750 0 : call store_mix_type_str(max_corr_mix_type_str, integer_string, 2, k-1)
1751 0 : call store_mix_type_str(max_corr_mix_type_str, integer_string, 3, k)
1752 0 : call store_mix_type_str(max_corr_mix_type_str, integer_string, 4, k+1)
1753 0 : call store_mix_type_str(max_corr_mix_type_str, integer_string, 5, k+2)
1754 :
1755 0 : k = max_resid_k
1756 0 : call store_mix_type_str(max_resid_mix_type_str, integer_string, 1, k-2)
1757 0 : call store_mix_type_str(max_resid_mix_type_str, integer_string, 2, k-1)
1758 0 : call store_mix_type_str(max_resid_mix_type_str, integer_string, 3, k)
1759 0 : call store_mix_type_str(max_resid_mix_type_str, integer_string, 4, k+1)
1760 0 : call store_mix_type_str(max_resid_mix_type_str, integer_string, 5, k+2)
1761 :
1762 : 111 format(i6, i3, 2x, a, f7.4, &
1763 : 1x, a, 1x, e10.3, 2x, a18, 1x, i5, e13.5, a, &
1764 : 1x, a, 1x, e10.3, 2x, a16, 1x, i5, e13.5, a, &
1765 : 1x, a)
1766 : ! 111 format(i6, 2x, i3, 2x, a, f8.4, &
1767 : ! 2x, a, 1x, e10.3, 2x, a19, 1x, i5, e11.3, 2x, a, &
1768 : ! 2x, a, 1x, e10.3, 2x, a14, 1x, i5, e11.3, 2x, a, &
1769 : ! 2x, a)
1770 : write(*,111) &
1771 0 : s% model_number, iter, &
1772 0 : 'coeff', coeff, &
1773 0 : 'avg resid', residual_norm, &
1774 : ! ' avg resid', residual_norm, &
1775 0 : trim(max_resid_str), max_resid_k, max_residual, &
1776 0 : ' mix type ' // trim(max_resid_mix_type_str), &
1777 0 : 'avg corr', correction_norm, &
1778 : ! 'mix type ' // trim(max_resid_mix_type_str), &
1779 : ! ' avg corr', correction_norm, &
1780 0 : trim(max_corr_str), max_corr_k, max_correction, &
1781 0 : ' mix type ' // trim(max_corr_mix_type_str), &
1782 0 : ' ' // trim(msg)
1783 : ! 'mix type ' // trim(max_corr_mix_type_str), &
1784 : ! ' ' // trim(msg)
1785 :
1786 0 : if (is_bad(slope)) call mesa_error(__FILE__,__LINE__,'write_msg')
1787 :
1788 : end subroutine write_msg
1789 :
1790 :
1791 11 : subroutine pointers(ierr)
1792 : use utils_lib, only: fill_with_NaNs, fill_with_NaNs_2D
1793 :
1794 : integer, intent(out) :: ierr
1795 :
1796 : integer :: i
1797 :
1798 11 : ierr = 0
1799 :
1800 11 : i = 1
1801 :
1802 11 : if(allocated(s% equ1)) then
1803 10 : if(size(s% equ1,dim=1) /= neq) then
1804 4 : deallocate(s% equ1)
1805 4 : allocate(s% equ1(1:neq))
1806 : end if
1807 : else
1808 1 : allocate(s% equ1(1:neq))
1809 : end if
1810 :
1811 11 : s% equ(1:nvar,1:nz) => s% equ1(1:neq)
1812 11 : equ1 => s% equ1
1813 11 : equ => s% equ
1814 :
1815 11 : allocate(dxsave1(1:neq))
1816 11 : allocate(ddxsave1(1:neq))
1817 11 : allocate(B1(1:neq))
1818 11 : allocate(soln1(1:neq))
1819 11 : allocate(grad_f1(1:neq))
1820 22 : allocate(rhs(1:nvar,1:nz))
1821 22 : allocate(xder(1:nvar,1:nz))
1822 22 : allocate(ddx(1:nvar,1:nz))
1823 11 : allocate(row_scale_factors1(1:neq))
1824 11 : allocate(col_scale_factors1(1:neq))
1825 11 : allocate(save_ublk1(1:nvar*neq))
1826 11 : allocate(save_dblk1(1:nvar*neq))
1827 11 : allocate(save_lblk1(1:nvar*neq))
1828 :
1829 11 : if (s% fill_arrays_with_NaNs) then
1830 0 : call fill_with_NaNs(dxsave1)
1831 0 : call fill_with_NaNs(ddxsave1)
1832 0 : call fill_with_NaNs(B1)
1833 0 : call fill_with_NaNs(soln1)
1834 0 : call fill_with_NaNs_2D(rhs)
1835 0 : call fill_with_NaNs_2D(xder)
1836 0 : call fill_with_NaNs_2D(ddx)
1837 0 : call fill_with_NaNs(row_scale_factors1)
1838 0 : call fill_with_NaNs(col_scale_factors1)
1839 0 : call fill_with_NaNs(save_ublk1)
1840 0 : call fill_with_NaNs(save_dblk1)
1841 0 : call fill_with_NaNs(save_lblk1)
1842 : end if
1843 :
1844 11 : dxsave(1:nvar,1:nz) => dxsave1(1:neq)
1845 11 : ddxsave(1:nvar,1:nz) => ddxsave1(1:neq)
1846 11 : B(1:nvar,1:nz) => B1(1:neq)
1847 11 : soln(1:nvar,1:nz) => soln1(1:neq)
1848 11 : grad_f(1:nvar,1:nz) => grad_f1(1:neq)
1849 :
1850 11 : allocate(ipiv1(1:neq))
1851 :
1852 11 : ipiv_blk1(1:neq) => ipiv1(1:neq)
1853 :
1854 11 : allocate(ublk1(1:nvar*neq),dblk1(1:nvar*neq),lblk1(1:nvar*neq))
1855 :
1856 11 : lblk(1:nvar,1:nvar,1:nz) => lblk1(1:nvar*neq)
1857 11 : dblk(1:nvar,1:nvar,1:nz) => dblk1(1:nvar*neq)
1858 11 : ublk(1:nvar,1:nvar,1:nz) => ublk1(1:nvar*neq)
1859 :
1860 11 : allocate(ublkF1(1:nvar*neq),dblkF1(1:nvar*neq),lblkF1(1:nvar*neq))
1861 :
1862 11 : lblkF(1:nvar,1:nvar,1:nz) => lblkF1(1:nvar*neq)
1863 11 : dblkF(1:nvar,1:nvar,1:nz) => dblkF1(1:nvar*neq)
1864 11 : ublkF(1:nvar,1:nvar,1:nz) => ublkF1(1:nvar*neq)
1865 :
1866 11 : if (s% fill_arrays_with_NaNs) then
1867 0 : call fill_with_NaNs(lblk1)
1868 0 : call fill_with_NaNs(dblk1)
1869 0 : call fill_with_NaNs(ublk1)
1870 0 : call fill_with_NaNs(lblkF1)
1871 0 : call fill_with_NaNs(dblkF1)
1872 0 : call fill_with_NaNs(ublkF1)
1873 : end if
1874 :
1875 11 : end subroutine pointers
1876 :
1877 11 : subroutine cleanup()
1878 :
1879 11 : if(associated(dxsave1)) deallocate(dxsave1)
1880 11 : if(associated(ddxsave)) deallocate(ddxsave)
1881 11 : if(associated(ddxsave)) deallocate(ddxsave)
1882 11 : if(associated(B1)) deallocate(B1)
1883 11 : if(associated(soln1)) deallocate(soln1)
1884 11 : if(associated(grad_f1)) deallocate(grad_f1)
1885 11 : if(associated(rhs)) deallocate(rhs)
1886 11 : if(associated(xder)) deallocate(xder)
1887 11 : if(associated(ddx)) deallocate(ddx)
1888 11 : if(associated(row_scale_factors1)) deallocate(row_scale_factors1)
1889 11 : if(associated(col_scale_factors1)) deallocate(col_scale_factors1)
1890 11 : if(associated(save_ublk1)) deallocate(save_ublk1)
1891 11 : if(associated(save_dblk1)) deallocate(save_dblk1)
1892 11 : if(associated(save_lblk1)) deallocate(save_lblk1)
1893 11 : if(associated(ipiv1)) deallocate(ipiv1)
1894 :
1895 11 : if(associated(ublk1)) deallocate(ublk1)
1896 11 : if(associated(dblk1)) deallocate(dblk1)
1897 11 : if(associated(lblk1)) deallocate(lblk1)
1898 11 : if(associated(ublkF1)) deallocate(ublkF1)
1899 11 : if(associated(dblkF1)) deallocate(dblkF1)
1900 11 : if(associated(lblkF1)) deallocate(lblkF1)
1901 :
1902 :
1903 11 : nullify(equ, equ1, dxsave1,dxsave, ddxsave, B1, &
1904 11 : soln1, grad_f1, rhs, xder, ddx, row_scale_factors1,&
1905 11 : col_scale_factors1, save_ublk1, save_dblk1, save_lblk1,&
1906 11 : B, soln, grad_f,ipiv1, ublk1, dblk1, lblk1, ublkF1,dblkF1, lblkF1)
1907 :
1908 11 : end subroutine cleanup
1909 :
1910 :
1911 0 : real(dp) function eval_slope(nvar, nz, grad_f, B)
1912 : integer, intent(in) :: nvar, nz
1913 : real(dp), intent(in), dimension(:,:) :: grad_f, B
1914 : integer :: i
1915 0 : eval_slope = 0
1916 0 : do i=1,nvar
1917 0 : eval_slope = eval_slope + dot_product(grad_f(i,1:nz),B(i,1:nz))
1918 : end do
1919 0 : end function eval_slope
1920 :
1921 :
1922 0 : real(dp) function eval_f(nvar, nz, equ)
1923 : integer, intent(in) :: nvar, nz
1924 : real(dp), intent(in), dimension(:,:) :: equ
1925 : integer :: k, i
1926 0 : real(dp) :: q
1927 : include 'formats'
1928 0 : eval_f = 0
1929 0 : do k = 1, nz
1930 0 : do i = 1, nvar
1931 0 : q = equ(i,k)
1932 0 : eval_f = eval_f + q*q
1933 : end do
1934 : end do
1935 0 : eval_f = eval_f/2
1936 0 : end function eval_f
1937 :
1938 : end subroutine do_solver
1939 :
1940 : end module star_solver
|
