Line data Source code
1 : ! ***********************************************************************
2 : !
3 : ! Copyright (C) 2012-2019 The MESA Team
4 : !
5 : ! This program is free software: you can redistribute it and/or modify
6 : ! it under the terms of the GNU Lesser General Public License
7 : ! as published by the Free Software Foundation,
8 : ! either version 3 of the License, or (at your option) any later version.
9 : !
10 : ! This program is distributed in the hope that it will be useful,
11 : ! but WITHOUT ANY WARRANTY; without even the implied warranty of
12 : ! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
13 : ! See the GNU Lesser General Public License for more details.
14 : !
15 : ! You should have received a copy of the GNU Lesser General Public License
16 : ! along with this program. If not, see <https://www.gnu.org/licenses/>.
17 : !
18 : ! ***********************************************************************
19 :
20 : module diffusion
21 :
22 : use const_def, only: dp, i8, amu, me, msun
23 : use chem_def
24 : use star_private_def
25 : use diffusion_support
26 : use diffusion_procs
27 :
28 : implicit none
29 :
30 : private
31 : public :: do_solve_diffusion
32 : public :: set_diffusion_classes
33 : public :: diffusion_min_nc
34 :
35 : integer, parameter :: diffusion_min_nc = 4 ! minimum number of classes
36 : logical, parameter :: use_dcoeff_dX = .true.
37 :
38 : contains
39 :
40 0 : subroutine do_solve_diffusion( &
41 0 : s, nz, species, nc, m, class, class_chem_id, net_iso, chem_id, &
42 0 : abar, ye, free_e, dm_bar_in, dm_in, &
43 0 : T, lnT, rho, lnd, r_mid, dlnP_dm, dlnT_dm, dlnRho_dm, &
44 0 : L_face, r_face_in, dlnP_dm_face, dlnT_dm_face, dlnRho_dm_face, &
45 : pure_Coulomb, total_time, dt_div_timescale, &
46 : max_steps, max_iters_total, max_iters_per_substep, &
47 0 : calculate_ionization, typical_charge, nsmooth_typical_charge, &
48 : min_T_for_radaccel, max_T_for_radaccel, &
49 : min_Z_for_radaccel, max_Z_for_radaccel, &
50 : screening_for_radaccel, op_mono_data_path, op_mono_data_cache_filename, &
51 : v_advection_max, R_center, &
52 0 : gamma, gamma_full_on, gamma_full_off, &
53 0 : T_full_on, T_full_off, diffusion_factor, &
54 0 : xa, steps_used, total_num_iters, total_num_retries, &
55 0 : nzlo, nzhi, X_init, X_final, &
56 0 : D_self_face, v_advection_face, v_total_face, &
57 0 : vlnP_face, vlnT_face, v_rad_face, g_rad_face, &
58 0 : E_field_face, g_field_face, &
59 : ierr)
60 :
61 : use diffusion_procs, only: fixup
62 : use kap_lib, only: load_op_mono_data
63 :
64 : type (star_info), pointer :: s
65 : integer, intent(in) :: nz, species, nc, m, &
66 : min_Z_for_radaccel, max_Z_for_radaccel
67 : ! nc = number of classes of species
68 : ! m = nc+1
69 : integer, intent(in) :: class(:) ! (species)
70 : integer, intent(in) :: class_chem_id(:) ! (nc)
71 : integer, intent(in) :: net_iso(:), chem_id(:)
72 : ! class(i) = class number for species i. class numbers from 1 to nc
73 : ! class_chem_id(j) = isotope id number from chem_def for "typical" member of class j
74 : real(dp), intent(in), dimension(:) :: &
75 : abar, ye, free_e, gamma, dm_bar_in, dm_in, &
76 : T, lnT, rho, lnd, r_mid, dlnP_dm, dlnT_dm, dlnRho_dm, &
77 : L_face, r_face_in, dlnP_dm_face, dlnT_dm_face, dlnRho_dm_face
78 : logical, intent(in) :: pure_Coulomb, screening_for_radaccel
79 : character (len=*), intent(in) :: &
80 : op_mono_data_path, op_mono_data_cache_filename
81 : real(dp), intent(in) :: &
82 : total_time, dt_div_timescale, &
83 : min_T_for_radaccel, max_T_for_radaccel, &
84 : v_advection_max, R_center, gamma_full_on, gamma_full_off, &
85 : T_full_on, T_full_off, diffusion_factor(nc)
86 : integer, intent(in) :: max_steps, max_iters_total, max_iters_per_substep
87 : logical, intent(in) :: calculate_ionization
88 : real(dp), intent(inout) :: typical_charge(:,:) ! (nc,nz)
89 : integer, intent(in) :: nsmooth_typical_charge
90 : real(dp), intent(inout) :: xa(:,:) ! (species,nz) ! mass fractions
91 : real(dp), intent(out), dimension(:,:) :: X_init, X_final, &
92 : D_self_face, v_advection_face, v_total_face, &
93 : vlnP_face, vlnT_face, v_rad_face, g_rad_face ! (nc,nz)
94 : real(dp), intent(inout) :: E_field_face(:) ! (nz)
95 : real(dp), intent(inout) :: g_field_face(:) ! (nz)
96 : integer, intent(out) :: steps_used, total_num_iters, total_num_retries, ierr
97 : integer, intent(inout) :: nzlo, nzhi !upper and lower bounds on region
98 :
99 : integer :: i, j, k, num_iters, idiag, n, neqs, ku, kl, &
100 : ldab, ldafb, ldb, ldx, h1, he4, nbound, bad_j, bad_k, &
101 : retry_count, max_retries, ierr_dealloc
102 : real(dp) :: &
103 0 : dt, min_dt, time, mstar, mtotal, sum_mass_nzlo_nzhi, xtotal_init(nc), &
104 0 : mass_init(nc), bad_X, bad_sum, bad_Xsum, &
105 0 : tol_correction_max, tol_correction_norm, remaining_time, &
106 0 : sumx, tmp, timescale, min_lambda, upwind_limit, lambda, max_del, avg_del
107 :
108 : integer, parameter :: min_nz_lo = 5
109 : real(dp), parameter :: &
110 : dt_retry_factor = 0.5d0, dt_max_factor = 2d0, dt_min_factor = 0.75d0, &
111 : tiny_X = 1d-50, tiny_C = 1d-50, max_sum_abs = 10, X_cleanup_tol = 1d-2, &
112 : max_flow_frac = 1d0, max_flow_X_limit = 1d-5, dx_avg_target = 0.975d0
113 :
114 0 : real(dp), dimension(:), allocatable :: cell_dm, dm_bar, sum_new_mass
115 : real(dp), dimension(:,:), allocatable :: &
116 0 : Z, X, ending_dX_dm, C, dC_dr, C_div_X, &
117 0 : new_mass, X_0, X_1, dX, dX_dt
118 :
119 : real(dp), dimension(:), allocatable :: &
120 0 : r_face, alfa_face, rho_face, T_face, four_pi_r2_rho_face, &
121 0 : dlnP_dr_face, dlnT_dr_face, dlnRho_dr_face, sum_starting_mass, &
122 0 : limit_coeffs_face, &
123 0 : e_ap, e_at, e_ar, &
124 0 : g_ap, g_at, g_ar
125 :
126 : real(dp), dimension(:), allocatable :: &
127 0 : AD_face, xm_face, sum_ending_mass
128 : real(dp), dimension(:,:), allocatable :: &
129 0 : X_face, C_div_X_face, GT_face, &
130 0 : rad_accel_face, log10_g_rad, &
131 0 : ending_mass, starting_mass, starting_dX_dm, &
132 0 : X_theta, X_minus_theta, e_ax, g_ax
133 : real(dp), dimension(:,:,:), allocatable :: &
134 0 : SIG_face, sigma_lnC
135 :
136 0 : real(dp), pointer, dimension(:) :: del1, rhs1, lblk1, dblk1, ublk1, b1
137 0 : real(dp), dimension(:,:), pointer :: rhs, del
138 0 : real(dp), dimension(:,:,:), pointer :: em1, e00, ep1
139 0 : integer, pointer :: ipiv1(:)
140 : integer :: lrd, lid, kmax_rad_accel, min_num_substeps, &
141 : iter_dbg, j_dbg, k_dbg, k_max
142 : integer(i8) :: time0, time1, clock_rate
143 0 : integer, pointer :: ipar_decsol(:)
144 0 : real(dp), pointer :: rpar_decsol(:)
145 0 : real(dp), dimension(species) :: xa_total_before
146 0 : real(dp), dimension(m) :: A
147 0 : real(dp) :: X_total_atol, X_total_rtol, vc_old, dt_old
148 : logical :: have_changed_matrix_coeffs, trace, last_step, &
149 : solved, do_timing, use_isolve
150 :
151 :
152 : logical, parameter :: dbg = .false.
153 : !logical, parameter :: dbg = .true.
154 :
155 :
156 : include 'formats'
157 :
158 0 : steps_used = 0
159 0 : total_num_iters = 0
160 0 : total_num_retries = 0
161 :
162 0 : use_isolve = s% diffusion_use_isolve
163 :
164 0 : min_num_substeps = s% diffusion_min_num_substeps
165 0 : trace = s% show_diffusion_substep_info
166 0 : do_timing = s% show_diffusion_timing
167 :
168 : if (dbg) then
169 : write(*,2) 'nzlo', nzlo
170 : write(*,2) 'nzhi', nzhi
171 : write(*,2) 'nz', nz
172 : write(*,'(A)')
173 : end if
174 :
175 :
176 0 : min_lambda = 1d-2
177 :
178 0 : if (do_timing) then
179 0 : call system_clock(time0,clock_rate)
180 : else
181 : time0 = 0
182 : end if
183 :
184 0 : ierr = 0
185 0 : ierr_dealloc = 0
186 0 : if (m /= nc+1) then
187 0 : ierr = -1
188 0 : write(*,*) 'm /= nc+1'
189 0 : return
190 : end if
191 :
192 0 : if (T(1) <= max_T_for_radaccel .and. s% op_mono_method == 'hu') then
193 : if (dbg) write(*,*) 'call load_op_mono_data'
194 : call load_op_mono_data( &
195 0 : op_mono_data_path, op_mono_data_cache_filename, ierr)
196 : if (dbg) write(*,*) 'done load_op_mono_data'
197 0 : if (ierr /= 0) then
198 0 : write(*,*) 'error while loading OP data, ierr = ',ierr
199 0 : return
200 : end if
201 : end if
202 :
203 0 : X_total_atol = s% diffusion_X_total_atol
204 0 : X_total_rtol = s% diffusion_X_total_rtol
205 :
206 0 : h1 = net_iso(ih1)
207 0 : if (h1 == 0) then
208 0 : ierr = -1; write(*,*) 'isos must include h1 for diffusion'; return
209 : end if
210 :
211 0 : he4 = net_iso(ihe4)
212 0 : if (he4 == 0) then
213 0 : ierr = -1; write(*,*) 'isos must include he4 for diffusion'; return
214 : end if
215 :
216 : !reset nzlo if necessary
217 0 : nbound=nzlo
218 0 : do k=nzlo,nzhi
219 0 : if (T(k) > T_full_off) then
220 : nbound=k; exit
221 : end if
222 : end do
223 0 : nzlo=nbound
224 :
225 : allocate( &
226 0 : e_ap(nz), e_at(nz), e_ar(nz), e_ax(m,nz), & ! for e field
227 0 : g_ap(nz), g_at(nz), g_ar(nz), g_ax(m,nz), & ! for g field
228 : new_mass(nc,nz), starting_dX_dm(nc,nz), starting_mass(nc,nz), &
229 0 : ending_mass(nc,nz), ending_dX_dm(nc,nz), GT_face(nc,nz), AD_face(nz), &
230 0 : X_theta(nc,nz), X_minus_theta(nc,nz), SIG_face(nc,nc,nz), sigma_lnC(nc,nc,nz), &
231 0 : r_face(nz+1), alfa_face(nz), rho_face(nz), T_face(nz), four_pi_r2_rho_face(nz), &
232 0 : dlnP_dr_face(nz), dlnT_dr_face(nz), dlnRho_dr_face(nz), &
233 0 : limit_coeffs_face(nz), sum_starting_mass(nz), &
234 0 : cell_dm(nz), dm_bar(nz), sum_new_mass(nz), &
235 0 : Z(m,nz), X(m,nz), C(m,nz), dC_dr(m,nz), C_div_X(m,nz), &
236 0 : X_0(nc,nz), X_1(nc,nz), dX(nc,nz), dX_dt(nc,nz), &
237 0 : xm_face(nz), sum_ending_mass(nz), &
238 0 : X_face(m,nz), C_div_X_face(m,nz), &
239 0 : rad_accel_face(m,nz), log10_g_rad(m,nz))
240 :
241 : call get_limit_coeffs( &
242 : s, nz, nzlo, nzhi, &
243 : gamma_full_on, gamma_full_off, &
244 : T_full_on, T_full_off, &
245 0 : gamma, T, limit_coeffs_face, k_max)
246 : if (dbg) write(*,4) 'k_max for limit_coeffs_face /= 0', &
247 : k_max, nzhi, nzhi - k_max
248 0 : nzhi = k_max ! max k s.t. limit_coeffs_face(k) > 0
249 :
250 0 : n = nzhi-nzlo+1
251 :
252 : if (dbg) write(*,2) ' nz', nz
253 : if (dbg) write(*,2) 'nzlo, q', nzlo, s% q(nzlo)
254 : if (dbg) write(*,2) 'nzhi, q', nzhi, s% q(nzhi)
255 : if (dbg) write(*,2) ' n', n
256 0 : if (n <= 1) return
257 :
258 0 : neqs = n*nc
259 0 : idiag = 2*nc
260 0 : ku = 2*nc - 1
261 0 : kl = ku
262 0 : ldab = ku + kl + 1
263 0 : ldafb = kl + ldab
264 0 : ldb = neqs
265 0 : ldx = neqs
266 :
267 0 : upwind_limit = s% diffusion_upwind_abs_v_limit
268 :
269 0 : mstar = sum(dm_in(1:nz))
270 0 : do i=1,species
271 : xa_total_before(i) = &
272 0 : dot_product(dm_in(1:nz),xa(i,1:nz))/mstar
273 : end do
274 :
275 0 : call do_alloc(ierr)
276 0 : if (ierr /= 0) then
277 0 : if (dbg .or. s% report_ierr) write(*,*) 'diffusion failed in do_alloc'
278 0 : return
279 : end if
280 :
281 0 : do k=2,nz
282 0 : alfa_face(k) = dm_in(k-1)/(dm_in(k) + dm_in(k-1))
283 : end do
284 0 : alfa_face(1) = 1d0
285 : ! e.g., T_face(k) = alfa_face(k)*T(k) + (1d0-alfa_face(k))*T(k-1)
286 :
287 0 : do k=nzlo+1,nzhi
288 0 : r_face(k) = r_face_in(k)
289 : end do
290 0 : r_face(nzlo) = r_face_in(1)
291 0 : r_face(nzhi+1) = R_center
292 :
293 : ! create classes
294 0 : A(1:nc) = chem_isos% Z_plus_N(class_chem_id(1:nc))
295 0 : A(m) = me/amu
296 0 : do k=1,nz
297 0 : cell_dm(k) = dm_in(k)
298 0 : dm_bar(k) = dm_bar_in(k)
299 0 : do j=1,nc
300 0 : X(j,k) = 0
301 : end do
302 0 : do j=1,species
303 0 : i = class(j)
304 0 : X(i,k) = X(i,k) + max(tiny_X, xa(j,k))
305 : end do
306 0 : do i=1,nc
307 0 : if (X(i,k) <= 0) X(i,k) = tiny_X
308 : end do
309 0 : tmp = sum(X(1:nc,k))
310 0 : do j=1,nc
311 0 : X(j,k) = X(j,k) / tmp
312 0 : X_init(j,k) = X(j,k)
313 : end do
314 0 : X(m,k) = 0d0 ! this will be set later
315 : end do
316 :
317 0 : if (nzlo > 1) then ! combine cells 1 to nzlo
318 0 : mtotal = sum(cell_dm(1:nzlo))
319 0 : do j=1,species ! change to species average of 1:nzlo
320 0 : xa(j,nzlo) = dot_product(cell_dm(1:nzlo),xa(j,1:nzlo))/mtotal
321 : end do
322 0 : do j=1,nc ! change to class average
323 0 : X(j,nzlo) = dot_product(cell_dm(1:nzlo),X(j,1:nzlo))/mtotal
324 : end do
325 0 : cell_dm(nzlo) = mtotal
326 0 : sumx = sum(X(1:nc,nzlo))
327 0 : do j=1,nc
328 0 : X(j,nzlo) = X(j,nzlo)/sumx
329 0 : X_init(j,nzlo) = X(j,nzlo)
330 : end do
331 : end if
332 0 : dm_bar(nzlo+1) = 0.5d0*(cell_dm(nzlo) + cell_dm(nzlo+1))
333 :
334 0 : do i=1,nc ! save for species conservation checks
335 0 : mass_init(i) = dot_product(cell_dm(nzlo:nzhi),X(i,nzlo:nzhi))
336 : end do
337 0 : sum_mass_nzlo_nzhi = sum(cell_dm(nzlo:nzhi))
338 :
339 : ! get info that is held constant during solve
340 : if (dbg) write(*,*) 'call setup_struct_info'
341 : call setup_struct_info( &
342 : s, nz, nzlo, nzhi, species, nc, m, X, A, tiny_X, &
343 : dlnP_dm_face, dlnT_dm_face, dlnRho_dm_face, cell_dm, dm_in, &
344 : abar, free_e, T, lnT, rho, lnd, L_face, r_face, alfa_face, &
345 : class, class_chem_id, calculate_ionization, nsmooth_typical_charge, &
346 : min_T_for_radaccel, max_T_for_radaccel, &
347 : min_Z_for_radaccel, max_Z_for_radaccel, &
348 : screening_for_radaccel, rho_face, T_face, four_pi_r2_rho_face, &
349 : dlnP_dr_face, dlnT_dr_face, dlnRho_dr_face, &
350 : Z, typical_charge, xm_face, &
351 : rad_accel_face, log10_g_rad, g_rad_face, &
352 0 : kmax_rad_accel, ierr)
353 : if (dbg) write(*,*) 'done setup_struct_info'
354 0 : if (failed('setup_struct_info')) then
355 0 : call dealloc
356 0 : return
357 : end if
358 :
359 0 : mtotal = sum(cell_dm(nzlo:nzhi))
360 0 : do j=1,nc
361 : xtotal_init(j) = &
362 0 : dot_product(cell_dm(nzlo:nzhi),X(j,nzlo:nzhi))/mtotal
363 : end do
364 :
365 0 : tol_correction_max = s% diffusion_tol_correction_max
366 0 : tol_correction_norm = s% diffusion_tol_correction_norm
367 0 : dt = 0
368 0 : min_dt = total_time/(1000*max_steps)
369 0 : time = 0
370 0 : steps_used = 0
371 0 : total_num_retries = 0
372 0 : total_num_iters = 0
373 :
374 : if (dbg) write(*,*) '1st time call fixup'
375 : call fixup(s, &
376 : nz, nc, m, nzlo, nzhi, &
377 : total_num_iters, s% diffusion_min_X_hard_limit, X_total_atol, X_total_rtol, &
378 : cell_dm, mtotal, xtotal_init, X, &
379 : lnT, sum_ending_mass, ending_mass, ending_dX_dm, &
380 0 : bad_j, bad_k, bad_X, bad_sum, bad_Xsum, ierr)
381 0 : if (failed('fixup')) then
382 0 : call dealloc
383 0 : return
384 : end if
385 0 : have_changed_matrix_coeffs = .false.
386 :
387 0 : vc_old = 0
388 0 : dt_old = 0
389 0 : avg_del = 0
390 :
391 0 : if (use_isolve) then
392 0 : call solve_with_isolve(ierr)
393 : else
394 0 : call do_step_loop(ierr)
395 : end if
396 0 : if (ierr /= 0) then
397 0 : call dealloc
398 0 : return
399 : end if
400 :
401 0 : if (total_time - time > 1d-10*total_time .or. steps_used == 0) then
402 0 : ierr = -1 ! failed to finish
403 : end if
404 :
405 0 : if (ierr == 0) then
406 :
407 0 : call set_new_xa(nz, nzlo, nzhi, species, nc, m, class, X_init, X, cell_dm, xa)
408 :
409 0 : do k=nzlo,nzhi
410 0 : do j=1,nc
411 0 : X_final(j,k) = X(j,k)
412 : end do
413 : end do
414 :
415 0 : do k=1,nzlo-1 ! fully mixed
416 0 : do i=1,species
417 0 : xa(i,k) = xa(i,nzlo)
418 : end do
419 0 : do j=1,nc
420 0 : X_final(j,k) = X(j,nzlo)
421 : end do
422 : end do
423 :
424 : end if
425 :
426 0 : call dealloc
427 :
428 : if (dbg .and. ierr /= 0) then
429 : call mesa_error(__FILE__,__LINE__,'failed in diffusion')
430 : end if
431 :
432 0 : if (do_timing) then
433 0 : call system_clock(time1,clock_rate)
434 0 : write(*,2) 'diffusion model timing', s% model_number, &
435 0 : dble(time1 - time0)/clock_rate
436 : end if
437 :
438 :
439 : contains
440 :
441 :
442 0 : subroutine solve_with_isolve(ierr)
443 0 : use num_def
444 : use num_lib
445 : use mtx_def
446 : use mtx_lib
447 :
448 : integer, intent(out) :: ierr
449 :
450 : integer :: &
451 : i, k, lout, iout, idid, ijac, max_steps, &
452 : imas, mlmas, mumas, itol, &
453 : nzmax, lrd, lid, isparse, &
454 : liwork, lwork, caller_id, which_solver, which_decsol
455 0 : real(dp) :: h, max_step_size, x_min, x_max
456 :
457 0 : integer, pointer :: iwork(:) !(liwork)
458 0 : real(dp), pointer :: work(:) !(lwork)
459 0 : integer, pointer :: ipar_decsol(:) !(lid)
460 0 : real(dp), pointer :: rpar_decsol(:) !(lrd)
461 :
462 0 : real(dp) :: atol(1), rtol(1), t, tend
463 0 : real(dp), dimension(:), pointer :: lblk, dblk, ublk ! =(nc,nc,n)
464 0 : real(dp), dimension(:), pointer :: uf_lblk, uf_dblk, uf_ublk ! =(nc,nc,n)
465 :
466 : integer, parameter :: lrpar = 1, lipar = 1
467 0 : real(dp), target :: rpar_ary(lrpar)
468 : integer, target :: ipar_ary(lipar)
469 : real(dp), pointer :: rpar(:)
470 : integer, pointer :: ipar(:)
471 :
472 0 : real(dp), pointer :: vars1(:), vars(:,:)
473 :
474 : include 'formats'
475 : ierr = 0
476 :
477 0 : rpar => rpar_ary
478 0 : ipar => ipar_ary
479 :
480 0 : which_decsol = bcyclic_dble
481 :
482 0 : which_solver = solver_option(s% diffusion_isolve_solver, ierr)
483 0 : if (ierr /= 0) then
484 : write(*,*) 'unknown solver name for diffusion_isolve_solver' // &
485 0 : trim(s% diffusion_isolve_solver)
486 : return
487 : end if
488 :
489 0 : allocate(vars1(neqs))
490 : if (which_decsol == bcyclic_dble) then
491 0 : allocate(lblk(nc*nc*n), dblk(nc*nc*n), ublk(nc*nc*n))
492 0 : allocate(uf_lblk(nc*nc*n), uf_dblk(nc*nc*n), uf_ublk(nc*nc*n))
493 : else
494 : nullify(lblk, dblk, ublk, uf_lblk, uf_dblk, uf_ublk)
495 : end if
496 :
497 0 : itol = 0
498 0 : rtol(1) = s% diffusion_rtol_for_isolve
499 0 : atol(1) = s% diffusion_atol_for_isolve
500 :
501 0 : iout = 1
502 0 : max_steps = s% diffusion_maxsteps_for_isolve
503 0 : isparse = 0
504 0 : lout = 6
505 0 : caller_id = 0
506 :
507 0 : x_min = s% diffusion_min_X_hard_limit
508 0 : x_max = 1d0 - s% diffusion_min_X_hard_limit
509 :
510 0 : ijac = 1 ! analytic jacobian
511 :
512 0 : imas = 0
513 0 : mlmas = 0
514 0 : mumas = 0
515 :
516 0 : ipar = 0
517 0 : rpar = 0
518 :
519 : lid = 0; lrd = 0
520 : if (which_decsol == lapack) then
521 : call mesa_error(__FILE__,__LINE__,'not ready for lapack in diff')
522 : nzmax = 0
523 : call lapack_work_sizes(neqs,lrd,lid)
524 0 : else if (which_decsol == bcyclic_dble) then
525 0 : nzmax = 0
526 0 : call bcyclic_dble_work_sizes(nc,n,lrd,lid)
527 : else
528 : write(*,*) 'unknown value for which_decsol', which_decsol
529 : call mesa_error(__FILE__,__LINE__)
530 : end if
531 :
532 0 : call isolve_work_sizes(neqs,nzmax,imas,kl,ku,mlmas,mumas,liwork,lwork)
533 :
534 0 : allocate(iwork(liwork),work(lwork),ipar_decsol(lid),rpar_decsol(lrd),stat=ierr)
535 0 : if (ierr /= 0) then
536 0 : write(*,*) 'allocate ierr', ierr
537 0 : call mesa_error(__FILE__,__LINE__) ! test_int_support
538 : end if
539 :
540 0 : iwork = 0
541 0 : work = 0
542 :
543 0 : t = 0
544 0 : tend = total_time
545 0 : h = total_time ! initial step size
546 0 : max_step_size = total_time
547 :
548 : ! set vars
549 0 : vars(1:nc,1:n) => vars1(1:neqs)
550 0 : do i=1,n
551 0 : k = i+nzlo-1
552 0 : do j=1,nc
553 0 : vars(j,i) = X(j,k)
554 : end do
555 : end do
556 :
557 : if (which_decsol == lapack) then
558 : call mesa_error(__FILE__,__LINE__,'not supported')
559 : call isolve( &
560 : which_solver, neqs, fcn, t, vars1, tend, &
561 : h, max_step_size, max_steps, &
562 : rtol, atol, itol, x_min, x_max, &
563 : jac, ijac, null_sjac, nzmax, isparse, kl, ku, &
564 : null_mas, imas, mlmas, mumas, &
565 : solout, iout, &
566 : lapack_decsol, null_decsols, null_decsolblk, &
567 : lrd, rpar_decsol, lid, ipar_decsol, &
568 : caller_id, 0, 0, &
569 : lblk, dblk, ublk, uf_lblk, uf_dblk, uf_ublk, &
570 : null_fcn_blk_dble, null_jac_blk_dble, &
571 : work, lwork, iwork, liwork, &
572 : lrpar, rpar, lipar, ipar, &
573 : lout, idid)
574 : else if (which_decsol == bcyclic_dble) then
575 : call isolve( &
576 : which_solver, neqs, null_fcn, t, vars1, tend, &
577 : h, max_step_size, max_steps, &
578 : rtol, atol, itol, x_min, x_max, &
579 : null_jac, ijac, null_sjac, nzmax, isparse, kl, ku, &
580 : null_mas, imas, mlmas, mumas, &
581 : solout, iout, &
582 : null_decsol, null_decsols, bcyclic_dble_decsolblk, &
583 : lrd, rpar_decsol, lid, ipar_decsol, &
584 : caller_id, nc, n, &
585 : lblk, dblk, ublk, uf_lblk, uf_dblk, uf_ublk, &
586 : fcn_blk_dble, jac_blk_dble, &
587 : work, lwork, iwork, liwork, &
588 : lrpar, rpar, lipar, ipar, &
589 0 : lout, idid)
590 : else
591 : call mesa_error(__FILE__,__LINE__,'bad which_decsol in mod_diffusion')
592 : end if
593 :
594 0 : if (idid /= 1) ierr = -1
595 : !if (ierr /= 0) then
596 : ! write(*,*) 'solver returned ierr /= 0', idid
597 : ! call mesa_error(__FILE__,__LINE__)
598 : !end if
599 0 : steps_used = iwork(17) ! number of accepted steps
600 0 : total_num_retries = iwork(16) - iwork(17) ! num computed - num accepted
601 0 : time = tend
602 :
603 0 : do i=1,n
604 0 : k = i+nzlo-1
605 0 : do j=1,nc
606 0 : X(j,k) = vars(j,i)
607 : end do
608 : end do
609 :
610 : if (which_decsol == bcyclic_dble) &
611 0 : deallocate(lblk,dblk,ublk,uf_lblk,uf_dblk,uf_ublk)
612 0 : deallocate(vars1,iwork,work,ipar_decsol,rpar_decsol)
613 :
614 :
615 0 : end subroutine solve_with_isolve
616 :
617 :
618 0 : subroutine fcn_blk_dble(&
619 : n_blk_dble,caller_id,nvar_blk_dble,nz_blk_dble,&
620 : time,h,vars,f,lrpar,rpar,lipar,ipar,ierr)
621 0 : use const_def, only: dp
622 : integer, intent(in) :: n_blk_dble, caller_id, &
623 : nvar_blk_dble, nz_blk_dble, lrpar, lipar
624 : real(dp), intent(in) :: time,h
625 : real(dp), intent(inout), pointer :: vars(:) ! (n)
626 : real(dp), intent(inout), pointer :: f(:) ! (n) dvars/dt
627 : integer, intent(inout), pointer :: ipar(:) ! (lipar)
628 : real(dp), intent(inout), pointer :: rpar(:) ! (lrpar)
629 : integer, intent(out) :: ierr ! nonzero means retry with smaller timestep.
630 :
631 : integer :: i, j, k
632 0 : real(dp), pointer :: f2(:,:)
633 :
634 : include 'formats'
635 :
636 0 : ierr = 0
637 :
638 0 : if (n_blk_dble /= neqs .or. nvar_blk_dble /= nc .or. nz_blk_dble /= n) &
639 0 : call mesa_error(__FILE__,__LINE__,'bad args for fcn_blk_dble')
640 :
641 0 : call update_for_new_vars(vars, time, ierr)
642 0 : if (ierr /= 0) then
643 : !call mesa_error(__FILE__,__LINE__,'fcn_blk_dble')
644 : return
645 : end if
646 :
647 : call get_dX_dt( &
648 : nz, nzlo, nzhi, m, nc, X, X_face, C_div_X, GT_face, AD_face, SIG_face, &
649 0 : alfa_face, cell_dm, v_advection_face, upwind_limit, dX_dt)
650 :
651 0 : f2(1:nc,1:n) => f(1:neqs)
652 0 : do i=1,n
653 0 : k = i+nzlo-1
654 0 : do j=1,nc
655 0 : f2(j,i) = dX_dt(j,k)
656 : end do
657 : end do
658 :
659 0 : end subroutine fcn_blk_dble
660 :
661 :
662 0 : subroutine jac_blk_dble( &
663 : n_blk_dble,caller_id,nvar_blk_dble,nz_blk_dble,&
664 : time,h,vars,f,lblk1,dblk1,ublk1,lrpar,rpar,lipar,ipar,ierr)
665 : use const_def,only: dp
666 : integer,intent(in) :: n_blk_dble,caller_id,nvar_blk_dble,nz_blk_dble,lrpar,lipar
667 : real(dp),intent(in) :: time,h
668 : real(dp),intent(inout), pointer :: vars(:) ! (n)
669 : real(dp),intent(inout), pointer :: f(:) ! (n) dvars/dt
670 : real(dp),dimension(:),pointer,intent(inout) :: lblk1,dblk1,ublk1
671 : integer,intent(inout),pointer :: ipar(:) ! (lipar)
672 : real(dp),intent(inout),pointer :: rpar(:) ! (lrpar)
673 : integer,intent(out) :: ierr ! nonzero means terminate integration
674 :
675 : integer :: i, j, jj, k
676 0 : real(dp), pointer :: f2(:,:)
677 0 : real(dp),dimension(:,:,:),pointer :: lblk,dblk,ublk
678 :
679 : include 'formats'
680 :
681 0 : ierr = 0
682 :
683 0 : if (n_blk_dble /= neqs .or. nvar_blk_dble /= nc .or. nz_blk_dble /= n) &
684 0 : call mesa_error(__FILE__,__LINE__,'bad args for jac_blk_dble')
685 :
686 0 : call update_for_new_vars(vars, time, ierr)
687 0 : if (ierr /= 0) then
688 : !call mesa_error(__FILE__,__LINE__,'jac_blk_dble')
689 : return
690 : end if
691 :
692 : call get_eqn_matrix_entries( &
693 : nz, nzlo, nzhi, m, nc, .true., X, X_face, C_div_X, &
694 : GT_face, AD_face, SIG_face, &
695 : alfa_face, cell_dm, v_advection_face, upwind_limit, &
696 0 : tiny_X, 0d0, dX, dX_dt, del, em1, e00, ep1)
697 :
698 0 : f2(1:nc,1:n) => f(1:neqs)
699 0 : lblk(1:nc,1:nc,1:n) => lblk1(1:nc*nc*n)
700 0 : dblk(1:nc,1:nc,1:n) => dblk1(1:nc*nc*n)
701 0 : ublk(1:nc,1:nc,1:n) => ublk1(1:nc*nc*n)
702 0 : do i=1,n
703 0 : k = i+nzlo-1
704 0 : do j=1,nc
705 0 : f2(j,i) = dX_dt(j,k)
706 0 : do jj=1,nc
707 0 : lblk(jj,j,i) = em1(jj,j,k)
708 0 : dblk(jj,j,i) = e00(jj,j,k)
709 0 : ublk(jj,j,i) = ep1(jj,j,k)
710 : end do
711 : end do
712 : end do
713 :
714 0 : end subroutine jac_blk_dble
715 :
716 :
717 0 : subroutine update_for_new_vars(vars, time, ierr)
718 : real(dp),intent(inout), pointer :: vars(:)
719 : real(dp), intent(in) :: time
720 : integer,intent(out) :: ierr
721 :
722 : integer :: i, j, k
723 0 : real(dp), pointer :: vars2(:,:)
724 :
725 :
726 : include 'formats'
727 :
728 0 : ierr = 0
729 0 : vars2(1:nc,1:n) => vars(1:neqs)
730 :
731 : ! copy vars to X
732 0 : do i=1,n
733 0 : k = i+nzlo-1
734 0 : do j=1,nc
735 0 : X(j,k) = vars2(j,i)
736 : end do
737 : end do
738 :
739 : if (.true.) then
740 :
741 : call fixup(s, &
742 : nz, nc, m, nzlo, nzhi, total_num_iters, &
743 : s% diffusion_min_X_hard_limit, X_total_atol, X_total_rtol, &
744 : cell_dm, mtotal, xtotal_init, X, &
745 : lnT, sum_ending_mass, ending_mass, ending_dX_dm, &
746 0 : bad_j, bad_k, bad_X, bad_sum, bad_Xsum, ierr)
747 0 : if (ierr /= 0) then
748 0 : s% retry_message = 'element diffusion failed in fixup'
749 0 : if (s% report_ierr) write(*, *) s% retry_message
750 0 : return
751 : end if
752 0 : if (failed('fixup')) return
753 :
754 : ! copy X to vars
755 0 : do i=1,n
756 0 : k = i+nzlo-1
757 0 : do j=1,nc
758 0 : vars2(j,i) = X(j,k)
759 : end do
760 : end do
761 :
762 : end if
763 :
764 : call get_matrix_coeffs( &
765 : s, nz, nc, m, nzlo, nzhi, class(h1), class(he4), &
766 : pure_Coulomb, dt, v_advection_max, tiny_C, diffusion_factor, &
767 : A, X, Z, rho_face, T_face, four_pi_r2_rho_face, &
768 : xm_face, cell_dm, dm_bar, dlnP_dr_face, dlnT_dr_face, dlnRho_dr_face, &
769 : r_face, r_mid, limit_coeffs_face, alfa_face, &
770 : rad_accel_face, log10_g_rad, g_rad_face, &
771 : min_T_for_radaccel, max_T_for_radaccel, &
772 : X_init, X_face, C, C_div_X, C_div_X_face, &
773 : e_ap, e_at, e_ar, e_ax, E_field_face, &
774 : g_ap, g_at, g_ar, g_ax, g_field_face, &
775 : v_advection_face, v_total_face, vlnP_face, vlnT_face, v_rad_face, &
776 0 : GT_face, D_self_face, AD_face, SIG_face, sigma_lnC, ierr)
777 0 : if (ierr /= 0) then
778 0 : write(*,1) 'failed in get_matrix_coeffs', time
779 0 : return
780 : end if
781 0 : if (failed('get_matrix_coeffs')) return
782 :
783 0 : end subroutine update_for_new_vars
784 :
785 :
786 : subroutine fcn(n, x, h, y, f, lrpar, rpar, lipar, ipar, ierr)
787 : use const_def, only: dp
788 : integer, intent(in) :: n, lrpar, lipar
789 : real(dp), intent(in) :: x, h
790 : real(dp), intent(inout) :: y(:)
791 : real(dp), intent(inout) :: f(:) ! dy/dx
792 : integer, intent(inout), pointer :: ipar(:) ! (lipar)
793 : real(dp), intent(inout), pointer :: rpar(:) ! (lrpar)
794 : integer, intent(out) :: ierr ! nonzero means retry with smaller timestep.
795 : ierr = 0
796 : f = 0
797 : end subroutine fcn
798 :
799 :
800 : subroutine jac(n, x, h, y, f, dfdy, ldfy, lrpar, rpar, lipar, ipar, ierr)
801 : use const_def, only: dp
802 : integer, intent(in) :: n, ldfy, lrpar, lipar
803 : real(dp), intent(in) :: x, h
804 : real(dp), intent(inout) :: y(:)
805 : real(dp), intent(inout) :: f(:) ! dy/dx
806 : real(dp), intent(inout) :: dfdy(:,:)
807 : ! dense: dfdy(i, j) = partial f(i) / partial y(j)
808 : ! banded: dfdy(i-j+mujac+1, j) = partial f(i) / partial y(j)
809 : ! uses rows 1 to mljac+mujac+1 of dfdy.
810 : ! The j-th column of the square matrix is stored in the j-th column of the
811 : ! array dfdy as follows:
812 : ! dfdy(mujac+1+i-j, j) = partial f(i) / partial y(j)
813 : ! for max(1, j-mujac)<=i<=min(N, j+mljac)
814 : integer, intent(inout), pointer :: ipar(:) ! (lipar)
815 : real(dp), intent(inout), pointer :: rpar(:) ! (lrpar)
816 : integer, intent(out) :: ierr ! nonzero means terminate integration
817 : ierr = 0
818 : f = 0
819 : dfdy = 0
820 : end subroutine jac
821 :
822 :
823 0 : subroutine solout(nr, xold, time, n, y, rwork_y, iwork_y, interp_y, lrpar, rpar, lipar, ipar, irtrn)
824 : ! nr is the step number.
825 : ! x is the current x value; xold is the previous x value.
826 : ! y is the current y value.
827 : ! irtrn negative means terminate integration.
828 : ! rwork_y and iwork_y hold info for interp_y
829 : ! note that these are not the same as the rwork and iwork arrays for the solver.
830 : use const_def, only: dp
831 : integer, intent(in) :: nr, n, lrpar, lipar
832 : real(dp), intent(in) :: xold, time
833 : real(dp), intent(inout) :: y(:)
834 : ! y can be modified if necessary to keep it in valid range of possible solutions.
835 : real(dp), intent(inout), target :: rwork_y(*)
836 : integer, intent(inout), target :: iwork_y(*)
837 : integer, intent(inout), pointer :: ipar(:) ! (lipar)
838 : real(dp), intent(inout), pointer :: rpar(:) ! (lrpar)
839 : interface
840 : include 'num_interp_y.dek'
841 : end interface
842 : integer, intent(out) :: irtrn ! < 0 causes solver to return to calling program.
843 : include 'formats'
844 0 : if (s% show_diffusion_substep_info .and. mod(nr,10) == 0) write(*,2) 'nr', nr, time
845 0 : irtrn = 0
846 0 : end subroutine solout
847 :
848 :
849 :
850 0 : subroutine do_step_loop(ierr)
851 : integer, intent(out) :: ierr
852 :
853 : include 'formats'
854 0 : ierr = 0
855 :
856 : step_loop: do while &
857 0 : (total_time - time > 1d-10*total_time .and. &
858 0 : steps_used < max_steps)
859 :
860 0 : steps_used = steps_used + 1
861 :
862 : if (dbg) write(*,*) 'call update_coeffs'
863 0 : call update_coeffs((kmax_rad_accel > 0 .and. steps_used > 1), ierr)
864 0 : if (failed('update_coeffs')) return
865 :
866 0 : remaining_time = total_time - time
867 :
868 : if (dbg) write(*,*) 'call get_timescale'
869 : call get_timescale( &
870 : s, nz, nzlo, nzhi, m, nc, -s% diffusion_min_X_hard_limit*0.5d0, v_advection_face, &
871 : upwind_limit, X, X_face, C_div_X, SIG_face, GT_face, AD_face, cell_dm, r_face, &
872 : steps_used, total_num_iters, dbg, iter_dbg, j_dbg, k_dbg, &
873 0 : j, k, class_chem_id, total_time, timescale)
874 :
875 0 : dt = dt_div_timescale*timescale
876 0 : dt = max(dt, 1d-4*remaining_time)
877 0 : if (min_num_substeps > 0) dt = min(dt, total_time/min_num_substeps)
878 :
879 : if (dbg) write(*,3) 'timescale remaining_time dt', &
880 : steps_used, total_num_iters, timescale, &
881 : remaining_time, dt, dt_div_timescale
882 :
883 0 : if (dt >= remaining_time) then
884 : if (dbg) write(*,3) 'dt >= remaining_time', &
885 : steps_used, total_num_iters, dt/remaining_time
886 0 : dt = remaining_time
887 0 : else if (dt > 0.5d0*remaining_time) then
888 : if (dbg) write(*,3) 'dt >= 0.5d0*remaining_time', &
889 : steps_used, total_num_iters, dt/remaining_time
890 0 : dt = 0.5d0*remaining_time
891 : else
892 : if (dbg) write(*,3) 'dt < 0.5d0*remaining_time', &
893 : steps_used, total_num_iters, dt/remaining_time
894 : end if
895 :
896 0 : last_step = time + dt >= total_time*(1d0 - 1d-10)
897 :
898 : if (dbg) write(*,3) &
899 : 'dt, ../total, ../remaining, remaining/total', &
900 : steps_used, total_num_iters, dt, dt/total_time, &
901 : dt/remaining_time, remaining_time/total_time
902 : if (dbg) write(*,*)
903 :
904 0 : max_retries = s% diffusion_max_retries_per_substep
905 0 : retry_loop: do retry_count = 1, max_retries
906 :
907 0 : if (dt < min_dt) then
908 : if (dbg) write(*,1) 'quit because dt < min_dt: min_dt', min_dt
909 : exit step_loop
910 : end if
911 :
912 0 : if (retry_count == 1) then ! save starting X in X_0
913 0 : if (trace) then
914 0 : write(*,'(A)')
915 0 : write(*,2) '1st try for substep', steps_used
916 : end if
917 0 : do k=nzlo,nzhi
918 0 : do j=1,nc
919 0 : X_0(j,k) = X(j,k)
920 0 : X_1(j,k) = X_0(j,k)
921 0 : dX(j,k) = 0d0
922 : end do
923 : end do
924 : else ! prepare for retry of solve_loop
925 0 : if (trace) then
926 0 : write(*,'(A)')
927 0 : write(*,2) 'retry substep', steps_used
928 : end if
929 0 : dt = dt*dt_retry_factor
930 0 : total_num_retries = total_num_retries+1
931 0 : do k=nzlo,nzhi ! restore X from X_0
932 0 : do j=1,nc
933 0 : X(j,k) = X_0(j,k)
934 0 : X_1(j,k) = X_0(j,k)
935 0 : dX(j,k) = 0d0
936 : end do
937 : end do
938 0 : if (have_changed_matrix_coeffs) then ! must recalculate them
939 : if (dbg) write(*,*) 'call get_matrix_coeffs'
940 : call get_matrix_coeffs( &
941 : s, nz, nc, m, nzlo, nzhi, class(h1), class(he4), &
942 : pure_Coulomb, dt, v_advection_max, tiny_C, diffusion_factor, &
943 : A, X, Z, rho_face, T_face, four_pi_r2_rho_face, &
944 : xm_face, cell_dm, dm_bar, dlnP_dr_face, dlnT_dr_face, dlnRho_dr_face, &
945 : r_face, r_mid, limit_coeffs_face, alfa_face, &
946 : rad_accel_face, log10_g_rad, g_rad_face, &
947 : min_T_for_radaccel, max_T_for_radaccel, &
948 : X_init, X_face, C, C_div_X, C_div_X_face, &
949 : e_ap, e_at, e_ar, e_ax, E_field_face, &
950 : g_ap, g_at, g_ar, g_ax, g_field_face, &
951 : v_advection_face, v_total_face, vlnP_face, vlnT_face, v_rad_face, &
952 0 : GT_face, D_self_face, AD_face, SIG_face, sigma_lnC, ierr)
953 0 : if (failed('get_matrix_coeffs')) return
954 : end if
955 : end if
956 :
957 0 : have_changed_matrix_coeffs = .false.
958 :
959 : if (dbg) write(*,*) 'call BE_solve'
960 0 : solved = BE_solve(1d0, last_step, ierr)
961 0 : if (ierr /= 0) exit step_loop
962 :
963 0 : if (solved) then ! this step is done
964 0 : time = time + dt
965 : cycle step_loop
966 : end if
967 :
968 : end do retry_loop
969 :
970 0 : exit step_loop
971 :
972 : end do step_loop
973 :
974 : end subroutine do_step_loop
975 :
976 :
977 0 : logical function BE_solve(theta, do_smooth, ierr) ! Backwards Euler
978 : real(dp), intent(in) :: theta ! multiplier for dt
979 : logical, intent(in) :: do_smooth
980 : integer, intent(out) :: ierr
981 :
982 : include 'formats'
983 :
984 0 : BE_solve = .false.
985 :
986 0 : solve_loop: do num_iters = 1, max_iters_per_substep
987 :
988 0 : if (total_num_iters >= max_iters_total) then
989 0 : ierr = -1
990 0 : return
991 : end if
992 :
993 0 : total_num_iters = total_num_iters+1
994 :
995 : call get_eqn_matrix_entries( &
996 : nz, nzlo, nzhi, m, nc, .false., X, X_face, C_div_X, &
997 : GT_face, AD_face, SIG_face, &
998 : alfa_face, cell_dm, v_advection_face, upwind_limit, &
999 0 : tiny_X, theta*dt, dX, dX_dt, del, em1, e00, ep1)
1000 :
1001 : call factor_and_solve_matrix( &
1002 : nz, nzlo, nzhi, nc, n, neqs, &
1003 : lblk1, dblk1, ublk1, ipiv1, del1, &
1004 0 : lrd, rpar_decsol, lid, ipar_decsol, ierr)
1005 0 : if (ierr /= 0) then
1006 0 : if (dbg .or. trace) write(*,*) 'retry because failed to solve matrix eqn'
1007 0 : ierr = 0; return
1008 : end if
1009 :
1010 : ! set lambda for positivity
1011 : call positivity( &
1012 : nc, nzlo, nzhi, X_1, -s% diffusion_min_X_hard_limit, &
1013 0 : del, lambda, num_iters)
1014 : if (dbg) write(*,2) 'lambda', num_iters, lambda
1015 :
1016 0 : if (lambda < min_lambda) then
1017 0 : lambda = min_lambda
1018 : if (dbg) write(*,2) 'min lambda', num_iters, lambda
1019 : end if
1020 :
1021 0 : do k=nzlo,nzhi ! apply the correction
1022 0 : do j=1,nc
1023 0 : dX(j,k) = dX(j,k) + lambda*del(j,k)
1024 0 : X_1(j,k) = X_0(j,k) + dX(j,k)
1025 0 : X(j,k) = X_1(j,k)
1026 : end do
1027 : end do
1028 :
1029 : call fixup(s, &
1030 : nz, nc, m, nzlo, nzhi, total_num_iters, &
1031 : s% diffusion_min_X_hard_limit, X_total_atol, X_total_rtol, &
1032 : cell_dm, mtotal, xtotal_init, X, &
1033 : lnT, sum_ending_mass, ending_mass, ending_dX_dm, &
1034 0 : bad_j, bad_k, bad_X, bad_sum, bad_Xsum, ierr)
1035 0 : if (ierr /= 0) then
1036 0 : if (dbg .or. trace) then
1037 : write(*,'(a,6i5,3f16.8,1p,99e16.8)') &
1038 0 : '(fixup failed) retry: step try iter total_iters j k X sum Xsum xm', &
1039 0 : steps_used, retry_count, num_iters, total_num_iters, &
1040 0 : bad_j, bad_k, bad_X, bad_sum, bad_Xsum, xm_face(bad_k)/Msun
1041 : end if
1042 :
1043 : if (dbg .and. bad_sum /= 0) then
1044 : call mesa_error(__FILE__,__LINE__,'bad sum from fixup')
1045 :
1046 : end if
1047 0 : ierr = 0; return
1048 : end if
1049 :
1050 0 : if (lambda == 1d0) then
1051 : ! if magnitude of correction is small enough, consider converged.
1052 0 : max_del = maxval(abs(del(1:nc,nzlo:nzhi)))
1053 0 : avg_del = sum(abs(del(1:nc,nzlo:nzhi)))/(nc*n)
1054 0 : if (max_del <= tol_correction_max .and. avg_del <= tol_correction_norm) then
1055 0 : if (dbg .or. trace) &
1056 0 : write(*,3) 'substep converged: iters max_del avg_del dt/total remain/total', &
1057 0 : steps_used, num_iters, max_del, avg_del, dt/total_time, &
1058 0 : remaining_time/total_time
1059 0 : BE_solve = .true.
1060 0 : return
1061 : end if
1062 0 : if (dbg .or. trace) then
1063 0 : if (max_del > tol_correction_max) &
1064 0 : write(*,2) 'max_del > tol_correction_max', num_iters, max_del, tol_correction_max
1065 0 : if (avg_del > tol_correction_norm) &
1066 0 : write(*,2) 'avg_del > tol_correction_norm', num_iters, avg_del, tol_correction_norm
1067 : end if
1068 : end if
1069 :
1070 0 : if (num_iters == max_iters_per_substep) then
1071 0 : if (dbg .or. trace) then
1072 0 : write(*,3) 'retry because num_iters == max_iters_per_substep', &
1073 0 : num_iters, max_iters_per_substep
1074 : end if
1075 0 : return
1076 : end if
1077 :
1078 : ! prepare for next iteration of solve_loop
1079 : call get_matrix_coeffs( &
1080 : s, nz, nc, m, nzlo, nzhi, class(h1), class(he4), &
1081 : pure_Coulomb, dt, v_advection_max, tiny_C, diffusion_factor, &
1082 : A, X, Z, rho_face, T_face, four_pi_r2_rho_face, &
1083 : xm_face, cell_dm, dm_bar, dlnP_dr_face, dlnT_dr_face, dlnRho_dr_face, &
1084 : r_face, r_mid, limit_coeffs_face, alfa_face, &
1085 : rad_accel_face, log10_g_rad, g_rad_face, &
1086 : min_T_for_radaccel, max_T_for_radaccel, &
1087 : X_init, X_face, C, C_div_X, C_div_X_face, &
1088 : e_ap, e_at, e_ar, e_ax, E_field_face, &
1089 : g_ap, g_at, g_ar, g_ax, g_field_face, &
1090 : v_advection_face, v_total_face, vlnP_face, vlnT_face, v_rad_face, &
1091 0 : GT_face, D_self_face, AD_face, SIG_face, sigma_lnC, ierr)
1092 0 : if (failed('get_matrix_coeffs')) return
1093 0 : have_changed_matrix_coeffs = .true.
1094 :
1095 : end do solve_loop
1096 :
1097 0 : call mesa_error(__FILE__,__LINE__,'diffusion bug -- should not fall through solve_loop')
1098 :
1099 0 : end function BE_solve
1100 :
1101 :
1102 0 : subroutine update_coeffs(update_g_rad, ierr)
1103 : logical, intent(in) :: update_g_rad
1104 : integer, intent(out) :: ierr
1105 :
1106 : include 'formats'
1107 :
1108 0 : ierr = 0
1109 :
1110 0 : if (update_g_rad) then
1111 : call update_rad_accel_face( &
1112 : nzlo, nzhi, nc, m, A, X_init, X, &
1113 : log10_g_rad, g_rad_face, &
1114 0 : rad_accel_face, kmax_rad_accel)
1115 : end if
1116 :
1117 : call get_matrix_coeffs( &
1118 : s, nz, nc, m, nzlo, nzhi, class(h1), class(he4), &
1119 : pure_Coulomb, dt, v_advection_max, tiny_C, diffusion_factor, &
1120 : A, X, Z, rho_face, T_face, four_pi_r2_rho_face, &
1121 : xm_face, cell_dm, dm_bar, dlnP_dr_face, dlnT_dr_face, dlnRho_dr_face, &
1122 : r_face, r_mid, limit_coeffs_face, alfa_face, &
1123 : rad_accel_face, log10_g_rad, g_rad_face, &
1124 : min_T_for_radaccel, max_T_for_radaccel, &
1125 : X_init, X_face, C, C_div_X, C_div_X_face, &
1126 : e_ap, e_at, e_ar, e_ax, E_field_face, &
1127 : g_ap, g_at, g_ar, g_ax, g_field_face, &
1128 : v_advection_face, v_total_face, vlnP_face, vlnT_face, v_rad_face, &
1129 0 : GT_face, D_self_face, AD_face, SIG_face, sigma_lnC, ierr)
1130 :
1131 0 : end subroutine update_coeffs
1132 :
1133 0 : subroutine do_alloc(ierr)
1134 : use mtx_lib, only: bcyclic_dble_work_sizes
1135 : use alloc, only: get_integer_work_array
1136 : integer, intent(out) :: ierr
1137 0 : ierr = 0
1138 : call bcyclic_dble_work_sizes(nc, nz, lrd, lid)
1139 0 : call get_integer_work_array(s, ipiv1, nc*nz, nc*nz_alloc_extra, ierr)
1140 0 : if (ierr /= 0) return
1141 0 : call get_integer_work_array(s, ipar_decsol, lid, 0, ierr)
1142 0 : if (ierr /= 0) return
1143 0 : call do_work_arrays(.true.,ierr)
1144 0 : if (ierr /= 0) return
1145 0 : em1(1:nc,1:nc,1:nz) => lblk1(1:nc*nc*nz)
1146 0 : e00(1:nc,1:nc,1:nz) => dblk1(1:nc*nc*nz)
1147 0 : ep1(1:nc,1:nc,1:nz) => ublk1(1:nc*nc*nz)
1148 0 : rhs(1:nc,1:nz) => rhs1(1:nc*nz)
1149 0 : del(1:nc,1:nz) => del1(1:nc*nz)
1150 0 : end subroutine do_alloc
1151 :
1152 :
1153 0 : subroutine dealloc
1154 0 : use alloc, only: return_integer_work_array
1155 0 : call return_integer_work_array(s, ipiv1)
1156 0 : call return_integer_work_array(s, ipar_decsol)
1157 0 : call do_work_arrays(.false.,ierr_dealloc)
1158 0 : end subroutine dealloc
1159 :
1160 :
1161 0 : subroutine do_work_arrays(alloc_flag, ierr)
1162 0 : use alloc, only: work_array
1163 : logical, intent(in) :: alloc_flag
1164 : integer, intent(out) :: ierr
1165 : logical, parameter :: crit = .false.
1166 : ierr = 0
1167 :
1168 : call work_array(s, alloc_flag, crit, &
1169 0 : del1, nc*nz, nc*nz_alloc_extra, 'mod_diffusion', ierr)
1170 : call work_array(s, alloc_flag, crit, &
1171 0 : rhs1, nc*nz, nc*nz_alloc_extra, 'mod_diffusion', ierr)
1172 0 : if (ierr /= 0) return
1173 : call work_array(s, alloc_flag, crit, &
1174 0 : b1, nc*nz, nc*nz_alloc_extra, 'mod_diffusion', ierr)
1175 0 : if (ierr /= 0) return
1176 :
1177 : call work_array(s, alloc_flag, crit, &
1178 0 : lblk1, nc*nc*nz, nc*nc*nz_alloc_extra, 'mod_diffusion', ierr)
1179 0 : if (ierr /= 0) return
1180 : call work_array(s, alloc_flag, crit, &
1181 0 : dblk1, nc*nc*nz, nc*nc*nz_alloc_extra, 'mod_diffusion', ierr)
1182 0 : if (ierr /= 0) return
1183 : call work_array(s, alloc_flag, crit, &
1184 0 : ublk1, nc*nc*nz, nc*nc*nz_alloc_extra, 'mod_diffusion', ierr)
1185 0 : if (ierr /= 0) return
1186 :
1187 : call work_array(s, alloc_flag, crit, &
1188 0 : rpar_decsol, lrd, 0, 'mod_diffusion', ierr)
1189 0 : if (ierr /= 0) return
1190 :
1191 0 : end subroutine do_work_arrays
1192 :
1193 :
1194 0 : logical function failed(str)
1195 : character (len=*) :: str
1196 0 : failed = .false.
1197 0 : if (ierr == 0) return
1198 0 : if (dbg .or. s% report_ierr) write(*,*) 'failed in ' // trim(str), ierr
1199 0 : call dealloc
1200 0 : failed = .true.
1201 0 : end function failed
1202 :
1203 : end subroutine do_solve_diffusion
1204 :
1205 :
1206 0 : subroutine get_timescale( &
1207 0 : s, nz, nzlo, nzhi, m, nc, eps, v_advection_face, upwind_limit, &
1208 0 : X, X_face, C_div_X, SIG_face, GT_face, AD_face, cell_dm, r_face, &
1209 : steps_used, iter, dbg, iter_dbg, i_dbg, k_dbg, i_t, k_t, &
1210 0 : class_chem_id, total_time, dt)
1211 : type (star_info), pointer :: s
1212 : integer, intent(in) :: nz, nzlo, nzhi, m, nc
1213 : real(dp), intent(in) :: eps, upwind_limit, total_time
1214 : real(dp), intent(in), dimension(:,:) :: X, X_face, C_div_X, v_advection_face
1215 : real(dp), intent(in) :: SIG_face(:,:,:) ! (nc,nc,nz)
1216 : real(dp), intent(in) :: GT_face(:,:) ! (nc,nz)
1217 : real(dp), intent(in) :: AD_face(:) ! (nz)
1218 : real(dp), intent(in) :: cell_dm(:) ! (nz)
1219 : real(dp), intent(in) :: r_face(:) ! (nz)
1220 : integer, intent(in) :: steps_used, iter, class_chem_id(:)
1221 : logical, intent(in) :: dbg
1222 : integer, intent(in) :: iter_dbg, i_dbg, k_dbg
1223 : integer, intent(out) :: i_t, k_t
1224 : real(dp), intent(out) :: dt
1225 : integer :: i, j, k
1226 0 : real(dp) :: f, flow_in, flow_out, dt1, dxdt, coeff, &
1227 0 : flow_in_GT, flow_out_GT, flow_in_SIG, flow_out_SIG, &
1228 0 : flow_in_max, flow_out_max, &
1229 0 : flow_in_GT_max, flow_out_GT_max, &
1230 0 : flow_in_SIG_max, flow_out_SIG_max
1231 :
1232 : include 'formats'
1233 :
1234 0 : dt = 1d99
1235 0 : flow_in_max = 0
1236 0 : flow_out_max = 0
1237 0 : flow_in_GT_max = 0
1238 0 : flow_out_GT_max = 0
1239 0 : flow_in_SIG_max = 0
1240 0 : flow_out_SIG_max = 0
1241 :
1242 0 : do k=nzlo+1,nzhi-1
1243 :
1244 0 : do i=1,nc
1245 :
1246 0 : if (X(i,k) < eps) cycle
1247 :
1248 0 : if (abs(v_advection_face(i,k)) >= upwind_limit) then
1249 0 : if (GT_face(i,k) > 0) then
1250 0 : flow_out = GT_face(i,k)*X(i,k)
1251 : else
1252 0 : flow_out = GT_face(i,k)*X(i,k-1)
1253 : end if
1254 : else
1255 0 : flow_out = GT_face(i,k)*X_face(i,k)
1256 : end if
1257 :
1258 0 : if (abs(v_advection_face(i,k+1)) >= upwind_limit) then
1259 0 : if (GT_face(i,k+1) > 0) then
1260 0 : flow_in = GT_face(i,k+1)*X(i,k+1)
1261 : else
1262 0 : flow_in = GT_face(i,k+1)*X(i,k)
1263 : end if
1264 : else
1265 0 : flow_in = GT_face(i,k+1)*X_face(i,k+1)
1266 : end if
1267 :
1268 0 : flow_in_GT = flow_in
1269 0 : flow_out_GT = flow_out
1270 0 : flow_in_SIG = 0
1271 0 : flow_out_SIG = 0
1272 :
1273 0 : do j=1,nc
1274 0 : coeff = SIG_face(i,j,k+1)*X_face(i,k+1)/max(smallX,X_face(j,k+1))
1275 0 : if (j == i) coeff = coeff + AD_face(k+1)
1276 0 : f = coeff*(C_div_X(j,k+1)*X(j,k+1) - C_div_X(j,k)*X(j,k))
1277 0 : flow_in = flow_in + f
1278 0 : flow_in_SIG = flow_in_SIG + f
1279 0 : coeff = SIG_face(i,j,k)*X_face(i,k)/max(smallX,X_face(j,k))
1280 0 : if (j == i) coeff = coeff + AD_face(k)
1281 0 : f = coeff*(C_div_X(j,k)*X(j,k) - C_div_X(j,k-1)*X(j,k-1))
1282 0 : flow_out = flow_out + f
1283 0 : flow_out_SIG = flow_out_SIG + f
1284 : end do
1285 :
1286 0 : if (dbg .and. i == i_dbg .and. k == k_dbg .and. iter == iter_dbg) then
1287 0 : write(*,'(A)')
1288 0 : write(*,*) 'dgb get_timescale'
1289 0 : call show_stuff
1290 0 : write(*,'(A)')
1291 : end if
1292 :
1293 0 : dxdt = (flow_in - flow_out)/cell_dm(k)
1294 0 : if (dxdt >= 0d0) cycle ! only consider decreasing abundances
1295 0 : dt1 = (X(i,k) + eps)/max(1d-50,-dxdt)
1296 :
1297 0 : if (dt1 < dt) then
1298 0 : dt = dt1
1299 0 : k_t = k
1300 0 : i_t = i
1301 0 : flow_in_max = flow_in
1302 0 : flow_out_max = flow_out
1303 0 : flow_in_GT_max = flow_in_GT
1304 0 : flow_out_GT_max = flow_out_GT
1305 0 : flow_in_SIG_max = flow_in_SIG
1306 0 : flow_out_SIG_max = flow_out_SIG
1307 : end if
1308 :
1309 : end do
1310 :
1311 : end do
1312 :
1313 0 : if (dbg) then
1314 :
1315 0 : i = i_t
1316 0 : k = k_t
1317 0 : flow_in = flow_in_max
1318 0 : flow_out = flow_out_max
1319 0 : flow_in_GT = flow_in_GT_max
1320 0 : flow_out_GT = flow_out_GT_max
1321 0 : flow_in_SIG = flow_in_SIG_max
1322 0 : flow_out_SIG = flow_out_SIG_max
1323 :
1324 0 : call show_stuff
1325 :
1326 0 : call mesa_error(__FILE__,__LINE__,'get_timescale')
1327 :
1328 : end if
1329 :
1330 :
1331 : contains
1332 :
1333 :
1334 0 : subroutine show_stuff
1335 : use chem_def, only: chem_isos
1336 :
1337 0 : real(dp) :: dr
1338 :
1339 : include 'formats'
1340 0 : dr = 0.5d0*(r_face(k-1) - r_face(k+1))
1341 :
1342 0 : write(*,'(A)')
1343 0 : write(*,'(A)')
1344 0 : write(*,'(a30,4i6,99e14.5)') 'timescale total_time', &
1345 0 : steps_used, iter, i, k, dt, total_time
1346 : write(*,'(a30,4i6)') trim(chem_isos% name(class_chem_id(i))) // &
1347 0 : ' nzlo, nzhi, nz', nzlo, nzhi, nz
1348 0 : write(*,'(A)')
1349 0 : write(*,'(a30,4i6,99e14.5)') 'dX', &
1350 0 : steps_used, iter, i, k, &
1351 0 : dt*(flow_in - flow_out)/cell_dm(k)
1352 0 : write(*,'(a30,4i6,99e14.5)') 'dX GT', &
1353 0 : steps_used, iter, i, k, &
1354 0 : dt*(flow_in_GT - flow_out_GT)/cell_dm(k)
1355 0 : write(*,'(a30,4i6,99e14.5)') 'dX SIG', &
1356 0 : steps_used, iter, i, k, &
1357 0 : dt*(flow_in_SIG - flow_out_SIG)/cell_dm(k)
1358 0 : write(*,'(A)')
1359 :
1360 0 : write(*,'(a30,4i6,99e14.5)') 'dX in', &
1361 0 : steps_used, iter, i, k, &
1362 0 : dt*flow_in/cell_dm(k)
1363 0 : write(*,'(a30,4i6,99e14.5)') 'dX out', &
1364 0 : steps_used, iter, i, k, &
1365 0 : dt*flow_out/cell_dm(k)
1366 0 : write(*,'(A)')
1367 0 : write(*,'(a30,4i6,99e14.5)') 'rho', &
1368 0 : steps_used, iter, i, k+1, s% rho(k+1)
1369 0 : write(*,'(a30,4i6,99e14.5)') 'rho', &
1370 0 : steps_used, iter, i, k, s% rho(k)
1371 0 : write(*,'(a30,4i6,99e14.5)') 'rho', &
1372 0 : steps_used, iter, i, k-1, s% rho(k-1)
1373 0 : write(*,'(A)')
1374 0 : write(*,'(a30,4i6,99e14.5)') 'C', &
1375 0 : steps_used, iter, i, k+1, X(i,k+1)*C_div_X(i,k+1)
1376 0 : write(*,'(a30,4i6,99e14.5)') 'C', &
1377 0 : steps_used, iter, i, k, X(i,k)*C_div_X(i,k)
1378 0 : write(*,'(a30,4i6,99e14.5)') 'C', &
1379 0 : steps_used, iter, i, k-1, X(i,k-1)*C_div_X(i,k-1)
1380 0 : write(*,'(A)')
1381 0 : write(*,'(a30,4i6,99e14.5)') 'X', &
1382 0 : steps_used, iter, i, k+1, X(i,k+1)
1383 0 : write(*,'(a30,4i6,99e14.5)') 'X', &
1384 0 : steps_used, iter, i, k, X(i,k)
1385 0 : write(*,'(a30,4i6,99e14.5)') 'X', &
1386 0 : steps_used, iter, i, k-1, X(i,k-1)
1387 0 : write(*,'(A)')
1388 0 : write(*,'(a30,4i6,99e14.5)') 'X_face', &
1389 0 : steps_used, iter, i, k+1, X_face(i,k+1)
1390 0 : write(*,'(a30,4i6,99e14.5)') 'X_face', &
1391 0 : steps_used, iter, i, k, X_face(i,k)
1392 0 : write(*,'(A)')
1393 0 : write(*,'(a30,4i6,99e14.5)') 'dt*v_advection_face/dr', &
1394 0 : steps_used, iter, i, k+1, dt*v_advection_face(i,k+1)/dr
1395 0 : write(*,'(a30,4i6,99e14.5)') 'dt*v_advection_face/dr', &
1396 0 : steps_used, iter, i, k, dt*v_advection_face(i,k)/dr
1397 0 : write(*,'(A)')
1398 0 : write(*,'(a30,4i6,99e14.5)') 'dt*v_advection_face', &
1399 0 : steps_used, iter, i, k+1, dt*v_advection_face(i,k+1)
1400 0 : write(*,'(a30,4i6,99e14.5)') 'dt*v_advection_face', &
1401 0 : steps_used, iter, i, k, dt*v_advection_face(i,k)
1402 0 : write(*,'(A)')
1403 0 : write(*,'(a30,4i6,99e14.5)') 'v_advection_face', &
1404 0 : steps_used, iter, i, k+1, v_advection_face(i,k+1)
1405 0 : write(*,'(a30,4i6,99e14.5)') 'v_advection_face', &
1406 0 : steps_used, iter, i, k, v_advection_face(i,k)
1407 0 : write(*,'(A)')
1408 0 : write(*,'(a30,4i6,99e14.5)') 'GT_face', &
1409 0 : steps_used, iter, i, k+1, GT_face(i,k+1)
1410 0 : write(*,'(a30,4i6,99e14.5)') 'GT_face', &
1411 0 : steps_used, iter, i, k, GT_face(i,k)
1412 0 : write(*,'(A)')
1413 :
1414 0 : write(*,2) 'i_dbg', i_dbg
1415 0 : write(*,2) 'k_dbg', k_dbg
1416 0 : write(*,2) 'iter_dbg', iter_dbg
1417 0 : write(*,*) 'i == i_dbg', i == i_dbg
1418 0 : write(*,*) 'k == k_dbg', k == k_dbg
1419 0 : write(*,*) 'iter == iter_dbg', iter == iter_dbg
1420 :
1421 0 : end subroutine show_stuff
1422 :
1423 :
1424 : end subroutine get_timescale
1425 :
1426 :
1427 0 : subroutine get_dX_dt( &
1428 0 : nz, nzlo, nzhi, m, nc, X, X_face, C_div_X, GT_face, AD_face, SIG_face, &
1429 0 : alfa_face, cell_dm, v_advection_face, upwind_limit, dX_dt)
1430 : integer, intent(in) :: nz, nzlo, nzhi, m, nc
1431 : real(dp), dimension(:,:), intent(in) :: X, X_face, C_div_X ! (nc,nz)
1432 : real(dp), intent(in) :: upwind_limit
1433 : real(dp), intent(in), dimension(:) :: alfa_face, cell_dm, AD_face
1434 : real(dp), intent(in), dimension(:,:) :: GT_face, v_advection_face
1435 : real(dp), intent(in), dimension(:,:,:) :: SIG_face
1436 : real(dp), intent(inout), dimension(:,:) :: dX_dt ! (nc,nz)
1437 : integer :: k
1438 0 : !$OMP PARALLEL DO PRIVATE(k) SCHEDULE(dynamic,2)
1439 : do k=nzlo,nzhi
1440 : call get1_dX_dt( &
1441 : k, nz, nzlo, nzhi, m, nc, X, X_face, C_div_X, &
1442 : GT_face, AD_face, SIG_face, alfa_face, cell_dm, v_advection_face, &
1443 : upwind_limit, dX_dt(:,k))
1444 : end do
1445 : !$OMP END PARALLEL DO
1446 0 : end subroutine get_dX_dt
1447 :
1448 :
1449 0 : subroutine get1_dX_dt( &
1450 0 : k, nz, nzlo, nzhi, m, nc, X, X_face, C_div_X, &
1451 0 : GT_face, AD_face, SIG_face, alfa_face, cell_dm, v_advection_face, &
1452 0 : upwind_limit, dX_dt)
1453 :
1454 : integer, intent(in) :: k, nz, nzlo, nzhi, m, nc
1455 : real(dp), dimension(:,:), intent(in) :: X, X_face, C_div_X ! (nc,nz)
1456 : real(dp), intent(in) :: upwind_limit
1457 : real(dp), intent(in), dimension(:) :: alfa_face, cell_dm, AD_face
1458 : real(dp), intent(in), dimension(:,:) :: GT_face, v_advection_face
1459 : real(dp), intent(in), dimension(:,:,:) :: SIG_face
1460 : real(dp), intent(inout), dimension(:) :: dX_dt ! (nc)
1461 : integer :: i, j
1462 0 : real(dp) :: alfa, beta, c, coeff, dC
1463 : include 'formats'
1464 :
1465 0 : dX_dt(:) = 0
1466 :
1467 0 : if (k > nzlo) then ! do face k
1468 :
1469 0 : alfa = alfa_face(k)
1470 0 : beta = 1d0 - alfa
1471 0 : do i=1,nc
1472 :
1473 0 : c = GT_face(i,k)
1474 0 : if (abs(v_advection_face(i,k)) >= upwind_limit) then
1475 0 : if (c > 0) then
1476 0 : dX_dt(i) = dX_dt(i) - X(i,k)*c
1477 : else
1478 0 : dX_dt(i) = dX_dt(i) - X(i,k-1)*c
1479 : end if
1480 : else
1481 0 : dX_dt(i) = dX_dt(i) - X_face(i,k)*c
1482 : end if
1483 :
1484 0 : do j=1,nc
1485 :
1486 0 : c = SIG_face(i,j,k)/max(smallX,X_face(j,k))
1487 0 : coeff = X_face(i,k)*c
1488 0 : if (j == i) coeff = coeff + AD_face(k)
1489 0 : dC = C_div_X(j,k)*X(j,k) - C_div_X(j,k-1)*X(j,k-1)
1490 0 : dX_dt(i) = dX_dt(i) - coeff*dC
1491 :
1492 : end do
1493 :
1494 : end do
1495 :
1496 : end if
1497 :
1498 0 : if (k < nzhi) then ! do face k+1
1499 :
1500 0 : alfa = alfa_face(k+1)
1501 0 : beta = 1d0 - alfa
1502 :
1503 0 : do i=1,nc
1504 :
1505 0 : c = GT_face(i,k+1)
1506 0 : if (abs(v_advection_face(i,k+1)) >= upwind_limit) then
1507 0 : if (c > 0) then
1508 0 : dX_dt(i) = dX_dt(i) + X(i,k+1)*c
1509 : else
1510 0 : dX_dt(i) = dX_dt(i) + X(i,k)*c
1511 : end if
1512 : else
1513 0 : dX_dt(i) = dX_dt(i) + X_face(i,k+1)*c
1514 : end if
1515 :
1516 0 : do j=1,nc
1517 :
1518 0 : c = SIG_face(i,j,k+1)/max(smallX,X_face(j,k+1))
1519 0 : coeff = X_face(i,k+1)*c
1520 0 : if (j == i) coeff = coeff + AD_face(k+1)
1521 0 : dC = C_div_X(j,k+1)*X(j,k+1) - C_div_X(j,k)*X(j,k)
1522 0 : dX_dt(i) = dX_dt(i) + coeff*dC
1523 :
1524 : end do
1525 :
1526 : end do
1527 :
1528 : end if
1529 :
1530 0 : do i=1,nc
1531 0 : dX_dt(i) = dX_dt(i)/cell_dm(k)
1532 : end do
1533 :
1534 0 : end subroutine get1_dX_dt
1535 :
1536 :
1537 0 : subroutine get_eqn_matrix_entries( &
1538 0 : nz, nzlo, nzhi, m, nc, jac_only, X, X_face, C_div_X, GT_face, AD_face, SIG_face, &
1539 0 : alfa_face, cell_dm, v_advection_face, upwind_limit, &
1540 0 : tiny_X, dt, dX, dX_dt, rhs, em1, e00, ep1)
1541 : integer, intent(in) :: nz, nzlo, nzhi, m, nc
1542 : logical, intent(in) :: jac_only
1543 : real(dp), dimension(:,:), intent(in) :: X, X_face, C_div_X ! (nc,nz)
1544 : real(dp), intent(in) :: upwind_limit, tiny_X, dt
1545 : real(dp), intent(in), dimension(:) :: alfa_face, cell_dm, AD_face
1546 : real(dp), intent(in), dimension(:,:) :: GT_face, v_advection_face
1547 : real(dp), intent(in), dimension(:,:,:) :: SIG_face
1548 : real(dp), intent(in), dimension(:,:) :: dX ! (nc,nz)
1549 : real(dp), intent(inout), dimension(:,:) :: dX_dt ! (nc,nz)
1550 : real(dp), intent(inout), dimension(:,:) :: rhs ! (nc,nz)
1551 : real(dp), intent(inout), dimension(:,:,:) :: em1, e00, ep1 ! (nc,nc,nz)
1552 : integer :: k
1553 : ! lhs(i,k) := X(i,k) - (flow(i,k) - flow(i,k-1))*dt/cell_dm(k)
1554 : ! rhs(i,k) := X_prev(i,k)
1555 : ! em1(i,j,k) = d(lhs(i,k))/d(X(j,k-1))
1556 : ! e00(i,j,k) = d(lhs(i,k))/d(X(j,k))
1557 : ! ep1(i,j,k) = d(lhs(i,k))/d(X(j,k+1))
1558 : include 'formats'
1559 0 : !$OMP PARALLEL DO PRIVATE(k) SCHEDULE(dynamic,2)
1560 : do k=nzlo,nzhi
1561 : call get1_eqn_matrix_entries( &
1562 : k, nz, nzlo, nzhi, m, nc, jac_only, X, X_face, C_div_X, &
1563 : GT_face, AD_face, SIG_face, alfa_face, cell_dm, v_advection_face, &
1564 : upwind_limit, tiny_X, dt, dX, dX_dt(:,k), rhs, em1, e00, ep1)
1565 : end do
1566 : !$OMP END PARALLEL DO
1567 0 : end subroutine get_eqn_matrix_entries
1568 :
1569 :
1570 0 : subroutine get1_eqn_matrix_entries( &
1571 0 : k, nz, nzlo, nzhi, m, nc, jac_only, X, X_face, C_div_X, &
1572 0 : GT_face, AD_face, SIG_face, alfa_face, cell_dm, v_advection_face, &
1573 0 : upwind_limit, tiny_X, dt, dX, dX_dt, rhs, em1, e00, ep1)
1574 :
1575 : integer, intent(in) :: k, nz, nzlo, nzhi, m, nc
1576 : logical, intent(in) :: jac_only
1577 : real(dp), dimension(:,:), intent(in) :: X, X_face, C_div_X ! (nc,nz)
1578 : real(dp), intent(in) :: upwind_limit, tiny_X, dt
1579 : real(dp), intent(in), dimension(:) :: alfa_face, cell_dm, AD_face
1580 : real(dp), intent(in), dimension(:,:) :: GT_face, v_advection_face
1581 : real(dp), intent(in), dimension(:,:,:) :: SIG_face
1582 : real(dp), intent(in), dimension(:,:) :: dX ! (nc,nz)
1583 : real(dp), intent(inout), dimension(:) :: dX_dt ! (nc)
1584 : real(dp), intent(inout), dimension(:,:) :: rhs ! (nc,nz)
1585 : real(dp), intent(inout), dimension(:,:,:) :: em1, e00, ep1 ! (nc,nc,nz)
1586 : integer :: i, j
1587 0 : real(dp) :: alfa, beta, c, coeff, dC_dXj00, dC_dXjm1, dC_dXjp1, &
1588 0 : dC, dcoeff_dXjm1, dcoeff_dXj00, dcoeff_dXjp1, dt_div_dm
1589 :
1590 : include 'formats'
1591 :
1592 0 : dX_dt(:) = 0; em1(:,:,k) = 0; e00(:,:,k) = 0; ep1(:,:,k) = 0
1593 :
1594 0 : if (k > nzlo) then ! do face k
1595 :
1596 0 : alfa = alfa_face(k)
1597 0 : beta = 1d0 - alfa
1598 0 : do i=1,nc
1599 :
1600 0 : c = GT_face(i,k)
1601 0 : if (abs(v_advection_face(i,k)) >= upwind_limit) then
1602 0 : if (c > 0) then
1603 0 : dX_dt(i) = dX_dt(i) - X(i,k)*c
1604 0 : e00(i,i,k) = e00(i,i,k) - c
1605 : else
1606 0 : dX_dt(i) = dX_dt(i) - X(i,k-1)*c
1607 0 : em1(i,i,k) = em1(i,i,k) - c
1608 : end if
1609 : else
1610 0 : dX_dt(i) = dX_dt(i) - X_face(i,k)*c
1611 0 : em1(i,i,k) = em1(i,i,k) - beta*c
1612 0 : e00(i,i,k) = e00(i,i,k) - alfa*c
1613 : end if
1614 :
1615 0 : do j=1,nc
1616 :
1617 0 : c = SIG_face(i,j,k)/max(smallX,X_face(j,k))
1618 0 : coeff = X_face(i,k)*c
1619 0 : if (j == i) coeff = coeff + AD_face(k)
1620 0 : dC = C_div_X(j,k)*X(j,k) - C_div_X(j,k-1)*X(j,k-1)
1621 0 : dX_dt(i) = dX_dt(i) - coeff*dC
1622 :
1623 0 : if (j == i) then
1624 0 : em1(i,j,k) = em1(i,j,k) - beta*c*dC
1625 0 : e00(i,j,k) = e00(i,j,k) - alfa*c*dC
1626 : end if
1627 :
1628 0 : dC_dXjm1 = -C_div_X(j,k-1)
1629 0 : dC_dXj00 = C_div_X(j,k)
1630 0 : em1(i,j,k) = em1(i,j,k) - coeff*dC_dXjm1
1631 0 : e00(i,j,k) = e00(i,j,k) - coeff*dC_dXj00
1632 :
1633 0 : if (use_dcoeff_dX .and. X_face(j,k) > smallX) then
1634 0 : dcoeff_dXjm1 = -beta*coeff/X_face(j,k)
1635 0 : dcoeff_dXj00 = -alfa*coeff/X_face(j,k)
1636 0 : em1(i,j,k) = em1(i,j,k) - dcoeff_dXjm1*dC
1637 0 : e00(i,j,k) = e00(i,j,k) - dcoeff_dXj00*dC
1638 : end if
1639 :
1640 : end do
1641 :
1642 : end do
1643 :
1644 : end if
1645 :
1646 0 : if (k < nzhi) then ! do face k+1
1647 :
1648 0 : alfa = alfa_face(k+1)
1649 0 : beta = 1d0 - alfa
1650 :
1651 0 : do i=1,nc
1652 :
1653 0 : c = GT_face(i,k+1)
1654 0 : if (abs(v_advection_face(i,k+1)) >= upwind_limit) then
1655 0 : if (c > 0) then
1656 0 : dX_dt(i) = dX_dt(i) + X(i,k+1)*c
1657 0 : ep1(i,i,k) = ep1(i,i,k) + c
1658 : else
1659 0 : dX_dt(i) = dX_dt(i) + X(i,k)*c
1660 0 : e00(i,i,k) = e00(i,i,k) + c
1661 : end if
1662 : else
1663 0 : dX_dt(i) = dX_dt(i) + X_face(i,k+1)*c
1664 0 : e00(i,i,k) = e00(i,i,k) + beta*c
1665 0 : ep1(i,i,k) = ep1(i,i,k) + alfa*c
1666 : end if
1667 :
1668 0 : do j=1,nc
1669 :
1670 0 : c = SIG_face(i,j,k+1)/max(smallX,X_face(j,k+1))
1671 0 : coeff = X_face(i,k+1)*c
1672 0 : if (j == i) coeff = coeff + AD_face(k+1)
1673 0 : dC = C_div_X(j,k+1)*X(j,k+1) - C_div_X(j,k)*X(j,k)
1674 0 : dX_dt(i) = dX_dt(i) + coeff*dC
1675 :
1676 0 : if (j == i) then
1677 0 : e00(i,j,k) = e00(i,j,k) + beta*c*dC
1678 0 : ep1(i,j,k) = ep1(i,j,k) + alfa*c*dC
1679 : end if
1680 :
1681 0 : dC_dXj00 = -C_div_X(j,k)
1682 0 : dC_dXjp1 = C_div_X(j,k+1)
1683 0 : e00(i,j,k) = e00(i,j,k) + coeff*dC_dXj00
1684 0 : ep1(i,j,k) = ep1(i,j,k) + coeff*dC_dXjp1
1685 :
1686 0 : if (use_dcoeff_dX .and. X_face(j,k+1) > smallX) then
1687 0 : dcoeff_dXj00 = -beta*coeff/X_face(j,k+1)
1688 0 : dcoeff_dXjp1 = -alfa*coeff/X_face(j,k+1)
1689 0 : e00(i,j,k) = e00(i,j,k) + dcoeff_dXj00*dC
1690 0 : ep1(i,j,k) = ep1(i,j,k) + dcoeff_dXjp1*dC
1691 : end if
1692 :
1693 : end do
1694 :
1695 : end do
1696 :
1697 : end if
1698 :
1699 0 : if (jac_only) then
1700 0 : do j=1,nc
1701 0 : dX_dt(j) = dX_dt(j)/cell_dm(k)
1702 0 : do i=1,nc
1703 0 : em1(i,j,k) = em1(i,j,k)/cell_dm(k)
1704 0 : e00(i,j,k) = e00(i,j,k)/cell_dm(k)
1705 0 : ep1(i,j,k) = ep1(i,j,k)/cell_dm(k)
1706 : end do
1707 : end do
1708 : return
1709 : end if
1710 :
1711 0 : dt_div_dm = dt/cell_dm(k)
1712 0 : do j=1,nc
1713 0 : dX_dt(j) = dX_dt(j)/cell_dm(k)
1714 0 : rhs(j,k) = dX_dt(j)*dt - dX(j,k)
1715 0 : do i=1,nc
1716 0 : em1(i,j,k) = -em1(i,j,k)*dt_div_dm
1717 0 : e00(i,j,k) = -e00(i,j,k)*dt_div_dm
1718 0 : ep1(i,j,k) = -ep1(i,j,k)*dt_div_dm
1719 : end do
1720 0 : e00(j,j,k) = e00(j,j,k) + 1d0
1721 : end do
1722 :
1723 : end subroutine get1_eqn_matrix_entries
1724 :
1725 :
1726 0 : subroutine factor_and_solve_matrix( &
1727 : nz, nzlo, nzhi, nc, n, neqs, &
1728 : lblk1_nz, dblk1_nz, ublk1_nz, ipiv1_nz, del1_nz, &
1729 : lrd, rpar_decsol, lid, ipar_decsol, ierr)
1730 :
1731 : use mtx_lib, only: bcyclic_dble_decsolblk, block_dble_mv
1732 : integer, intent(in) :: nz, nzlo, nzhi, nc, n, neqs
1733 : real(dp), pointer, dimension(:) :: &
1734 : lblk1_nz, dblk1_nz, ublk1_nz, del1_nz
1735 : integer, pointer :: ipiv1_nz(:)
1736 :
1737 : integer :: lrd, lid
1738 : integer, pointer :: ipar_decsol(:)
1739 : real(dp), pointer :: rpar_decsol(:)
1740 : integer, intent(out) :: ierr
1741 :
1742 : integer :: caller_id, ierr2
1743 : integer, pointer :: ipiv1_n(:)
1744 : real(dp), pointer, dimension(:) :: rhs1_n, lblk1_n, dblk1_n, ublk1_n
1745 0 : real(dp), pointer, dimension(:,:,:) :: lblk, dblk, ublk
1746 :
1747 : ierr2 = 0
1748 0 : caller_id = 0
1749 :
1750 0 : rhs1_n(1:nc*n) => del1_nz(1+nc*(nzlo-1):nc*nzhi)
1751 0 : ipiv1_n(1:nc*n) => ipiv1_nz(1+nc*(nzlo-1):nc*nzhi)
1752 :
1753 0 : lblk1_n(1:nc*nc*n) => lblk1_nz(1+nc*nc*(nzlo-1):nc*nc*nzhi)
1754 0 : dblk1_n(1:nc*nc*n) => dblk1_nz(1+nc*nc*(nzlo-1):nc*nc*nzhi)
1755 0 : ublk1_n(1:nc*nc*n) => ublk1_nz(1+nc*nc*(nzlo-1):nc*nc*nzhi)
1756 :
1757 0 : lblk(1:nc,1:nc,1:n) => lblk1_n(1:nc*nc*n)
1758 0 : dblk(1:nc,1:nc,1:n) => dblk1_n(1:nc*nc*n)
1759 0 : ublk(1:nc,1:nc,1:n) => ublk1_n(1:nc*nc*n)
1760 :
1761 :
1762 : ! factor
1763 : call bcyclic_dble_decsolblk( &
1764 : 0, caller_id, nc, n, lblk1_n, dblk1_n, ublk1_n, rhs1_n, ipiv1_n, &
1765 0 : lrd, rpar_decsol, lid, ipar_decsol, ierr)
1766 0 : if (ierr /= 0) then
1767 : call bcyclic_dble_decsolblk( &
1768 : 2, caller_id, nc, n, lblk1_n, dblk1_n, ublk1_n, rhs1_n, ipiv1_n, &
1769 0 : lrd, rpar_decsol, lid, ipar_decsol, ierr2)
1770 0 : return
1771 : end if
1772 :
1773 : ! solve
1774 : call bcyclic_dble_decsolblk( &
1775 : 1, caller_id, nc, n, lblk1_n, dblk1_n, ublk1_n, rhs1_n, ipiv1_n, &
1776 0 : lrd, rpar_decsol, lid, ipar_decsol, ierr)
1777 : ! deallocate
1778 : call bcyclic_dble_decsolblk( &
1779 : 2, caller_id, nc, n, lblk1_n, dblk1_n, ublk1_n, rhs1_n, ipiv1_n, &
1780 0 : lrd, rpar_decsol, lid, ipar_decsol, ierr2)
1781 :
1782 0 : end subroutine factor_and_solve_matrix
1783 :
1784 :
1785 0 : subroutine positivity(nc, nzlo, nzhi, X, eps, del, lambda, num_iters)
1786 : integer, intent(in) :: nc, nzlo, nzhi, num_iters
1787 : real(dp), intent(in) :: eps
1788 : real(dp), intent(in), dimension(:,:) :: X, del
1789 : real(dp), intent(out) :: lambda
1790 :
1791 : integer :: j, k, bad_j, bad_k
1792 0 : real(dp) :: alpha, new_x, old_x, dx
1793 :
1794 : include 'formats'
1795 :
1796 0 : lambda = 1d0
1797 0 : bad_j = 0
1798 0 : do k=nzlo,nzhi
1799 0 : do j=1,nc
1800 0 : old_x = X(j,k)
1801 0 : if (old_x < 1d-99) cycle
1802 0 : dx = del(j,k)
1803 0 : new_x = old_x + dx
1804 0 : if (new_x >= 0) cycle
1805 0 : alpha = -(old_x + eps)/dx
1806 : ! so dx*alpha = -old_x - eps,
1807 : ! and therefore old_x + alpha*dx = -eps
1808 0 : if (alpha < lambda) then
1809 0 : lambda = alpha
1810 0 : bad_k = k
1811 0 : bad_j = j
1812 : end if
1813 : end do
1814 : end do
1815 0 : if (lambda < 1) lambda = 0.8d0*lambda ! under correct
1816 :
1817 0 : end subroutine positivity
1818 :
1819 :
1820 : subroutine write_plot_data( &
1821 : nzlo, nzhi, nc, class_chem_id, &
1822 : X1, X2, v)
1823 : use chem_def, only: chem_isos
1824 : use utils_lib, only : mkdir
1825 : integer, intent(in) :: nzlo, nzhi, nc, class_chem_id(:)
1826 : real(dp), intent(in), dimension(:,:) :: X1, X2, v
1827 :
1828 : integer :: k, j, ierr
1829 : character (len=strlen) :: fname
1830 : integer, parameter :: io = 33
1831 :
1832 : include 'formats'
1833 :
1834 : write(*,2) 'nzlo', nzlo
1835 : write(*,2) 'nzhi', nzhi
1836 :
1837 : call mkdir('plot_data/')
1838 : fname = 'plot_data/diffusion_plot.data'
1839 : ierr = 0
1840 : open(io, file=trim(fname), iostat=ierr)
1841 : if (ierr /= 0) then
1842 : write(*,*) 'failed to open ' // trim(fname)
1843 : call mesa_error(__FILE__,__LINE__,'write_plot_data')
1844 : end if
1845 :
1846 : write(io,'(a6)',advance='no') 'k'
1847 : do j=1,nc
1848 : write(io,'(a20)',advance='no') &
1849 : 'X1_' // trim(chem_isos% name(class_chem_id(j)))
1850 : write(io,'(a20)',advance='no') &
1851 : 'X2_' // trim(chem_isos% name(class_chem_id(j)))
1852 : end do
1853 : do j=1,nc
1854 : write(io,'(a20)',advance='no') &
1855 : 'logX1_' // trim(chem_isos% name(class_chem_id(j)))
1856 : write(io,'(a20)',advance='no') &
1857 : 'logX2_' // trim(chem_isos% name(class_chem_id(j)))
1858 : end do
1859 : do j=1,nc
1860 : write(io,'(a20)',advance='no') &
1861 : 'v_' // trim(chem_isos% name(class_chem_id(j)))
1862 : end do
1863 : write(io,'(A)')
1864 :
1865 : do k=nzlo,nzhi
1866 : write(io,'(i6)',advance='no') k
1867 : do j=1,nc
1868 : write(io,'(1pe20.12)',advance='no') X1(j,k)
1869 : write(io,'(1pe20.12)',advance='no') X2(j,k)
1870 : end do
1871 : do j=1,nc
1872 : write(io,'(1pe20.12)',advance='no') safe_log10(X1(j,k))
1873 : write(io,'(1pe20.12)',advance='no') safe_log10(X2(j,k))
1874 : end do
1875 : do j=1,nc
1876 : write(io,'(1pe20.12)',advance='no') v(j,k)
1877 : end do
1878 : write(io,'(A)')
1879 : end do
1880 :
1881 : close(io)
1882 :
1883 : end subroutine write_plot_data
1884 :
1885 :
1886 :
1887 : ! this routine sets up tables for 4 classes: h, he, o, and fe.
1888 : subroutine get_min_number_of_classes( &
1889 : species, chem_id, class, class_chem_id, class_name)
1890 : use chem_def
1891 : integer, parameter :: nc = diffusion_min_nc
1892 : integer, intent(in) :: species
1893 : integer, intent(in) :: chem_id(:) ! (species)
1894 : integer, intent(out) :: class(:) ! (species)
1895 : integer, intent(out) :: class_chem_id(:) ! (nc)
1896 : character (len=8), intent(out) :: class_name(:) ! (nc)
1897 : real(dp) :: A
1898 : integer :: i
1899 : integer, parameter :: c_h = 1, c_he = 2, c_o = 3, c_fe = 4
1900 : class_name(c_h) = 'c_h'
1901 : class_name(c_he) = 'c_he'
1902 : class_name(c_o) = 'c_o'
1903 : class_name(c_fe) = 'c_fe'
1904 : class_chem_id(c_h) = ih1
1905 : class_chem_id(c_he) = ihe4
1906 : class_chem_id(c_o) = io16
1907 : class_chem_id(c_fe) = ife56
1908 : do i=1,species
1909 : A = chem_isos% Z_plus_N(chem_id(i))
1910 : if (A < 3) then
1911 : class(i) = c_h
1912 : else if (A < 12) then
1913 : class(i) = c_he
1914 : else if (A < 20) then
1915 : class(i) = c_o
1916 : else
1917 : class(i) = c_fe
1918 : end if
1919 : end do
1920 : end subroutine get_min_number_of_classes
1921 :
1922 :
1923 : ! this routine sets up tables with a separate class for each isotope.
1924 : subroutine get_max_number_of_classes( &
1925 : species, chem_id, class, class_chem_id, class_name)
1926 : use chem_def
1927 : integer, intent(in) :: species, chem_id(:) ! (species)
1928 : integer, intent(out) :: class(:) ! (species)
1929 : integer, intent(out) :: class_chem_id(:) ! (species)
1930 : character (len=8), intent(out) :: class_name(:) ! (species)
1931 : integer :: i, ci
1932 : do i=1,species
1933 : ci = chem_id(i)
1934 : class_name(i) = 'c_' // trim(chem_isos% name(ci))
1935 : class(i) = i
1936 : class_chem_id(i) = ci
1937 : end do
1938 : end subroutine get_max_number_of_classes
1939 :
1940 :
1941 : subroutine get_diffusion_classes( &
1942 : nc, species, chem_id, class_chem_id, class_A_cutoff, &
1943 : class, class_name)
1944 : use chem_def
1945 : integer, intent(in) :: nc, species, chem_id(:) ! (species)
1946 : integer, intent(in) :: class_chem_id(:) ! (nc)
1947 : real(dp), intent(in) :: class_A_cutoff(:) ! (nc)
1948 : integer, intent(out) :: class(:) ! (species)
1949 : character (len=8), intent(out) :: class_name(:) ! (nc)
1950 : real(dp) :: A
1951 : integer :: i, j
1952 : do i=1,species
1953 : A = chem_isos% Z_plus_N(chem_id(i))
1954 : class(i) = nc
1955 : do j=1,nc-1
1956 : if (A < class_A_cutoff(j)) then
1957 : class(i) = j
1958 : exit
1959 : end if
1960 : end do
1961 : end do
1962 : do j=1,nc
1963 : class_name(j) = 'c_' // trim(chem_isos% name(class_chem_id(j)))
1964 : end do
1965 : end subroutine get_diffusion_classes
1966 :
1967 :
1968 0 : subroutine set_diffusion_classes( &
1969 0 : nc, species, chem_id, class_chem_id, class_A_max, use_full_net, &
1970 0 : class, class_name)
1971 : use chem_def
1972 : integer, intent(in) :: nc, species, chem_id(:) ! (species)
1973 : integer, intent(in) :: class_chem_id(:) ! (nc)
1974 : real(dp), intent(in) :: class_A_max(:) ! (nc)
1975 : logical, intent(in) :: use_full_net
1976 : integer, intent(out) :: class(:) ! (species)
1977 : character (len=8), intent(out) :: class_name(:) ! (nc)
1978 0 : real(dp) :: A
1979 : integer :: i, j
1980 0 : if(use_full_net) then
1981 0 : do i=1,species
1982 0 : class(i) = i
1983 : end do
1984 : else
1985 0 : do i=1,species
1986 0 : A = chem_isos% Z(chem_id(i)) + chem_isos% N(chem_id(i))
1987 0 : class(i) = nc
1988 0 : do j=1,nc-1
1989 0 : if (A <= class_A_max(j)) then
1990 0 : class(i) = j
1991 0 : exit
1992 : end if
1993 : end do
1994 : end do
1995 : end if
1996 0 : do j=1,nc
1997 0 : class_name(j) = 'c_' // trim(chem_isos% name(class_chem_id(j)))
1998 : end do
1999 0 : end subroutine set_diffusion_classes
2000 :
2001 :
2002 : real(dp) function adjust_timestep( &
2003 : vc_target, vc, vc_old, dt, dt_old, &
2004 : max_timestep_factor, min_timestep_factor) result(dt_next)
2005 : ! H211b "low pass" controller.
2006 : ! Soderlind Wang, J of Computational and Applied Math 185 (2006) 225 – 243.
2007 0 : use num_def
2008 : real(dp), intent(in) :: vc_target, vc, vc_old, dt, dt_old, &
2009 : max_timestep_factor, min_timestep_factor
2010 :
2011 : real(dp), parameter :: order = 1d0
2012 : real(dp), parameter :: beta1 = 0.25d0/order
2013 : real(dp), parameter :: beta2 = 0.25d0/order
2014 : real(dp), parameter :: alpha2 = 0.25d0
2015 : real(dp) :: vcratio, vcratio_prev, limtr
2016 :
2017 : include 'formats'
2018 :
2019 : dt_next = 1d99
2020 :
2021 : if (vc_old > 0 .and. dt_old > 0) then ! use 2 values to do "low pass" controller
2022 : vcratio = limiter(vc_target/vc)
2023 : vcratio_prev = limiter(vc_target/vc_old)
2024 : limtr = limiter(pow(vcratio,beta1) * pow(vcratio_prev,beta2) * pow(dt/dt_old,-alpha2))
2025 : dt_next = min(dt_next, dt*limtr)
2026 : else ! no history available, so fall back to the basic controller
2027 : dt_next = min(dt_next, dt*vc_target/vc)
2028 : end if
2029 :
2030 : if (max_timestep_factor > 0 .and. dt_next > max_timestep_factor*dt) then
2031 : dt_next = max_timestep_factor*dt
2032 : end if
2033 :
2034 : if (min_timestep_factor > 0 .and. dt_next < min_timestep_factor*dt) then
2035 : dt_next = min_timestep_factor*dt
2036 : end if
2037 :
2038 : contains
2039 :
2040 : real(dp) function limiter(x)
2041 : real(dp), intent(in) :: x
2042 : real(dp), parameter :: kappa = 2
2043 : ! for x >= 0 and kappa = 2, limiter value is between 0.07 and 4.14
2044 : ! for x = 1, limiter = 1
2045 : limiter = 1 + kappa*atan((x-1)/kappa)
2046 : end function limiter
2047 :
2048 : end function adjust_timestep
2049 :
2050 : end module diffusion
|
