Line data Source code
1 : ! ***********************************************************************
2 : !
3 : ! Copyright (C) 2015 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_alpha_rti_eqns
21 :
22 : use star_private_def
23 : use const_def, only: dp
24 : use auto_diff_support
25 :
26 : implicit none
27 :
28 : private
29 : public :: do1_dalpha_RTI_dt_eqn
30 :
31 : contains
32 :
33 0 : subroutine do1_dalpha_RTI_dt_eqn(s, k, nvar, ierr)
34 : use star_utils, only: em1, e00, ep1
35 :
36 : type (star_info), pointer :: s
37 : integer, intent(in) :: k, nvar
38 : integer, intent(out) :: ierr
39 :
40 : type(auto_diff_real_4var_order1) :: a_m1, a_00, a_p1
41 : type(auto_diff_real_4var_order1) :: da00, flux00, dap1, fluxp1, RTI_D, A_plus_B_div_rho
42 : type(auto_diff_real_4var_order1) :: source_minus, source_plus, dadt_source, dadt_actual
43 : type(auto_diff_real_4var_order1) :: dadt_mix, dadt_expected, resid
44 :
45 : integer :: nz, i_alpha_RTI, i_dalpha_RTI_dt
46 0 : real(dp), pointer, dimension(:) :: sig
47 : real(dp) :: &
48 0 : dq, dm, sig00, sigp1, &
49 0 : eqn_scale, dPdr_drhodr, instability2, instability, RTI_B, &
50 0 : r00, rho, rmid, cs, fac
51 : logical :: test_partials
52 :
53 : include 'formats'
54 :
55 : ! Debugging setup
56 0 : ierr = 0
57 : !test_partials = (k == s% solver_test_partials_k)
58 0 : test_partials = .false.
59 :
60 : ! Equation setup
61 0 : i_alpha_RTI = s% i_alpha_RTI
62 0 : i_dalpha_RTI_dt = s% i_dalpha_RTI_dt
63 0 : nz = s% nz
64 :
65 : ! Starting/fixed quantities
66 0 : sig => s% sig_RTI
67 0 : dq = s% dq(k)
68 0 : dm = s% dm(k)
69 0 : rho = s% rho_start(k)
70 0 : r00 = s% r_start(k)
71 0 : fac = s% alpha_RTI_diffusion_factor
72 :
73 0 : sig00 = fac*sig(k)
74 0 : if (k < nz) then
75 0 : sigp1 = fac*sig(k+1)
76 : else
77 0 : sigp1 = 0
78 : end if
79 :
80 0 : a_00 = s% alpha_RTI(k)
81 0 : a_00%d1val2 = 1d0
82 :
83 : ! alpha's and fluxes
84 0 : if (k > 1) then
85 0 : a_m1 = s% alpha_RTI(k-1)
86 0 : a_m1%d1val1 = 1d0
87 :
88 0 : da00 = a_m1 - a_00
89 0 : flux00 = -sig00*da00
90 : else
91 0 : a_m1 = 0d0
92 0 : flux00 = 0d0
93 : end if
94 :
95 0 : if (k < nz) then
96 0 : a_p1 = s% alpha_RTI(k+1)
97 0 : a_p1%d1val3 = 1d0
98 :
99 0 : dap1 = a_00 - a_p1
100 0 : fluxp1 = -sigp1*dap1
101 : else
102 0 : a_p1 = 0d0
103 0 : fluxp1 = 0d0
104 : end if
105 :
106 : ! Flux divergence
107 0 : dadt_mix = (fluxp1 - flux00)/dm
108 :
109 : ! Sources and sink s
110 0 : dPdr_drhodr = s% dPdr_dRhodr_info(k)
111 :
112 0 : if (a_00 <= 0d0 .or. s% RTI_D <= 0d0) then
113 0 : source_minus = 0d0
114 : else
115 0 : cs = s% csound_start(k)
116 0 : if (k < nz) then
117 0 : rmid = 0.5d0*(s% r_start(k) + s% r_start(k+1))
118 : else
119 0 : rmid = 0.5d0*(s% r_start(k) + s% R_center)
120 : end if
121 0 : RTI_D = s% RTI_D*max(1d0,a_00/s% RTI_max_alpha)
122 0 : source_minus = RTI_D*a_00*cs/rmid
123 : end if
124 :
125 0 : instability2 = -dPdr_drhodr ! > 0 means Rayleigh-Taylor unstable
126 : if (instability2 <= 0d0 .or. &
127 0 : s% q(k) > s% alpha_RTI_src_max_q .or. &
128 : s% q(k) < s% alpha_RTI_src_min_q) then
129 0 : source_plus = 0d0
130 0 : instability2 = 0d0
131 : instability = 0d0
132 0 : A_plus_B_div_rho = 0d0
133 : else
134 0 : RTI_B = s% RTI_B
135 0 : instability = sqrt(instability2)
136 0 : if (s% alpha_RTI_start(k) < s% RTI_max_alpha) then
137 0 : A_plus_B_div_rho = (s% RTI_A + RTI_B*a_00)/rho
138 : else ! turn off source when reach max
139 0 : A_plus_B_div_rho = 0d0
140 : end if
141 0 : source_plus = A_plus_B_div_rho*instability
142 : end if
143 :
144 0 : s% source_minus_alpha_RTI(k) = source_minus%val
145 0 : s% source_plus_alpha_RTI(k) = source_plus%val
146 :
147 0 : dadt_source = source_plus - source_minus
148 :
149 0 : dadt_expected = dadt_mix + dadt_source
150 :
151 : ! We use dxh to avoid subtraction errors.
152 : ! At least that's what I assume. I just preserved this
153 : ! choice when re-writing it... - Adam S. Jermyn 6/15/2021
154 0 : dadt_actual = s% dxh_alpha_RTI(k)/s% dt
155 0 : dadt_actual%d1val2 = 1d0 / s% dt
156 :
157 0 : eqn_scale = max(1d0, s% x_scale(i_dalpha_RTI_dt,k)/s% dt)
158 0 : resid = (dadt_expected - dadt_actual)/eqn_scale
159 :
160 : ! Unpack
161 0 : s% equ(i_dalpha_RTI_dt,k) = resid%val
162 0 : call em1(s, i_dalpha_RTI_dt, i_alpha_RTI, k, nvar, resid%d1val1)
163 0 : call e00(s, i_dalpha_RTI_dt, i_alpha_RTI, k, nvar, resid%d1val2)
164 0 : call ep1(s, i_dalpha_RTI_dt, i_alpha_RTI, k, nvar, resid%d1val3)
165 :
166 : ! Debugging
167 : if (test_partials) then
168 : write(*,2) 'a_00', k, a_00%val
169 : write(*,2) 'dadt_mix', k, dadt_mix
170 : write(*,2) 'dadt_source', k, dadt_source
171 : write(*,2) 'source_plus', k, source_plus
172 : write(*,2) 'source_minus', k, source_minus
173 : write(*,2) 'A_plus_B_div_rho', k, A_plus_B_div_rho
174 : write(*,2) 'instability', k, instability
175 : write(*,2) 's% q(k)', k, s% q(k)
176 : write(*,2) 's% Peos(k)', k, s% Peos(k)
177 : write(*,2) 's% rho(k)', k, s% rho(k)
178 : write(*,2) 's% alpha_RTI_src_max_q', k, s% alpha_RTI_src_max_q
179 : write(*,2) 'dadt_expected', k, dadt_expected
180 : write(*,2) 'dadt_actual', k, dadt_actual
181 : write(*,2) 'eqn_scale', k, eqn_scale
182 : write(*,2) 's% dt', k, s% dt
183 : write(*,2) 'equ(i,k)', k, s% equ(i_dalpha_RTI_dt,k)
184 : s% solver_test_partials_val = s% equ(i_dalpha_RTI_dt,k)
185 : end if
186 :
187 : if (test_partials) then
188 : s% solver_test_partials_var = i_alpha_RTI
189 : s% solver_test_partials_dval_dx = resid%d1val2
190 : write(*,*) 'do1_dalpha_RTI_dt_eqn', s% solver_test_partials_var
191 : end if
192 :
193 0 : end subroutine do1_dalpha_RTI_dt_eqn
194 :
195 : end module hydro_alpha_rti_eqns
|
