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