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