Line data Source code
1 : ! ***********************************************************************
2 : !
3 : ! Copyright (C) 2012 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 hydro_chem_eqns
21 :
22 : use star_private_def
23 : use const_def, only: dp
24 : use utils_lib
25 :
26 : implicit none
27 :
28 : private
29 : public :: do_chem_eqns, do1_chem_eqns
30 :
31 : contains
32 :
33 0 : subroutine do_chem_eqns(s, nvar, ierr)
34 : type (star_info), pointer :: s
35 : integer, intent(in) :: nvar
36 : integer, intent(out) :: ierr
37 : integer :: k, op_err
38 : include 'formats'
39 0 : ierr = 0
40 0 : !$OMP PARALLEL DO PRIVATE(k,op_err) SCHEDULE(dynamic,2)
41 : do k = 1, s% nz
42 : if (ierr /= 0) cycle
43 : call do1_chem_eqns(s, k, nvar, op_err)
44 : if (op_err /= 0) ierr = op_err
45 : end do
46 : !$OMP END PARALLEL DO
47 0 : end subroutine do_chem_eqns
48 :
49 :
50 52348 : subroutine do1_chem_eqns(s, k, nvar, ierr)
51 :
52 : use chem_def
53 : use net_lib, only: show_net_reactions, show_net_params
54 : use rates_def, only: i_rate
55 : use star_utils, only: em1, e00, ep1
56 :
57 : type (star_info), pointer :: s
58 : integer, intent(in) :: k, nvar
59 : integer, intent(out) :: ierr
60 :
61 : !integer, pointer :: reaction_id(:) ! maps net reaction number to reaction id
62 : integer :: nz, j, i, jj, ii, species
63 : real(dp) :: &
64 52348 : dxdt_expected_dxa, dxdt_expected, dxdt_actual, &
65 52348 : dq, dm, dequ, dxdt_nuc, dxdt_mix, max_abs_residual, &
66 52348 : sum_dxdt_nuc, &
67 52348 : d_dxdt_mix_dx00, d_dxdt_mix_dxm1, d_dxdt_mix_dxp1, &
68 52348 : sum_dx_burning, sum_dx_mixing, residual, &
69 52348 : dxdt_factor, eqn_scale, &
70 52348 : dequ_dlnd, dequ_dlnT
71 : logical :: test_partials, doing_op_split_burn
72 : logical, parameter :: checking = .false.
73 :
74 : include 'formats'
75 :
76 52348 : ierr = 0
77 :
78 52348 : species = s% species
79 52348 : nz = s% nz
80 :
81 52348 : dq = s% dq(k)
82 52348 : dm = s% dm(k)
83 :
84 52348 : max_abs_residual = 0
85 52348 : sum_dxdt_nuc = 0
86 :
87 52348 : if (s% do_mix) then
88 :
89 52348 : d_dxdt_mix_dxm1 = s% d_dxdt_mix_dxm1(k)
90 52348 : d_dxdt_mix_dx00 = s% d_dxdt_mix_dx00(k)
91 52348 : d_dxdt_mix_dxp1 = s% d_dxdt_mix_dxp1(k)
92 :
93 : else
94 :
95 : d_dxdt_mix_dxm1 = 0
96 : d_dxdt_mix_dx00 = 0
97 : d_dxdt_mix_dxp1 = 0
98 :
99 : end if
100 :
101 52348 : sum_dx_burning = 0
102 52348 : sum_dx_mixing = 0
103 :
104 523480 : do j=1,species ! composition equation for species j in cell k
105 :
106 : !test_partials = (k == s% solver_test_partials_k .and. s% net_iso(ihe4) == j)
107 418784 : test_partials = .false.
108 :
109 418784 : i = s%nvar_hydro+j
110 :
111 418784 : dxdt_actual = s% xa_sub_xa_start(j,k)/s% dt
112 :
113 418784 : doing_op_split_burn = s% op_split_burn .and. s% T_start(k) >= s% op_split_burn_min_T
114 418784 : if (s% do_burn .and. .not. doing_op_split_burn) then
115 418784 : dxdt_nuc = s% dxdt_nuc(j,k)
116 : else
117 0 : dxdt_nuc = 0
118 : end if
119 :
120 418784 : if (s% do_mix) then
121 418784 : dxdt_mix = s% dxdt_mix(j,k)
122 : else
123 0 : dxdt_mix = 0
124 : end if
125 :
126 418784 : dxdt_expected = dxdt_mix + dxdt_nuc
127 :
128 418784 : dxdt_factor = 1d0
129 :
130 418784 : eqn_scale = max(s% min_chem_eqn_scale, s% x_scale(i,k)/s% dt)
131 418784 : residual = (dxdt_expected - dxdt_actual)/eqn_scale
132 418784 : s% equ(i,k) = residual
133 :
134 : if (abs(residual) > max_abs_residual) &
135 : max_abs_residual = abs(s% equ(i,k))
136 :
137 418784 : if (is_bad(s% equ(i,k))) then
138 0 : s% retry_message = 'bad residual for do1_chem_eqns'
139 0 : !$OMP critical (star_chem_eqns_bad_num)
140 0 : if (s% report_ierr) then
141 0 : write(*,3) 'do1_chem_eqns: equ ' // trim(s% nameofequ(i)), &
142 0 : i, k, s% equ(i,k)
143 0 : write(*,2) 'dxdt_expected', k, dxdt_expected
144 0 : write(*,2) 'dxdt_actual', k, dxdt_actual
145 0 : write(*,2) 'eqn_scale', k, eqn_scale
146 0 : write(*,2) 'dxdt_mix', k, dxdt_mix
147 0 : write(*,2) 'dxdt_nuc', k, dxdt_nuc
148 : end if
149 0 : if (s% stop_for_bad_nums) call mesa_error(__FILE__,__LINE__,'do1_chem_eqns')
150 : !$OMP end critical (star_chem_eqns_bad_num)
151 : end if
152 :
153 418784 : call e00(s, i, i, k, nvar, -1d0/eqn_scale/s% dt)
154 :
155 : ! all the rest are jacobian terms for dxdt_expected/eqn_scale
156 :
157 418784 : if (s% do_burn) then
158 :
159 3769056 : do jj=1,species
160 3350272 : ii = s% nvar_hydro+jj
161 3350272 : dxdt_expected_dxa = s% d_dxdt_nuc_dx(j,jj,k)
162 3350272 : dequ = dxdt_expected_dxa/eqn_scale
163 : if (checking) call check_dequ(dequ,'d_dxdt_nuc_dx')
164 3769056 : call e00(s, i, ii, k, nvar, dxdt_factor*dequ)
165 : end do
166 :
167 418784 : dequ_dlnd = s% d_dxdt_nuc_drho(j,k)*s% rho(k)/eqn_scale
168 418784 : call e00(s, i, s% i_lnd, k, nvar, dxdt_factor*dequ_dlnd)
169 :
170 418784 : dequ_dlnT = s% d_dxdt_nuc_dT(j,k)*s% T(k)/eqn_scale
171 418784 : call e00(s, i, s% i_lnT, k, nvar, dxdt_factor*dequ_dlnT)
172 :
173 : end if
174 :
175 418784 : if (s% do_mix) then
176 :
177 418784 : dxdt_expected_dxa = d_dxdt_mix_dx00
178 418784 : dequ = dxdt_expected_dxa/eqn_scale
179 : if (checking) call check_dequ(dequ,'d_dxdt_mix_dx00')
180 418784 : call e00(s, i, i, k, nvar, dxdt_factor*dequ)
181 418784 : if (k > 1) then
182 418432 : dxdt_expected_dxa = d_dxdt_mix_dxm1
183 418432 : dequ = dxdt_expected_dxa/eqn_scale
184 : if (checking) call check_dequ(dequ,'d_dxdt_mix_dxm1')
185 418432 : call em1(s, i, i, k, nvar, dxdt_factor*dequ)
186 : end if
187 418784 : if (k < nz) then
188 418432 : dxdt_expected_dxa = d_dxdt_mix_dxp1
189 418432 : dequ = dxdt_expected_dxa/eqn_scale
190 : if (checking) call check_dequ(dequ,'d_dxdt_mix_dxp1')
191 418432 : call ep1(s, i, i, k, nvar, dxdt_factor*dequ)
192 : end if
193 :
194 : end if
195 :
196 52348 : if (test_partials) then
197 : s% solver_test_partials_dx_sink = s% net_iso(img24)
198 : s% solver_test_partials_val = s% dxdt_nuc(j,k)
199 : s% solver_test_partials_var = s% nvar_hydro + j
200 : s% solver_test_partials_dval_dx = s% d_dxdt_nuc_dx(j,j,k)
201 : write(*,*) 'do1_chem_eqns', s% solver_test_partials_var
202 : end if
203 :
204 : end do
205 :
206 : contains
207 :
208 : subroutine check_dequ(dequ, str)
209 : real(dp), intent(in) :: dequ
210 : character (len=*), intent(in) :: str
211 : include 'formats'
212 : if (is_bad(dequ)) then
213 : ierr = -1
214 : if (s% report_ierr) then
215 : write(*,2) 'do1_chem_eqns: bad ' // trim(str), k
216 : end if
217 : if (s% stop_for_bad_nums) call mesa_error(__FILE__,__LINE__,'do1_chem_eqns')
218 : return
219 : end if
220 52348 : end subroutine check_dequ
221 :
222 : end subroutine do1_chem_eqns
223 :
224 : end module hydro_chem_eqns
|
