Line data Source code
1 1 : program test_turb
2 :
3 1 : use math_lib
4 : use auto_diff
5 : use const_def, only: dp, pi, rsun, lsun, msun, kerg, mp, boltz_sigma, standard_cgrav
6 : use turb
7 :
8 : implicit none
9 :
10 1 : call check_efficient_MLT_scaling()
11 1 : call check_TDC()
12 1 : call compare_TDC_and_Cox_MLT()
13 :
14 : contains
15 :
16 3 : subroutine header(text)
17 : character(len=*), intent(in) :: text
18 :
19 3 : write (*, '(a)') ' ----------------------------------------------------------------'
20 3 : write (*, '(a)') ' '
21 3 : write (*, '(a)') ' '//text
22 3 : write (*, '(a)') ' '
23 3 : write (*, '(a)') ' ----------------------------------------------------------------'
24 :
25 4 : end subroutine header
26 :
27 4 : subroutine check_efficient_MLT_scaling()
28 : type(auto_diff_real_star_order1) :: chiT, chiRho, Cp, grav, Lambda, rho, P, T, opacity, gradr, grada, gradL
29 : character(len=3) :: MLT_option
30 1 : real(dp) :: mixing_length_alpha, Henyey_MLT_nu_param, Henyey_MLT_y_param, max_conv_vel
31 : type(auto_diff_real_star_order1) :: Gamma, gradT, Y_face, conv_vel, conv_vel2, D, r, L
32 : integer :: mixing_type, ierr
33 :
34 : include 'formats'
35 :
36 1 : call header('Test MLT')
37 :
38 1 : MLT_option = 'Cox'
39 1 : mixing_length_alpha = 1d0
40 1 : Henyey_MLT_nu_param = 0d0
41 1 : Henyey_MLT_y_param = 0d0
42 1 : chiT = 1d0
43 1 : chiRho = 1d0
44 :
45 1 : T = 1d5
46 1 : rho = 1d-5
47 1 : grav = 1d4
48 1 : r = Rsun
49 1 : Cp = 2.5d0*kerg/mp
50 1 : P = rho*T*kerg/mp
51 1 : Lambda = P/(rho*grav)
52 1 : opacity = 1d0
53 1 : grada = 0.4d0
54 1 : gradL = grada
55 :
56 1 : L = 1d5*Lsun
57 1 : gradr = 3d0*P*opacity*L/(64*pi*boltz_sigma*pow4(T)*grav*pow2(r))
58 :
59 1 : max_conv_vel = 1d99
60 :
61 : call set_MLT(MLT_option, mixing_length_alpha, Henyey_MLT_nu_param, Henyey_MLT_y_param, &
62 : chiT, chiRho, Cp, grav, Lambda, rho, P, T, opacity, &
63 : gradr, grada, gradL, &
64 1 : Gamma, gradT, Y_face, conv_vel, D, mixing_type, max_conv_vel, ierr)
65 :
66 1 : write (*, 1) 'vc at 1d5 Lsun', conv_vel%val
67 :
68 1 : L = 1d8*Lsun
69 1 : gradr = 3d0*P*opacity*L/(64*pi*boltz_sigma*pow4(T)*grav*pow2(r))
70 :
71 : call set_MLT(MLT_option, mixing_length_alpha, Henyey_MLT_nu_param, Henyey_MLT_y_param, &
72 : chiT, chiRho, Cp, grav, Lambda, rho, P, T, opacity, &
73 : gradr, grada, gradL, &
74 1 : Gamma, gradT, Y_face, conv_vel2, D, mixing_type, max_conv_vel, ierr)
75 :
76 1 : write (*, 1) 'vc at 1d8 Lsun', conv_vel2%val
77 :
78 1 : write (*, 1) 'Ratio:', conv_vel2%val/conv_vel%val
79 1 : write (*, '(a)') 'Expected ~10 because in the efficient limit vc ~ L^{1/3}'
80 1 : end subroutine check_efficient_MLT_scaling
81 :
82 4 : subroutine compare_TDC_and_Cox_MLT()
83 1 : real(dp) :: mixing_length_alpha, conv_vel_start, &
84 1 : alpha_TDC_DAMP, alpha_TDC_DAMPR, alpha_TDC_PtdVdt, dt, cgrav, m, scale, L_start, alpha_TDC_C, alpha_TDC_S
85 : type(auto_diff_real_star_order1) :: &
86 : r, L, T, P, opacity, rho, dV, chiRho, chiT, Cp, gradr, grada, scale_height, gradL, grav, Lambda
87 : type(auto_diff_real_star_order1) :: gradT, Y_face, conv_vel, D, Gamma, Eq_div_w
88 1 : real(dp) :: Henyey_MLT_nu_param, Henyey_MLT_y_param, max_conv_vel
89 :
90 : character(len=3) :: MLT_option
91 : integer :: mixing_type, ierr, tdc_num_iters
92 : logical :: report, include_mlt_corr_to_TDC
93 :
94 : include 'formats'
95 :
96 1 : call header('Compare TDC with MLT Cox')
97 :
98 : ! For limiting the conv_vel coming out of mlt/TDC with Csound.
99 1 : max_conv_vel = 1d99 ! we don't limit the conv_vel for testing.
100 :
101 : ! General
102 1 : mixing_length_alpha = 2.0d0
103 1 : chiT = 1d0
104 1 : chiRho = 1d0
105 1 : T = 1d5
106 1 : rho = 1d-5
107 1 : r = Rsun
108 1 : m = Msun
109 1 : cgrav = standard_cgrav
110 1 : grav = m*cgrav/pow2(r)
111 1 : Cp = 2.5d0*kerg/mp
112 1 : P = rho*T*kerg/mp
113 1 : scale_height = P/(rho*grav)
114 1 : Lambda = mixing_length_alpha*scale_height
115 1 : opacity = 1d0
116 1 : grada = 0.4d0
117 1 : gradL = grada
118 :
119 1 : L = 70*Lsun
120 1 : gradr = 3d0*P*opacity*L/(64*pi*boltz_sigma*pow4(T)*grav*pow2(r))
121 :
122 : ! Adjust L down to get just slightly superadiabatic gradR)
123 1 : L = L*(1d0 + 1d-5)*(grada/gradr)
124 1 : gradr = 3d0*P*opacity*L/(64*pi*boltz_sigma*pow4(T)*grav*pow2(r))
125 :
126 : ! TDC
127 1 : alpha_TDC_DAMP = 1.0d0
128 1 : alpha_TDC_DAMPR = 0.0d0
129 1 : alpha_TDC_PtdVdt = 0.0d0
130 1 : alpha_TDC_C = 1.0d0
131 1 : alpha_TDC_S = 1.0d0
132 1 : dV = 0d0
133 1 : conv_vel_start = 0d0 !1d10
134 1 : scale = L%val*1d-3
135 1 : report = .false.
136 1 : dt = 1d40 ! Long time-step so we get into equilibrium
137 1 : Eq_div_w = 0d0
138 1 : include_mlt_corr_to_TDC = .true.
139 :
140 : ! MLT
141 1 : MLT_option = 'Cox'
142 1 : Henyey_MLT_nu_param = 0d0
143 1 : Henyey_MLT_y_param = 0d0
144 :
145 1 : write (*, 1) 'gradR - gradA', gradr%val - grada%val
146 :
147 : call set_TDC( &
148 : conv_vel_start, mixing_length_alpha, alpha_TDC_DAMP, alpha_TDC_DAMPR, alpha_TDC_PtdVdt, dt, cgrav, m, report, &
149 : mixing_type, scale, chiT, chiRho, gradr, r, P, T, rho, dV, Cp, opacity, &
150 : scale_height, gradL, grada, conv_vel, D, Y_face, gradT, tdc_num_iters, max_conv_vel, &
151 1 : Eq_div_w, grav, include_mlt_corr_to_TDC, alpha_TDC_C, alpha_TDC_S, ierr)
152 :
153 :
154 1 : write (*, 1) 'TDC: Y, conv_vel_start, conv_vel, dt ', Y_face%val, conv_vel_start, conv_vel%val, dt
155 :
156 : call set_MLT(MLT_option, mixing_length_alpha, Henyey_MLT_nu_param, Henyey_MLT_y_param, &
157 : chiT, chiRho, Cp, grav, Lambda, rho, P, T, opacity, &
158 : gradr, grada, gradL, &
159 1 : Gamma, gradT, Y_face, conv_vel, D, mixing_type, max_conv_vel, ierr)
160 :
161 1 : write (*, 1) 'MLT: Y, conv_vel_start, conv_vel, Gamma', Y_face%val, conv_vel_start, conv_vel%val, Gamma%val
162 :
163 1 : end subroutine compare_TDC_and_Cox_MLT
164 :
165 1 : subroutine check_TDC()
166 1 : real(dp) :: mixing_length_alpha, conv_vel_start
167 1 : real(dp) :: alpha_TDC_DAMP, alpha_TDC_DAMPR, alpha_TDC_PtdVdt, dt, cgrav, m, scale, max_conv_vel, L_start, &
168 1 : alpha_TDC_C, alpha_TDC_S
169 : type(auto_diff_real_star_order1) :: &
170 : r, L, T, P, opacity, rho, dV, chiRho, chiT, Cp, gradr, grada, scale_height, gradL
171 : type(auto_diff_real_star_order1) :: gradT, Y_face, conv_vel, D, Eq_div_w, grav
172 : integer :: mixing_type, ierr, tdc_num_iters
173 : logical :: report, include_mlt_corr_to_TDC
174 : integer :: j
175 :
176 : include 'formats'
177 :
178 1 : call header('Test TDC')
179 :
180 : ! For limiting the conv_vel coming out of mlt/TDC with Csound.
181 1 : max_conv_vel = 1d99 ! we don't limit the conv_vel for testing.
182 :
183 1 : conv_vel_start = 52320587.415154047d0
184 :
185 1 : mixing_length_alpha = 2.0d0
186 1 : alpha_TDC_DAMP = 1.0d0
187 1 : alpha_TDC_DAMPR = 0.0d0
188 1 : alpha_TDC_PtdVdt = 0.0d0
189 1 : alpha_TDC_C = 1.0d0
190 1 : alpha_TDC_S = 1.0d0
191 1 : cgrav = 6.6743000000000004d-8
192 1 : m = 5.8707400456875664d34
193 1 : scale = 5.0386519362246294d45
194 1 : L = 1.0941528815883500015d0*(-5.0386519362246294d45)
195 1 : r = 10314294541.567163d0
196 1 : P = 5.0581587249808894d20
197 1 : T = 613193666.51783681d0
198 1 : rho = 5204.5732574745753d0
199 1 : dV = 3.8256494463482604d-7
200 1 : Cp = 6628075118.4606590d0
201 1 : opacity = 9.0750171231469945d-2
202 1 : scale_height = 2638686602.0063782d0
203 1 : gradL = 0.25207587267343501d0
204 1 : grada = 0.25204697256872738d0
205 1 : report = .false.
206 1 : chiT = 1d0
207 1 : chiRho = 1d0
208 :
209 1 : gradr = 3d0 * P * opacity * L / (64 * pi * boltz_sigma * pow4(T) * cgrav * m)
210 1 : grav = m * cgrav / pow2(r)
211 1 : Eq_div_w = 0d0
212 1 : include_mlt_corr_to_TDC = .true.
213 :
214 1 : write (*, *) "####################################"
215 1 : write (*, *) "Running dt test"
216 :
217 32 : do j = 0, 30
218 31 : dt = 500d0*pow(1.02d0, j)
219 : call set_TDC( &
220 : conv_vel_start, mixing_length_alpha, alpha_TDC_DAMP, alpha_TDC_DAMPR, alpha_TDC_PtdVdt, dt, cgrav, m, report, &
221 : mixing_type, scale, chiT, chiRho, gradr, r, P, T, rho, dV, Cp, opacity, &
222 : scale_height, gradL, grada, conv_vel, D, Y_face, gradT, tdc_num_iters, max_conv_vel, &
223 31 : Eq_div_w, grav, include_mlt_corr_to_TDC, alpha_TDC_C, alpha_TDC_S, ierr)
224 :
225 :
226 31 : write (*, 1) 'dt, gradT, conv_vel_start, conv_vel', dt, gradT%val, conv_vel_start, conv_vel%val
227 63 : if (report) stop
228 : end do
229 :
230 1 : end subroutine check_TDC
231 :
232 : end program test_turb
|
