Line data Source code
1 : ! ***********************************************************************
2 : !
3 : ! Copyright (C) 2010-2019 The MESA Team
4 : !
5 : ! This program is free software: you can redistribute it and/or modify
6 : ! it under the terms of the GNU Lesser General Public License
7 : ! as published by the Free Software Foundation,
8 : ! either version 3 of the License, or (at your option) any later version.
9 : !
10 : ! This program is distributed in the hope that it will be useful,
11 : ! but WITHOUT ANY WARRANTY; without even the implied warranty of
12 : ! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
13 : ! See the GNU Lesser General Public License for more details.
14 : !
15 : ! You should have received a copy of the GNU Lesser General Public License
16 : ! along with this program. If not, see <https://www.gnu.org/licenses/>.
17 : !
18 : ! ***********************************************************************
19 :
20 : module element_diffusion
21 :
22 : use star_private_def
23 : use const_def, only: dp, i8, pi4, amu, qe, secyer, no_mixing, phase_separation_mixing
24 :
25 : implicit none
26 :
27 : private
28 : public :: do_element_diffusion, finish_element_diffusion
29 :
30 : logical, parameter :: dbg = .false.
31 :
32 : contains
33 :
34 0 : subroutine do_element_diffusion(s, dt_in, ierr)
35 : ! return ierr /= 0 if cannot satisfy accuracy requirements
36 : use chem_def, only: chem_isos
37 : use chem_lib, only: chem_get_iso_id
38 : use star_utils, only: start_time, update_time
39 : use diffusion, only: &
40 : do_solve_diffusion, set_diffusion_classes, diffusion_min_nc
41 : type (star_info), pointer :: s
42 : real(dp), intent(in) :: dt_in
43 : integer, intent(out) :: ierr
44 :
45 : integer :: i, j, k, kk, nc, m, nzlo, nzhi, nz, species, iounit, &
46 : steps_used, total_num_iters, total_num_retries, cid
47 : integer(i8) :: time0
48 0 : real(dp) :: s1, s2, dqsum, dt, total, &
49 0 : gradT_mid, gradRho_mid, alfa, gradRho_face, chiRho_face, chiT_face
50 0 : real(dp) :: Amass, Zcharge, min_D_mix
51 :
52 : integer, dimension(:), allocatable :: &
53 0 : class, class_chem_id, mixing_type
54 : real(dp), dimension(:), allocatable :: &
55 0 : gamma, free_e, &
56 0 : dlnPdm_face, dlnT_dm_face, dlnRho_dm_face, &
57 0 : dlnPdm, dlnT_dm, dlnPdm_mid, dlnT_dm_mid, dlnRho_dm_mid, dm_hat
58 : real(dp), dimension(:,:), allocatable :: &
59 0 : X_init, X_final, typical_charge, &
60 0 : D_self, v_advection, v_total, &
61 0 : vlnP, vlnT, v_rad, g_rad, GT, xa_save
62 0 : real(dp), dimension(:,:,:), allocatable :: CD
63 0 : character (len=8), allocatable :: class_name(:)
64 :
65 :
66 : logical :: dumping
67 :
68 : include 'formats'
69 :
70 0 : ierr = 0
71 0 : dt = dt_in
72 0 : nz = s% nz
73 :
74 0 : s% num_diffusion_solver_steps = 0
75 :
76 0 : s% eps_diffusion(1:nz) = 0d0
77 0 : do k = 1, nz
78 0 : s% energy_start(k) = s% energy(k)
79 : ! This is just to track energy changes due to diffusion.
80 : ! s% energy_start gets reset later in do_struct_burn_mix before structure solve.
81 : end do
82 :
83 0 : if ((.not. s% do_element_diffusion) .or. dt < s% diffusion_dt_limit) then
84 0 : s% num_diffusion_solver_iters = 0 ! Flush diff iters to avoid crashing the timestep.
85 0 : s% edv(:,1:nz) = 0
86 0 : s% eps_WD_sedimentation(1:nz) = 0d0
87 0 : if (s% do_element_diffusion .and. s% report_ierr .and. dt < s% diffusion_dt_limit) &
88 0 : write(*,2) 'skip diffusion this step: dt < s% diffusion_dt_limit', &
89 0 : s% model_number, dt, s% diffusion_dt_limit
90 0 : return
91 : end if
92 :
93 0 : s% need_to_setvars = .true.
94 :
95 0 : if (s% use_other_diffusion_factor) then
96 0 : call s% other_diffusion_factor(s% id, ierr)
97 0 : if (ierr /= 0) then
98 0 : if (s% report_ierr) write(*,*) 'do_element_diffusion failed in other_diffusion_factor'
99 0 : return
100 : end if
101 : end if
102 :
103 0 : if (s% use_other_diffusion) then
104 0 : call s% other_diffusion(s% id, dt_in, ierr)
105 0 : if (ierr /= 0 .and. s% report_ierr) &
106 0 : write(*,*) 'do_element_diffusion failed in other_diffusion_factor'
107 0 : return
108 : end if
109 :
110 0 : if (s% doing_timing) call start_time(s, time0, total)
111 :
112 0 : nz = s% nz
113 0 : nzlo = 1
114 0 : nzhi = nz
115 0 : min_D_mix = s% D_mix_ignore_diffusion
116 0 : species = s% species
117 0 : if ( s% diffusion_use_full_net ) then
118 0 : if( species > 100 ) then
119 0 : stop "Net is too large for diffusion_use_full_net. max_num_diffusion_classes = 100"
120 : end if
121 0 : nc = species
122 : else
123 0 : nc = s% diffusion_num_classes
124 : end if
125 0 : m = nc+1
126 :
127 : allocate( &
128 0 : class(species), class_chem_id(nc), class_name(nc), CD(nc,nc,nz), mixing_type(nz), &
129 0 : gamma(nz), free_e(nz), dlnPdm_face(nz), dlnT_dm_face(nz), dlnRho_dm_face(nz), &
130 0 : dlnPdm(nz), dlnT_dm(nz), dlnPdm_mid(nz), dlnT_dm_mid(nz), dlnRho_dm_mid(nz), dm_hat(nz), &
131 : X_init(nc,nz), X_final(nc,nz), typical_charge(nc,nz), &
132 : D_self(nc,nz), v_advection(nc,nz), v_total(nc,nz), xa_save(species,nz), &
133 0 : vlnP(nc,nz), vlnT(nc,nz), v_rad(nc,nz), g_rad(nc,nz), GT(nc,nz))
134 :
135 0 : call set_extras(ierr)
136 0 : if (ierr /= 0) then
137 0 : if (s% report_ierr) write(*,*) 'do_element_diffusion failed in set_extras'
138 0 : return
139 : end if
140 :
141 :
142 : !write(*,1) 's% diffusion_min_T_at_surface', s% diffusion_min_T_at_surface
143 : !reset nzlo if necessary; nzlo=1 above
144 0 : nzlo = 1
145 0 : if (s% diffusion_min_T_at_surface > 0) then
146 : k = nzlo
147 0 : do while (s% T(k) < s% diffusion_min_T_at_surface)
148 0 : k = k+1
149 0 : if (k > nz) exit
150 0 : nzlo = k
151 : end do
152 : end if
153 :
154 : !write(*,2) 'diffusion_min_T_at_surface', nzlo
155 :
156 0 : if (s% diffusion_min_dq_ratio_at_surface > 0) then
157 0 : dqsum = sum(s% dq(1:nzlo))
158 : k = nzlo
159 0 : do while (dqsum < s% diffusion_min_dq_ratio_at_surface*s% dq(k+1))
160 0 : k = k+1
161 0 : if (k >= nz) exit
162 0 : dqsum = dqsum + s% dq(k)
163 0 : nzlo = k
164 : end do
165 : end if
166 :
167 : !write(*,2) 'before diffusion_min_dq_ratio_at_surface', nzlo, s% diffusion_min_dq_at_surface
168 :
169 0 : if (s% diffusion_min_dq_at_surface > 0) then
170 0 : dqsum = sum(s% dq(1:nzlo))
171 : k = nzlo
172 0 : do while (dqsum < s% diffusion_min_dq_at_surface)
173 0 : k = k+1
174 0 : if (k > nz) exit
175 : ! don't go across composition transition
176 0 : if (maxloc(s% xa(:,k),dim=1) /= maxloc(s% xa(:,k-1),dim=1)) exit
177 0 : dqsum = dqsum + s% dq(k)
178 0 : nzlo = k
179 : end do
180 : !write(*,2) 'after diffusion_min_dq_ratio_at_surface', nzlo, dqsum
181 : end if
182 :
183 : !write(*,3) 'nzlo mixing_type', nzlo, s% mixing_type(nzlo)
184 :
185 0 : if (s% D_mix(nzlo) > min_D_mix) then
186 0 : kk = nzlo
187 0 : do k = nzlo, nz-1
188 0 : if (maxloc(s% xa(:,k),dim=1) /= maxloc(s% xa(:,k-1),dim=1) .or. &
189 0 : s% D_mix(k+1) <= min_D_mix) then
190 0 : nzlo=k; exit
191 : end if
192 : end do
193 0 : nzlo = kk + (nzlo - kk)*4/5 ! back up into the convection zone
194 : end if
195 :
196 : !write(*,3) 'nzlo mixing_type', nzlo, s% mixing_type(nzlo)
197 : !write(*,3) 'nzlo-1 mixing_type', nzlo-1, s% mixing_type(nzlo-1)
198 : !write(*,*)
199 : !write(*,*)
200 :
201 : !reset nzhi if necessary; nzhi=nz above
202 0 : if (s% D_mix(nzhi) > min_D_mix) then
203 0 : do k = nz, 2, -1
204 0 : if (s% D_mix(k-1) <= min_D_mix) then
205 0 : nzhi=k; exit
206 : end if
207 : end do
208 0 : nzhi = (3*nzhi+nz)/4 ! back up some into the convection zone
209 : end if
210 :
211 0 : if (s% do_phase_separation .and. s% phase_separation_no_diffusion) then
212 : ! check for phase separation boundary, don't do diffusion deeper than that
213 0 : do k = 1,nzhi
214 0 : if(s% mixing_type(k) == phase_separation_mixing) then
215 0 : nzhi = k
216 0 : exit
217 : end if
218 : end do
219 : end if
220 :
221 0 : if(s% diffusion_use_full_net) then
222 0 : do j=1,nc
223 0 : class_chem_id(j) = s% chem_id(j) ! Just a 1-1 map between classes and chem_ids.
224 : end do
225 : else
226 0 : do j=1,nc
227 0 : cid = chem_get_iso_id(s% diffusion_class_representative(j))
228 0 : if (cid <= 0) then
229 : write(*,'(a,3x,i3)') 'bad entry for diffusion_class_representative: ' // &
230 0 : trim(s% diffusion_class_representative(j)), j
231 0 : return
232 : end if
233 0 : class_chem_id(j) = cid
234 : end do
235 : end if
236 :
237 : call set_diffusion_classes( &
238 : nc, species, s% chem_id, class_chem_id, s% diffusion_class_A_max, s% diffusion_use_full_net, &
239 0 : class, class_name)
240 :
241 0 : s% diffusion_call_number = s% diffusion_call_number + 1
242 0 : dumping = (s% diffusion_call_number == s% diffusion_dump_call_number)
243 :
244 0 : if ( s% diffusion_use_full_net ) then
245 0 : s% diffusion_calculates_ionization = .true. ! class_typical_charges can't be used, so make sure they aren't.
246 : end if
247 :
248 0 : if (.not. s% diffusion_calculates_ionization) then
249 0 : do j=1,nc
250 0 : typical_charge(j,1:nz) = s% diffusion_class_typical_charge(j)
251 : end do
252 : end if
253 :
254 0 : do k=1,nz
255 0 : do j=1,species
256 0 : xa_save(j,k) = s% xa(j,k)
257 : end do
258 : end do
259 :
260 0 : mixing_type(1:nz) = no_mixing
261 :
262 0 : do k=1,nz-1
263 0 : s1 = dlnPdm(k)
264 0 : s2 = dlnPdm(k+1)
265 0 : if (s1*s2 <= 0) then
266 0 : dlnPdm_mid(k) = 0
267 : else
268 0 : dlnPdm_mid(k) = 2*s1*s2/(s1+s2)
269 : end if
270 0 : gradT_mid = 0.5d0*(s% gradT(k) + s% gradT(k+1))
271 0 : dlnT_dm_mid(k) = gradT_mid*dlnPdm_mid(k)
272 0 : gradRho_mid = (1 - s% chiT(k)*gradT_mid)/s% chiRho(k)
273 0 : dlnRho_dm_mid(k) = gradRho_mid*dlnPdm_mid(k)
274 : end do
275 0 : dlnPdm_mid(nz) = dlnPdm(nz)
276 0 : dlnT_dm_mid(nz) = dlnT_dm(nz)
277 0 : gradRho_mid = (1 - s% chiT(nz)*s% gradT(nz))/s% chiRho(nz)
278 0 : dlnRho_dm_mid(nz) = gradRho_mid*dlnPdm_mid(nz)
279 :
280 0 : do k=2,nz
281 0 : dlnPdm_face(k) = dlnPdm(k)
282 0 : dlnT_dm_face(k) = dlnT_dm(k)
283 0 : alfa = s% dm(k-1)/(s% dm(k) + s% dm(k-1))
284 0 : chiT_face = alfa*s% chiT(k) + (1-alfa)*s% chiT(k-1)
285 0 : chiRho_face = alfa*s% chiRho(k) + (1-alfa)*s% chiRho(k-1)
286 0 : gradRho_face = (1 - chiT_face*s% gradT(k))/chiRho_face
287 0 : dlnRho_dm_face(k) = gradRho_face*dlnPdm(k)
288 : end do
289 0 : dlnPdm_face(1) = 0
290 0 : dlnT_dm_face(1) = 0
291 0 : dlnRho_dm_face(1) = 0
292 :
293 0 : if (dumping) call dump_diffusion_info
294 :
295 : ! print *, "About to call do_solve_diffusion. Here's the class structure"
296 : ! print *, "chem_id: ", s% chem_id
297 : ! print *, "class_chem_id: ", class_chem_id
298 : ! print *, "class: ", class
299 : ! print *, "class name: ", class_name
300 :
301 : ! args are at cell center points.
302 : !if (s% show_diffusion_info) write(*,*) 'call solve_diffusion'
303 : !write(*,4) 'call do_solve_diffusion nzlo nzhi nz', nzlo, nzhi, nz, &
304 : ! sum(s% xa(1,1:nzlo))
305 : call do_solve_diffusion( &
306 : s, nz, species, nc, m, class, class_chem_id, s% net_iso, s% chem_id, &
307 : s% abar, s% ye, free_e, s% dm_bar, s% dm, &
308 : s% T, s% lnT, s% rho, s% lnd, s% rmid, &
309 : dlnPdm_mid, dlnT_dm_mid, dlnRho_dm_mid, &
310 : s% L, s% r, dlnPdm_face, dlnT_dm_face, dlnRho_dm_face, &
311 : s% diffusion_use_iben_macdonald, dt, s% diffusion_dt_div_timescale, &
312 : s% diffusion_steps_hard_limit, s% diffusion_iters_hard_limit, &
313 : s% diffusion_max_iters_per_substep, &
314 : s% diffusion_calculates_ionization, typical_charge, &
315 : s% diffusion_nsmooth_typical_charge, &
316 : s% diffusion_min_T_for_radaccel, s% diffusion_max_T_for_radaccel, &
317 : s% diffusion_min_Z_for_radaccel, s% diffusion_max_Z_for_radaccel, &
318 : s% diffusion_screening_for_radaccel, &
319 : s% op_mono_data_path, s% op_mono_data_cache_filename, &
320 : s% diffusion_v_max, s% R_center, &
321 : gamma, s% diffusion_gamma_full_on, s% diffusion_gamma_full_off, &
322 : s% diffusion_T_full_on, s% diffusion_T_full_off, &
323 : s% diffusion_class_factor, s% xa, &
324 : steps_used, total_num_iters, total_num_retries, nzlo, nzhi, X_init, X_final, &
325 : D_self, v_advection, v_total, vlnP, vlnT, v_rad, g_rad, &
326 0 : s% E_field, s% g_field_element_diffusion, ierr )
327 0 : s% num_diffusion_solver_steps = steps_used
328 0 : s% num_diffusion_solver_iters = total_num_iters
329 :
330 0 : if (dbg .or. s% show_diffusion_info .or. ierr /= 0) then
331 0 : if (ierr == 0) then
332 : write(*,'(a,f6.3,3x,a,1pe10.3,3x,99(a,i5,3x))') &
333 0 : 'log_dt', log10(s% dt/secyer), 'age', s% star_age, 'model', s% model_number, &
334 0 : 'iters', total_num_iters, 'steps', steps_used, 'retries', total_num_retries, &
335 0 : 'nzlo', nzlo, 'nzhi', nzhi, 'n', nzhi-nzlo+1, 'nz', nz, &
336 0 : 'diffusion_call_number', s% diffusion_call_number
337 : else
338 : write(*,'(a,2x,f10.3,3x,99(a,i5,3x))') &
339 0 : 'do_solve_diffusion FAILED: log_dt', log10(s% dt/secyer), 'model', s% model_number, &
340 0 : 'iters', total_num_iters, 'steps', steps_used, 'retries', total_num_retries, &
341 0 : 'nzlo', nzlo, 'nzhi', nzhi, 'n', nzhi-nzlo+1, 'nz', nz, &
342 0 : 'diffusion_call_number', s% diffusion_call_number
343 : end if
344 : end if
345 :
346 0 : if (ierr /= 0) then
347 0 : do k=1,nz
348 0 : do j=1,species
349 0 : s% xa(j,k) = xa_save(j,k)
350 : end do
351 : end do
352 0 : if (s% report_ierr) then
353 0 : write(*, *)
354 0 : write(*, *) 'solve_diffusion returned false'
355 0 : write(*, *) 's% model_number', s% model_number
356 0 : write(*, *) 's% diffusion_call_number', s% diffusion_call_number
357 0 : write(*, *)
358 : end if
359 : end if
360 :
361 0 : if (dumping) call mesa_error(__FILE__,__LINE__,'debug: dump_diffusion_info')
362 :
363 0 : do k=nzlo+1,nzhi
364 0 : do j=1,species
365 0 : i = class(j)
366 0 : s% diffusion_D_self(j,k) = D_self(i,k)
367 0 : s% edv(j,k) = v_total(i,k)
368 0 : s% v_rad(j,k) = v_rad(i,k)
369 0 : s% g_rad(j,k) = g_rad(i,k)
370 0 : s% typical_charge(j,k) = typical_charge(i,k)
371 0 : s% diffusion_dX(j,k) = s% xa(j,k) - xa_save(j,k)
372 : end do
373 : end do
374 :
375 0 : if(s% do_diffusion_heating .and. s% do_WD_sedimentation_heating) then
376 0 : write(*,*) "do_diffusion_heating is incompatible with do_WD_sedimentation_heating"
377 0 : write(*,*) "at least one of these options must be set to .false."
378 0 : call mesa_error(__FILE__,__LINE__,'do_element_diffusion')
379 : end if
380 :
381 0 : s% eps_WD_sedimentation(1:nz) = 0d0
382 :
383 0 : if(s% do_WD_sedimentation_heating) then
384 0 : do k=nzlo+1,nzhi
385 : ! loop over all ion species
386 0 : do i = 1,s% species
387 : ! limit heating to species with significant mass fractions
388 0 : if(s% xa(i,k) > s% min_xa_for_WD_sedimentation_heating) then
389 0 : Amass = chem_isos% Z_plus_N(s% chem_id(i))
390 0 : Zcharge = chem_isos% Z(s% chem_id(i))
391 : s% eps_WD_sedimentation(k) = s% eps_WD_sedimentation(k) + &
392 : s% eps_WD_sedimentation_factor * &
393 : ( Amass - Zcharge * qe * s% E_field(k)/(amu * s% g_field_element_diffusion(k)) ) * &
394 : sign(1d0,-1d0*s% edv(i,k)) * min( s% diffusion_v_max, abs(s% edv(i,k)) ) * s% xa(i,k) * &
395 0 : s% g_field_element_diffusion(k) / Amass
396 : end if
397 : end do
398 : ! For diagnostics:
399 0 : if( .false. .and. mod(k,100) == 0) then
400 : write(*,*) "Zone: ", k
401 : write(*,*) "eE/mg = ", (qe * s% E_field(k)/(amu * s% g_field_element_diffusion(k)))
402 : write(*,*) "Net Force on Ne22 (units of mp*g) = ", &
403 : (22d0 - 10d0 * qe * s% E_field(k)/(amu * s% g_field_element_diffusion(k)))
404 : write(*,*) "g_field_element_diffusion/grav = ", s% g_field_element_diffusion(k) / s% grav(k)
405 : end if
406 : end do
407 : end if
408 :
409 :
410 0 : do k=1,nzlo
411 0 : do j=1,species
412 0 : s% diffusion_D_self(j,k) = s% diffusion_D_self(j,nzlo+1)
413 0 : s% edv(j,k) = 0d0 ! s% edv(j,nzlo+1)
414 0 : s% v_rad(j,k) = s% v_rad(j,nzlo+1)
415 0 : s% g_rad(j,k) = s% g_rad(j,nzlo+1)
416 0 : s% typical_charge(j,k) = s% typical_charge(j,nzlo+1)
417 0 : s% diffusion_dX(j,k) = s% xa(j,k) - xa_save(j,k)
418 : end do
419 0 : s% E_field(k) = 0d0 ! s% E_field(nzlo+1)
420 0 : s% g_field_element_diffusion(k) = 0d0 ! s% g_field_element_diffusion(nzlo+1)
421 : end do
422 :
423 0 : do k=nzhi+1,nz
424 0 : do j=1,species
425 0 : s% diffusion_D_self(j,k) = s% diffusion_D_self(j,nzhi)
426 0 : s% edv(j,k) = 0d0 ! s% edv(j,nzhi)
427 0 : s% v_rad(j,k) = s% v_rad(j,nzhi)
428 0 : s% g_rad(j,k) = s% g_rad(j,nzhi)
429 0 : s% typical_charge(j,k) = s% typical_charge(j,nzhi)
430 0 : s% diffusion_dX(j,k) = s% xa(j,k) - xa_save(j,k)
431 : end do
432 0 : s% E_field(k) = 0d0 ! s% E_field(nzhi)
433 0 : s% g_field_element_diffusion(k) = 0d0 ! s% g_field_element_diffusion(nzhi)
434 : end do
435 :
436 0 : if (s% doing_timing) call update_time(s, time0, total, s% time_element_diffusion)
437 :
438 : contains
439 :
440 : subroutine check_xa_sums(ierr)
441 : integer, intent(out) :: ierr
442 : integer :: k
443 : include 'formats'
444 : do k=1, nz
445 : if (abs(sum(s% xa(1:species, k)) - 1d0) > 1d-3) then
446 : write(*,*) 'k', k
447 : write(*,1) 'sum', sum(s% xa(1:species, k))
448 : write(*,1) 'sum-1d0', sum(s% xa(1:species, k))-1d0
449 : write(*,1) 'abs(sum-1d0)', abs(sum(s% xa(1:species, k))-1d0)
450 : do j=1,species
451 : write(*,2) 's% xa(j,k)', j, s% xa(j, k)
452 : end do
453 : ierr = -1
454 : end if
455 : end do
456 0 : end subroutine check_xa_sums
457 :
458 0 : subroutine dump_diffusion_info
459 : use utils_lib
460 : use chem_def, only: chem_isos
461 : integer :: i, k, ierr
462 :
463 0 : ierr = 0
464 0 : write(*, *)
465 0 : write(*, *) 'write diffusion.data'
466 0 : write(*, *)
467 0 : open(newunit=iounit, file='diffusion.data', action='write', status='replace', iostat=ierr)
468 0 : if (ierr /= 0) then
469 0 : write(*, *) 'failed to open diffusion dump file'
470 : return
471 : end if
472 :
473 : ! args
474 0 : write(iounit, '(99i20)') nz, nzlo, nzhi, species, nc, &
475 0 : s% diffusion_steps_hard_limit, &
476 0 : s% diffusion_nsmooth_typical_charge
477 :
478 0 : do i=1,species
479 0 : write(iounit,*) trim(chem_isos% name(s% chem_id(i)))
480 : end do
481 :
482 0 : write(iounit, '(99i20)') class(1:species)
483 :
484 0 : write(iounit, '(99i20)') chem_isos% Z(s% chem_id(1:species))
485 :
486 0 : write(iounit, '(99i20)') chem_isos% Z_plus_N(s% chem_id(1:species))
487 :
488 0 : write(iounit, '(99e22.10)') chem_isos% W(s% chem_id(1:species))
489 :
490 0 : do i=1,nc
491 0 : write(iounit, '(a)') trim(chem_isos% name(class_chem_id(i)))
492 : end do
493 :
494 0 : do i=1,nc
495 0 : write(iounit, '(a)') trim(class_name(i))
496 : end do
497 :
498 0 : if (s% diffusion_calculates_ionization) then
499 0 : write(iounit,*) 1
500 : else
501 0 : write(iounit,*) 0
502 : end if
503 :
504 0 : if (s% diffusion_use_iben_macdonald) then
505 0 : write(iounit,*) 1
506 : else
507 0 : write(iounit,*) 0
508 : end if
509 :
510 0 : if (s% diffusion_screening_for_radaccel) then
511 0 : write(iounit,*) 1
512 : else
513 0 : write(iounit,*) 0
514 : end if
515 :
516 : write(iounit, '(99(1pd26.16))') &
517 0 : s% mstar, s% dt, &
518 0 : s% diffusion_v_max, s% diffusion_gamma_full_on, s% diffusion_gamma_full_off, &
519 0 : s% diffusion_T_full_on, s% diffusion_T_full_off, &
520 0 : s% diffusion_max_T_for_radaccel, s% diffusion_dt_div_timescale, &
521 0 : s% R_center
522 :
523 0 : do k=1, nz
524 : write(iounit, '(99(1pd26.16))') &
525 0 : gamma(k), s% abar(k), s% ye(k), free_e(k), s% dm_bar(k), s% dm(k), &
526 0 : s% T(k), s% lnT(k), s% Rho(k), s% lnd(k), s% rmid(k), &
527 0 : dlnPdm_mid(k), dlnT_dm_mid(k), dlnRho_dm_mid(k), &
528 0 : s% L(k), s% r(k), dlnPdm_face(k), dlnT_dm_face(k), dlnRho_dm_face(k)
529 0 : write(iounit, '(99(1pd26.16))') s% xa(1:species, k)
530 : end do
531 :
532 0 : close(iounit)
533 :
534 0 : end subroutine dump_diffusion_info
535 :
536 :
537 0 : subroutine set_extras(ierr)
538 : integer, intent(out) :: ierr
539 : integer :: k, op_err
540 0 : op_err = 0
541 0 : ierr = 0
542 0 : !$OMP PARALLEL DO PRIVATE(k, op_err) SCHEDULE(dynamic,2)
543 : do k=1,nz
544 : call set1_extras(k, op_err)
545 : if (op_err /= 0) ierr = op_err
546 : end do
547 : !$OMP END PARALLEL DO
548 0 : end subroutine set_extras
549 :
550 :
551 0 : subroutine set1_extras(k,ierr)
552 : integer, intent(in) :: k
553 : integer, intent(out) :: ierr
554 0 : real(dp) :: grav, area, P_face
555 0 : ierr = 0
556 0 : free_e(k) = exp(s% lnfree_e(k))
557 0 : gamma(k) = s% gam(k)
558 0 : if (k==1) then
559 0 : dlnPdm(k) = 0; dlnT_dm(k) = 0; return
560 : end if
561 0 : grav = -s% cgrav(k)*s% m(k)/s% r(k)**2
562 0 : area = pi4*s% r(k)**2
563 0 : P_face = 0.5d0*(s% Peos(k) + s% Peos(k-1))
564 0 : dlnPdm(k) = grav/(area*P_face) ! estimate based on QHSE
565 0 : dlnT_dm(k) = s% gradT(k)*dlnPdm(k)
566 : end subroutine set1_extras
567 :
568 :
569 : end subroutine do_element_diffusion
570 :
571 0 : subroutine finish_element_diffusion(s,dt)
572 : type (star_info), pointer :: s
573 : real(dp), intent(in) :: dt
574 : integer :: k
575 :
576 0 : do k=1,s% nz
577 0 : s% eps_diffusion(k) = (s% energy_start(k) - s% energy(k))/dt
578 : end do
579 :
580 0 : end subroutine finish_element_diffusion
581 :
582 : end module element_diffusion
|
