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 net_burn_const_P
21 : use const_def, only: dp, i8, ln10
22 : use chem_def
23 : use math_lib
24 : use net_def
25 : use rates_def, only: num_rvs
26 : use mtx_def
27 : use utils_lib, only: fill_with_NaNs,fill_with_NaNs_2D
28 :
29 : implicit none
30 :
31 : contains
32 :
33 2 : subroutine burn_1_zone_const_P( &
34 : net_handle, eos_handle, num_isos, num_reactions, &
35 : which_solver, starting_temp, starting_x, clip, &
36 : ntimes, times, log10Ps_f1, &
37 4 : rate_factors, weak_rate_factor, &
38 : reaction_Qs, reaction_neuQs, screening_mode, &
39 : h, max_step_size, max_steps, rtol, atol, itol, x_min, x_max, which_decsol, &
40 : caller_id, solout, iout, &
41 2 : ending_x, ending_temp, ending_rho, ending_lnS, initial_rho, initial_lnS, &
42 : nfcn, njac, nstep, naccpt, nrejct, time_doing_net, time_doing_eos, ierr)
43 : use num_def
44 : use num_lib
45 : use mtx_lib
46 : use rates_def, only: rates_reaction_id_max
47 :
48 : integer, intent(in) :: net_handle, eos_handle
49 : integer, intent(in) :: num_isos
50 : integer, intent(in) :: num_reactions
51 : real(dp), pointer, intent(in) :: starting_x(:) ! (num_isos)
52 : real(dp), intent(in) :: starting_temp
53 : logical, intent(in) :: clip ! if true, set negative x's to zero during burn.
54 :
55 : integer, intent(in) :: which_solver ! as defined in num_def.f
56 : integer, intent(in) :: ntimes ! ending time is times(num_times); starting time is 0
57 : real(dp), pointer, intent(in) :: times(:) ! (num_times)
58 : real(dp), pointer, intent(in) :: log10Ps_f1(:) ! =(4,numtimes) interpolant for log10P(time)
59 :
60 : real(dp), intent(in) :: rate_factors(:) ! (num_reactions)
61 : real(dp), intent(in) :: weak_rate_factor
62 : real(dp), pointer, intent(in) :: reaction_Qs(:) ! (rates_reaction_id_max)
63 : real(dp), pointer, intent(in) :: reaction_neuQs(:) ! (rates_reaction_id_max)
64 : integer, intent(in) :: screening_mode
65 :
66 : ! args to control the solver -- see num/public/num_isolve.dek
67 : real(dp), intent(inout) :: h
68 : real(dp), intent(in) :: max_step_size ! maximal step size.
69 : integer, intent(in) :: max_steps ! maximal number of allowed steps.
70 : ! absolute and relative error tolerances
71 : real(dp), intent(inout) :: rtol(*) ! relative error tolerance(s)
72 : real(dp), intent(inout) :: atol(*) ! absolute error tolerance(s)
73 : integer, intent(in) :: itol ! switch for rtol and atol
74 : real(dp), intent(in) :: x_min, x_max ! bounds on allowed values
75 : integer, intent(in) :: which_decsol ! from mtx_def
76 : integer, intent(in) :: caller_id
77 : interface ! subroutine called after each successful step
78 : include "num_solout.dek"
79 : end interface
80 : integer, intent(in) :: iout
81 : real(dp), intent(inout) :: ending_x(:)
82 : real(dp), intent(out) :: ending_temp, ending_rho, ending_lnS, initial_rho, initial_lnS
83 : integer, intent(out) :: nfcn ! number of function evaluations
84 : integer, intent(out) :: njac ! number of jacobian evaluations
85 : integer, intent(out) :: nstep ! number of computed steps
86 : integer, intent(out) :: naccpt ! number of accepted steps
87 : integer, intent(out) :: nrejct ! number of rejected steps
88 : real(dp), intent(inout) :: time_doing_net
89 : ! if < 0, then ignore
90 : ! else on return has input value plus time spent doing eval_net
91 : real(dp), intent(inout) :: time_doing_eos
92 : ! if < 0, then ignore
93 : ! else on return has input value plus time spent doing eos
94 : integer, intent(out) :: ierr
95 :
96 2 : type (Net_Info) :: n
97 : type (Net_General_Info), pointer :: g
98 : integer :: ijac, nzmax, isparse, mljac, mujac, imas, mlmas, mumas, lrd, lid, &
99 : lout, liwork, lwork, i, j, lrpar, lipar, idid, nvar
100 2 : integer, pointer :: ipar(:), iwork(:), ipar_burn_const_P_decsol(:)
101 2 : real(dp), pointer :: rpar(:), work(:), v(:), rpar_burn_const_P_decsol(:)
102 2 : real(dp) :: t, lgT, lgRho, tend
103 :
104 : include 'formats'
105 :
106 44 : ending_x = 0
107 2 : ending_temp = 0
108 2 : ending_rho = 0
109 2 : ending_lnS = 0
110 2 : nfcn = 0
111 2 : njac = 0
112 2 : nstep = 0
113 2 : naccpt = 0
114 2 : nrejct = 0
115 : ierr = 0
116 :
117 2 : nvar = num_isos + 1
118 :
119 2 : call get_net_ptr(net_handle, g, ierr)
120 2 : if (ierr /= 0) then
121 : return
122 : end if
123 :
124 2 : if (g% num_isos /= num_isos) then
125 0 : write(*,*) 'invalid num_isos', num_isos
126 0 : return
127 : end if
128 :
129 2 : if (g% num_reactions /= num_reactions) then
130 0 : write(*,*) 'invalid num_reactions', num_reactions
131 0 : return
132 : end if
133 :
134 2 : if (which_decsol == lapack) then
135 2 : nzmax = 0
136 2 : isparse = 0
137 : call lapack_work_sizes(nvar, lrd, lid)
138 : else
139 0 : write(*,1) 'net 1 zone burn const P: unknown value for which_decsol', which_decsol
140 0 : call mesa_error(__FILE__,__LINE__)
141 : end if
142 :
143 2 : ijac = 1
144 2 : mljac = nvar ! square matrix
145 2 : mujac = nvar
146 :
147 2 : imas = 0
148 2 : mlmas = 0
149 2 : mumas = 0
150 :
151 2 : lout = 0
152 :
153 : call isolve_work_sizes(nvar, nzmax, imas, mljac, mujac, mlmas, mumas, liwork, lwork)
154 :
155 2 : lipar = burn_lipar
156 2 : lrpar = burn_const_P_lrpar
157 :
158 : allocate(v(nvar), iwork(liwork), work(lwork), rpar(lrpar), ipar(lipar), &
159 2 : ipar_burn_const_P_decsol(lid), rpar_burn_const_P_decsol(lrd), stat=ierr)
160 2 : if (ierr /= 0) then
161 0 : write(*, *) 'allocate ierr', ierr
162 0 : return
163 : end if
164 :
165 2 : if (g% fill_arrays_with_nans) then
166 0 : call fill_with_NaNs(v)
167 0 : call fill_with_NaNs(work)
168 0 : call fill_with_NaNs(rpar)
169 0 : call fill_with_NaNs(rpar_burn_const_P_decsol)
170 : end if
171 :
172 :
173 2 : i = burn_const_P_lrpar
174 :
175 2 : ipar(i_burn_caller_id) = caller_id
176 2 : ipar(i_net_handle) = net_handle
177 2 : ipar(i_eos_handle) = eos_handle
178 2 : ipar(i_screening_mode) = screening_mode
179 2 : ipar(i_sparse_format) = isparse
180 2 : if (clip) then
181 0 : ipar(i_clip) = 1
182 : else
183 2 : ipar(i_clip) = 0
184 : end if
185 :
186 132 : iwork = 0
187 2642 : work = 0
188 :
189 2 : t = 0
190 2 : tend = times(ntimes)
191 :
192 2 : rpar(r_burn_const_P_rho) = 1d-99 ! dummy value, will be calculated later
193 2 : rpar(r_burn_const_P_pressure) = exp10(log10Ps_f1(1)) ! no interpolation yet
194 2 : rpar(r_burn_const_P_temperature) = starting_temp
195 2 : rpar(r_burn_const_P_init_rho) = -1d99
196 2 : rpar(r_burn_const_P_init_lnS) = -1d99
197 2 : rpar(r_burn_const_P_time_net) = time_doing_net
198 2 : rpar(r_burn_const_P_time_eos) = time_doing_eos
199 :
200 86 : v(1:num_isos) = starting_x(1:num_isos)
201 2 : v(nvar) = log(starting_temp)
202 :
203 2 : if (which_decsol == lapack) then
204 2 : call do_isolve(lapack_decsol, null_decsols, ierr)
205 : else
206 0 : write(*,*) 'unknown value for which_decsol', which_decsol
207 0 : call mesa_error(__FILE__,__LINE__)
208 : end if
209 2 : if (ierr /= 0) then
210 0 : call dealloc
211 0 : return
212 : end if
213 :
214 2 : nfcn = iwork(14)
215 2 : njac = iwork(15)
216 2 : nstep = iwork(16)
217 2 : naccpt = iwork(17)
218 2 : nrejct = iwork(18)
219 2 : time_doing_net = rpar(r_burn_const_P_time_net)
220 2 : time_doing_eos = rpar(r_burn_const_P_time_eos)
221 :
222 44 : ending_x(1:num_isos) = v(1:num_isos)
223 2 : ending_temp = exp(v(nvar))
224 2 : ending_rho = rpar(r_burn_const_P_rho)
225 2 : ending_lnS = rpar(r_burn_const_P_lnS)
226 2 : initial_rho = rpar(r_burn_const_P_init_rho)
227 2 : initial_lnS = rpar(r_burn_const_P_init_lnS)
228 :
229 : if (ierr /= 0) then
230 : write(*, '(a30,i10)') 'nfcn', nfcn
231 : write(*, '(a30,i10)') 'njac', njac
232 : write(*, '(a30,i10)') 'nstep', nstep
233 : write(*, '(a30,i10)') 'naccpt', naccpt
234 : write(*, '(a30,i10)') 'nrejct', nrejct
235 : call mesa_error(__FILE__,__LINE__)
236 : end if
237 :
238 6 : call dealloc
239 :
240 :
241 : contains
242 :
243 :
244 2 : subroutine dealloc
245 2 : deallocate(iwork, work, rpar, ipar, ipar_burn_const_P_decsol, rpar_burn_const_P_decsol)
246 2 : end subroutine dealloc
247 :
248 :
249 2 : subroutine do_isolve(decsol, decsols, ierr)
250 : interface
251 : include "mtx_decsol.dek"
252 : include "mtx_decsols.dek"
253 : end interface
254 : integer, intent(out) :: ierr
255 : integer :: caller_id, nvar_blk, nz_blk
256 : real(dp), dimension(:), pointer :: &
257 2 : lblk, dblk, ublk, uf_lblk, uf_dblk, uf_ublk
258 :
259 2 : nullify(lblk, dblk, ublk, uf_lblk, uf_dblk, uf_ublk)
260 2 : caller_id = 0
261 2 : nvar_blk = 0
262 2 : nz_blk = 0
263 : include 'formats'
264 2 : ierr = 0
265 : ! NOTE: don't use x_max for const_P since x includes other variables
266 : call isolve( &
267 : which_solver, nvar, burn_derivs, t, v, tend, &
268 : h, max_step_size, max_steps, &
269 : rtol, atol, itol, x_min, 1d99, &
270 : burn_jacob, ijac, null_sjac, nzmax, isparse, mljac, mujac, &
271 : null_mas, imas, mlmas, mumas, &
272 : solout, iout, &
273 : decsol, decsols, null_decsolblk, &
274 : lrd, rpar_burn_const_P_decsol, lid, ipar_burn_const_P_decsol, &
275 : caller_id, nvar_blk, nz_blk, &
276 : lblk, dblk, ublk, uf_lblk, uf_dblk, uf_ublk, &
277 : null_fcn_blk_dble, null_jac_blk_dble, &
278 : work, lwork, iwork, liwork, &
279 : lrpar, rpar, lipar, ipar, &
280 2 : lout, idid)
281 2 : if (idid < 0) ierr = -1
282 2 : end subroutine do_isolve
283 :
284 :
285 4246 : subroutine burn_derivs(nvar, t, h, v, f, lrpar, rpar, lipar, ipar, ierr)
286 : integer, intent(in) :: nvar, lrpar, lipar
287 : real(dp), intent(in) :: t, h
288 : real(dp), intent(inout) :: v(:) ! (nvar)
289 : real(dp), intent(inout) :: f(:) ! (nvar) ! dvdt
290 : integer, intent(inout), pointer :: ipar(:) ! (lipar)
291 : real(dp), intent(inout), pointer :: rpar(:) ! (lrpar)
292 : integer, intent(out) :: ierr
293 : integer, parameter :: ld_dfdv = 0
294 8492 : real(dp) :: dfdv(ld_dfdv,nvar)
295 : ierr = 0
296 4246 : call burn_jacob(nvar, t, h, v, f, dfdv, ld_dfdv, lrpar, rpar, lipar, ipar, ierr)
297 4246 : end subroutine burn_derivs
298 :
299 :
300 5088 : subroutine burn_jacob( &
301 5088 : nvar, time, h, v, f, dfdv, ld_dfdv, lrpar, rpar, lipar, ipar, ierr)
302 : use chem_lib, only: basic_composition_info
303 : use net_eval, only: eval_net
304 : use eos_def
305 : use eos_lib, only: Radiation_Pressure, eosPT_get
306 : use rates_def, only: rates_reaction_id_max
307 :
308 : integer, intent(in) :: nvar, ld_dfdv, lrpar, lipar
309 : real(dp), intent(in) :: time, h
310 : real(dp), intent(inout) :: v(:)
311 : real(dp), intent(inout) :: f(:), dfdv(:,:)
312 : integer, intent(inout), pointer :: ipar(:) ! (lipar)
313 : real(dp), intent(inout), pointer :: rpar(:) ! (lrpar)
314 : integer, intent(out) :: ierr
315 :
316 : integer :: net_handle, num_reactions, eos_handle
317 : real(dp) :: &
318 20352 : abar, zbar, z2bar, z53bar, ye, mass_correction, sumx, T, logT, &
319 5088 : rho, logRho, pressure, Pgas, Prad, lgPgas, &
320 111936 : dlnT_dt, x(nvar-1)
321 5088 : real(dp) :: eta, d_eta_dlnT, d_eta_dlnRho
322 5088 : real(dp) :: eps_neu_total, eps_nuc
323 5088 : real(dp) :: d_eps_nuc_dT
324 5088 : real(dp) :: d_eps_nuc_dRho
325 111936 : real(dp) :: d_eps_nuc_dx(nvar-1)
326 111936 : real(dp) :: dxdt(nvar-1)
327 111936 : real(dp) :: d_dxdt_dRho(nvar-1)
328 111936 : real(dp) :: d_dxdt_dT(nvar-1)
329 2355744 : real(dp) :: d_dxdt_dx(nvar-1, nvar-1)
330 127200 : real(dp), target :: eps_nuc_categories(num_categories)
331 : logical :: rates_only, skip_jacobian
332 : integer :: screening_mode, i, num_isos
333 : integer(i8) :: time0, time1, clock_rate
334 :
335 10176 : real(dp) :: xh, Y, z, Cp, rate_limit
336 5088 : real(dp) :: dlnRho_dlnPgas_const_T, dlnRho_dlnT_const_Pgas
337 5088 : real(dp) :: dlnRho_dlnT_const_P, d_epsnuc_dlnT_const_P, d_Cp_dlnT
338 137376 : real(dp) :: res(num_eos_basic_results)
339 137376 : real(dp) :: d_dlnRho_const_T(num_eos_basic_results)
340 137376 : real(dp) :: d_dlnT_const_Rho(num_eos_basic_results)
341 325632 : real(dp) :: d_dxa_const_TRho(num_eos_d_dxa_results, nvar-1)
342 5088 : integer, pointer :: net_iso(:), chem_id(:)
343 :
344 : type (Net_General_Info), pointer :: g
345 5088 : real(dp), pointer, dimension(:) :: actual_Qs, actual_neuQs
346 5088 : logical, pointer :: from_weaklib(:)
347 5088 : actual_Qs => null()
348 5088 : actual_neuQs => null()
349 5088 : from_weaklib => null()
350 :
351 : include 'formats'
352 :
353 5088 : num_isos = nvar-1
354 :
355 5088 : ierr = 0
356 117024 : f = 0
357 524552 : dfdv = 0
358 :
359 5088 : eos_handle = ipar(i_eos_handle)
360 :
361 5088 : net_handle = ipar(i_net_handle)
362 5088 : call get_net_ptr(net_handle, g, ierr)
363 5088 : if (ierr /= 0) then
364 0 : write(*,*) 'invalid handle for eval_net -- did you call alloc_net_handle?'
365 0 : return
366 : end if
367 :
368 111936 : v(1:num_isos) = max(1d-30, min(1d0, v(1:num_isos)))
369 : ! positive definite mass fractions
370 223872 : v(1:num_isos) = v(1:num_isos)/sum(v(1:num_isos))
371 111936 : x(1:num_isos) = v(1:num_isos)
372 :
373 5088 : num_reactions = g% num_reactions
374 :
375 5088 : i = burn_const_P_lrpar
376 :
377 5088 : if (ipar(i_clip) /= 0) then
378 0 : do i=1,num_isos
379 0 : x(i) = max(0d0, min(1d0, x(i)))
380 : end do
381 : end if
382 :
383 : call basic_composition_info( &
384 5088 : num_isos, g% chem_id, x, xh, Y, z, abar, zbar, z2bar, z53bar, ye, mass_correction, sumx)
385 :
386 5088 : logT = v(nvar)/ln10
387 5088 : T = exp10(logT)
388 :
389 5088 : pressure = rpar(r_burn_const_P_pressure)
390 5088 : Prad = Radiation_Pressure(T)
391 5088 : Pgas = pressure - Prad
392 5088 : if (Pgas <= 0) then
393 0 : write(*,1) 'Pgas <= 0 in burn at const P', Pgas
394 0 : ierr = -1
395 0 : return
396 : end if
397 5088 : lgPgas = log10(Pgas)
398 :
399 5088 : chem_id => g% chem_id
400 5088 : net_iso => g% net_iso
401 :
402 5088 : if (rpar(r_burn_const_P_time_eos) >= 0) then
403 0 : call system_clock(time0,clock_rate)
404 : else
405 : time0 = 0
406 : end if
407 :
408 : call eosPT_get( &
409 : eos_handle, &
410 : num_isos, chem_id, net_iso, x, &
411 : Pgas, lgPgas, T, logT, &
412 : Rho, logRho, dlnRho_dlnPgas_const_T, dlnRho_dlnT_const_Pgas, &
413 : res, d_dlnRho_const_T, d_dlnT_const_Rho, &
414 5088 : d_dxa_const_TRho, ierr)
415 :
416 5088 : Cp = res(i_Cp)
417 5088 : eta = res(i_eta)
418 5088 : d_eta_dlnRho = d_dlnRho_const_T(i_eta)
419 5088 : d_eta_dlnT = d_dlnT_const_Rho(i_eta)
420 5088 : screening_mode = ipar(i_screening_mode)
421 5088 : rpar(r_burn_const_P_rho) = Rho
422 5088 : rpar(r_burn_const_P_temperature) = T
423 5088 : rpar(r_burn_const_P_lnS) = res(i_lnS)
424 5088 : if (rpar(r_burn_const_P_init_rho) < -1d90) &
425 2 : rpar(r_burn_const_P_init_rho) = Rho
426 5088 : if (rpar(r_burn_const_P_init_lnS) < -1d90) &
427 2 : rpar(r_burn_const_P_init_lnS) = res(i_lnS)
428 5088 : if (ierr /= 0 .or. Cp <= 0) then
429 0 : ierr = -1
430 0 : return
431 :
432 : write(*,*) 'eosPT_get failed'
433 : write(*,1) 'xh', xh
434 : write(*,1) 'Y', Y
435 : write(*,1) 'Z', 1 - (xh + Y)
436 : write(*,1) 'abar', abar
437 : write(*,1) 'zbar', zbar
438 : write(*,1) 'pressure', pressure
439 : write(*,1) 'Prad', Prad
440 : write(*,1) 'Pgas', Pgas
441 : write(*,1) 'lgPgas', lgPgas
442 : write(*,1) 'T', T
443 : write(*,1) 'logT', logT
444 : write(*,1) 'Rho', Rho
445 : write(*,1) 'logRho', logRho
446 : write(*,1) 'Cp', Cp
447 : return
448 : end if
449 :
450 5088 : if (rpar(r_burn_const_P_time_eos) >= 0) then
451 0 : call system_clock(time1,clock_rate)
452 : rpar(r_burn_const_P_time_eos) = &
453 0 : rpar(r_burn_const_P_time_eos) + dble(time1 - time0) / clock_rate
454 0 : if (rpar(r_burn_const_P_time_net) >= 0) time0 = time1
455 5088 : else if (rpar(r_burn_const_P_time_net) >= 0) then
456 0 : call system_clock(time0,clock_rate)
457 : end if
458 :
459 5088 : rates_only = .false.
460 5088 : skip_jacobian = .false.
461 :
462 : call eval_net( &
463 : n, g, rates_only, skip_jacobian, &
464 : num_isos, num_reactions, g% num_wk_reactions, &
465 : x, T, logT, rho, logRho, &
466 : abar, zbar, z2bar, ye, eta, d_eta_dlnT, d_eta_dlnRho, &
467 : rate_factors, weak_rate_factor, &
468 : reaction_Qs, reaction_neuQs, &
469 : eps_nuc, d_eps_nuc_dRho, d_eps_nuc_dT, d_eps_nuc_dx, &
470 : dxdt, d_dxdt_dRho, d_dxdt_dT, d_dxdt_dx, &
471 : screening_mode, &
472 : eps_nuc_categories, eps_neu_total, &
473 5088 : actual_Qs, actual_neuQs, from_weaklib, .false., ierr)
474 :
475 5088 : if (rpar(r_burn_const_P_time_net) >= 0) then
476 0 : call system_clock(time1,clock_rate)
477 : rpar(r_burn_const_P_time_net) = &
478 0 : rpar(r_burn_const_P_time_net) + dble(time1 - time0) / clock_rate
479 : end if
480 :
481 5088 : if (ierr /= 0) then
482 : return
483 :
484 :
485 :
486 : write(*,*) 'eval_net failed'
487 : write(*,1) 'xh', xh
488 : write(*,1) 'Y', Y
489 : write(*,1) 'Z', 1 - (xh + Y)
490 : write(*,1) 'abar', abar
491 : write(*,1) 'zbar', zbar
492 : write(*,1) 'pressure', pressure
493 : write(*,1) 'Prad', Prad
494 : write(*,1) 'Pgas', Pgas
495 : write(*,1) 'lgPgas', lgPgas
496 : write(*,1) 'T', T
497 : write(*,1) 'logT', logT
498 : write(*,1) 'Rho', Rho
499 : write(*,1) 'logRho', logRho
500 : write(*,1) 'Cp', Cp
501 : ierr = -1
502 :
503 :
504 : call mesa_error(__FILE__,__LINE__,'net_burn_const_P')
505 :
506 : return
507 : end if
508 :
509 111936 : f(1:num_isos) = dxdt
510 5088 : dlnT_dt = eps_nuc/(Cp*T)
511 5088 : f(nvar) = dlnT_dt
512 :
513 5088 : if (ld_dfdv > 0) then
514 :
515 842 : dlnRho_dlnT_const_P = -res(i_chiT)/res(i_chiRho)
516 842 : d_epsnuc_dlnT_const_P = d_eps_nuc_dT*T + d_eps_nuc_dRho*Rho*dlnRho_dlnT_const_P
517 842 : d_Cp_dlnT = d_dlnT_const_Rho(i_Cp) + d_dlnRho_const_T(i_Cp)*dlnRho_dlnT_const_P
518 :
519 389846 : dfdv(1:num_isos,1:num_isos) = d_dxdt_dx
520 :
521 842 : dfdv(nvar,nvar) = d_epsnuc_dlnT_const_P/(Cp*T) - dlnT_dt*(1 + d_Cp_dlnT/Cp)
522 :
523 : ! d_dxdt_dlnT
524 : dfdv(1:num_isos,nvar) = &
525 18524 : d_dxdt_dT(1:num_isos)*T + d_dxdt_dRho(1:num_isos)*Rho*dlnRho_dlnT_const_P
526 :
527 : ! d_dlnTdt_dx
528 18524 : dfdv(nvar,1:num_isos) = d_eps_nuc_dx(1:num_isos)/(Cp*T)
529 :
530 : end if
531 :
532 10176 : end subroutine burn_jacob
533 :
534 :
535 : subroutine burn_sjac(n,time,h,y,f,nzmax,ia,ja,values,lrpar,rpar,lipar,ipar,ierr)
536 5088 : use mtx_lib, only: dense_to_sparse_with_diag
537 : integer, intent(in) :: n, nzmax, lrpar, lipar
538 : real(dp), intent(in) :: time, h
539 : real(dp), intent(inout) :: y(n)
540 : integer, intent(out) :: ia(:), ja(:)
541 : real(dp), intent(inout) :: f(:), values(:)
542 : integer, intent(inout), pointer :: ipar(:) ! (lipar)
543 : real(dp), intent(inout), pointer :: rpar(:) ! (lrpar)
544 : integer, intent(out) :: ierr ! nonzero means terminate integration
545 : real(dp), pointer :: dfdv(:,:) ! (n,n)
546 : integer :: ld_dfdv, nz, i, j, cnt, nnz
547 : include 'formats'
548 : !write(*,1) 'burn_sjac', x
549 : ierr = 0
550 : ld_dfdv = n
551 : allocate(dfdv(n,n),stat=ierr)
552 : if(g% fill_arrays_with_nans) then
553 : call fill_with_NaNs_2D(dfdv)
554 : end if
555 : if (ierr /= 0) return
556 : call burn_jacob(n,time,h,y,f,dfdv,ld_dfdv,lrpar,rpar,lipar,ipar,ierr)
557 : if (ierr /= 0) then
558 : deallocate(dfdv)
559 : return
560 : end if
561 : ! remove entries with abs(value) < 1d-16
562 : cnt = 0; nnz = 0
563 : do i=1,n
564 : do j=1,n
565 : if (dfdv(i,j) /= 0) then
566 : nnz = nnz + 1
567 : if (abs(dfdv(i,j)) < 1d-16) then
568 : cnt = cnt+1; dfdv(i,j) = 0
569 : end if
570 : end if
571 : end do
572 : end do
573 : call dense_to_sparse_with_diag( &
574 : ipar(i_sparse_format),n,n,dfdv,nzmax,nz,ia,ja,values,ierr)
575 : deallocate(dfdv)
576 : end subroutine burn_sjac
577 :
578 :
579 : end subroutine burn_1_zone_const_P
580 :
581 :
582 : end module net_burn_const_P
583 :
|
