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 1 : 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 : 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 : real(dp) :: mixing_length_alpha, conv_vel_start, &
84 : TDC_alpha_D, TDC_alpha_R, TDC_alpha_Pt, dt, cgrav, m, scale, L_start, TDC_alpha_C, TDC_alpha_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, energy
88 : 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, use_TDC_enthalpy_flux_limiter
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 : TDC_alpha_D = 1.0d0
128 1 : TDC_alpha_R = 0.0d0
129 1 : TDC_alpha_Pt = 0.0d0
130 1 : TDC_alpha_C = 1.0d0
131 1 : TDC_alpha_S = 1.0d0
132 1 : dV = 0d0
133 1 : energy = 0d0
134 1 : conv_vel_start = 0d0 !1d10
135 1 : scale = L%val*1d-3
136 1 : report = .false.
137 1 : dt = 1d40 ! Long time-step so we get into equilibrium
138 1 : Eq_div_w = 0d0 ! TDC_alpha_M is implicit in this term
139 1 : include_mlt_corr_to_TDC = .true.
140 1 : use_TDC_enthalpy_flux_limiter = .false.
141 :
142 : ! MLT
143 1 : MLT_option = 'Cox'
144 1 : Henyey_MLT_nu_param = 0d0
145 1 : Henyey_MLT_y_param = 0d0
146 :
147 1 : write (*, 1) 'gradR - gradA', gradr%val - grada%val
148 :
149 : call set_TDC( &
150 : conv_vel_start, mixing_length_alpha, TDC_alpha_D, TDC_alpha_R, TDC_alpha_Pt, dt, cgrav, m, report, &
151 : mixing_type, scale, chiT, chiRho, gradr, r, P, T, rho, dV, Cp, opacity, &
152 : scale_height, gradL, grada, conv_vel, D, Y_face, gradT, tdc_num_iters, max_conv_vel, &
153 : Eq_div_w, grav, include_mlt_corr_to_TDC, TDC_alpha_C, TDC_alpha_S, use_TDC_enthalpy_flux_limiter, &
154 1 : energy, ierr)
155 :
156 :
157 1 : write (*, 1) 'TDC: Y, conv_vel_start, conv_vel, dt ', Y_face%val, conv_vel_start, conv_vel%val, dt
158 :
159 : call set_MLT(MLT_option, mixing_length_alpha, Henyey_MLT_nu_param, Henyey_MLT_y_param, &
160 : chiT, chiRho, Cp, grav, Lambda, rho, P, T, opacity, &
161 : gradr, grada, gradL, &
162 1 : Gamma, gradT, Y_face, conv_vel, D, mixing_type, max_conv_vel, ierr)
163 :
164 1 : write (*, 1) 'MLT: Y, conv_vel_start, conv_vel, Gamma', Y_face%val, conv_vel_start, conv_vel%val, Gamma%val
165 :
166 1 : end subroutine compare_TDC_and_Cox_MLT
167 :
168 1 : subroutine check_TDC()
169 : real(dp) :: mixing_length_alpha, conv_vel_start
170 : real(dp) :: TDC_alpha_D, TDC_alpha_R, TDC_alpha_Pt, dt, cgrav, m, scale, max_conv_vel, L_start, TDC_alpha_C, TDC_alpha_S
171 : type(auto_diff_real_star_order1) :: &
172 : r, L, T, P, opacity, rho, dV, chiRho, chiT, Cp, gradr, grada, scale_height, gradL
173 : type(auto_diff_real_star_order1) :: gradT, Y_face, conv_vel, D, Eq_div_w, grav, energy
174 : integer :: mixing_type, ierr, tdc_num_iters
175 : logical :: report, include_mlt_corr_to_TDC, use_TDC_enthalpy_flux_limiter
176 : integer :: j
177 :
178 : include 'formats'
179 :
180 1 : call header('Test TDC')
181 :
182 : ! For limiting the conv_vel coming out of mlt/TDC with Csound.
183 1 : max_conv_vel = 1d99 ! we don't limit the conv_vel for testing.
184 :
185 1 : conv_vel_start = 52320587.415154047d0
186 :
187 1 : mixing_length_alpha = 2.0d0
188 1 : TDC_alpha_D = 1.0d0
189 1 : TDC_alpha_R = 0.0d0
190 1 : TDC_alpha_Pt = 0.0d0
191 1 : TDC_alpha_C = 1.0d0
192 1 : TDC_alpha_S = 1.0d0
193 1 : cgrav = 6.6743000000000004d-8
194 1 : m = 5.8707400456875664d34
195 1 : scale = 5.0386519362246294d45
196 1 : L = 1.0941528815883500015d0*(-5.0386519362246294d45)
197 1 : r = 10314294541.567163d0
198 1 : P = 5.0581587249808894d20
199 1 : T = 613193666.51783681d0
200 1 : rho = 5204.5732574745753d0
201 1 : dV = 3.8256494463482604d-7
202 1 : Cp = 6628075118.4606590d0
203 1 : opacity = 9.0750171231469945d-2
204 1 : scale_height = 2638686602.0063782d0
205 1 : gradL = 0.25207587267343501d0
206 1 : grada = 0.25204697256872738d0
207 1 : report = .false.
208 1 : chiT = 1d0
209 1 : chiRho = 1d0
210 :
211 1 : gradr = 3d0 * P * opacity * L / (64 * pi * boltz_sigma * pow4(T) * cgrav * m)
212 1 : grav = m * cgrav / pow2(r)
213 1 : Eq_div_w = 0d0 ! TDC_alpha_M is implicit in this term
214 1 : energy = 0d0
215 1 : include_mlt_corr_to_TDC = .true.
216 1 : use_TDC_enthalpy_flux_limiter = .false.
217 :
218 1 : write (*, *) "####################################"
219 1 : write (*, *) "Running dt test"
220 :
221 32 : do j = 0, 30
222 31 : dt = 500d0*pow(1.02d0, j)
223 : call set_TDC( &
224 : conv_vel_start, mixing_length_alpha, TDC_alpha_D, TDC_alpha_R, TDC_alpha_Pt, dt, cgrav, m, report, &
225 : mixing_type, scale, chiT, chiRho, gradr, r, P, T, rho, dV, Cp, opacity, &
226 : scale_height, gradL, grada, conv_vel, D, Y_face, gradT, tdc_num_iters, max_conv_vel, &
227 : Eq_div_w, grav, include_mlt_corr_to_TDC, TDC_alpha_C, TDC_alpha_S, use_TDC_enthalpy_flux_limiter, &
228 31 : energy, ierr)
229 :
230 :
231 31 : write (*, 1) 'dt, gradT, conv_vel_start, conv_vel', dt, gradT%val, conv_vel_start, conv_vel%val
232 32 : if (report) stop
233 : end do
234 :
235 133 : end subroutine check_TDC
236 :
237 : end program test_turb
|
