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, alpha_TDC_DAMP, alpha_TDC_DAMPR, alpha_TDC_PtdVdt, dt, cgrav, m, scale
84 : type(auto_diff_real_star_order1) :: &
85 : r, L, T, P, opacity, rho, dV, chiRho, chiT, Cp, gradr, grada, scale_height, gradL, grav, Lambda
86 : type(auto_diff_real_star_order1) :: gradT, Y_face, conv_vel, D, Gamma
87 1 : real(dp) :: Henyey_MLT_nu_param, Henyey_MLT_y_param, max_conv_vel
88 : character(len=3) :: MLT_option
89 : integer :: mixing_type, ierr, tdc_num_iters
90 : logical :: report
91 :
92 : include 'formats'
93 :
94 1 : call header('Compare TDC with MLT Cox')
95 :
96 : ! For limiting the conv_vel coming out of mlt/TDC with Csound.
97 1 : max_conv_vel = 1d99 ! we don't limit the conv_vel for testing.
98 :
99 : ! General
100 1 : mixing_length_alpha = 2.0d0
101 1 : chiT = 1d0
102 1 : chiRho = 1d0
103 1 : T = 1d5
104 1 : rho = 1d-5
105 1 : r = Rsun
106 1 : m = Msun
107 1 : cgrav = standard_cgrav
108 1 : grav = m*cgrav/pow2(r)
109 1 : Cp = 2.5d0*kerg/mp
110 1 : P = rho*T*kerg/mp
111 1 : scale_height = P/(rho*grav)
112 1 : Lambda = mixing_length_alpha*scale_height
113 1 : opacity = 1d0
114 1 : grada = 0.4d0
115 1 : gradL = grada
116 :
117 1 : L = 70*Lsun
118 1 : gradr = 3d0*P*opacity*L/(64*pi*boltz_sigma*pow4(T)*grav*pow2(r))
119 :
120 : ! Adjust L down to get just slightly superadiabatic gradR)
121 1 : L = L*(1d0 + 1d-5)*(grada/gradr)
122 1 : gradr = 3d0*P*opacity*L/(64*pi*boltz_sigma*pow4(T)*grav*pow2(r))
123 :
124 : ! TDC
125 1 : alpha_TDC_DAMP = 1.0d0
126 1 : alpha_TDC_DAMPR = 0.0d0
127 1 : alpha_TDC_PtdVdt = 0.0d0
128 1 : dV = 0d0
129 1 : conv_vel_start = 0d0 !1d10
130 1 : scale = L%val*1d-3
131 1 : report = .false.
132 1 : dt = 1d40 ! Long time-step so we get into equilibrium
133 :
134 : ! MLT
135 1 : MLT_option = 'Cox'
136 1 : Henyey_MLT_nu_param = 0d0
137 1 : Henyey_MLT_y_param = 0d0
138 :
139 1 : write (*, 1) 'gradR - gradA', gradr%val - grada%val
140 :
141 : call set_TDC( &
142 : conv_vel_start, mixing_length_alpha, alpha_TDC_DAMP, alpha_TDC_DAMPR, alpha_TDC_PtdVdt, dt, cgrav, m, report, &
143 : mixing_type, scale, chiT, chiRho, gradr, r, P, T, rho, dV, Cp, opacity, &
144 1 : scale_height, gradL, grada, conv_vel, D, Y_face, gradT, tdc_num_iters, max_conv_vel, ierr)
145 :
146 1 : write (*, 1) 'TDC: Y, conv_vel_start, conv_vel, dt ', Y_face%val, conv_vel_start, conv_vel%val, dt
147 :
148 : call set_MLT(MLT_option, mixing_length_alpha, Henyey_MLT_nu_param, Henyey_MLT_y_param, &
149 : chiT, chiRho, Cp, grav, Lambda, rho, P, T, opacity, &
150 : gradr, grada, gradL, &
151 1 : Gamma, gradT, Y_face, conv_vel, D, mixing_type, max_conv_vel, ierr)
152 :
153 1 : write (*, 1) 'MLT: Y, conv_vel_start, conv_vel, Gamma', Y_face%val, conv_vel_start, conv_vel%val, Gamma%val
154 :
155 1 : end subroutine compare_TDC_and_Cox_MLT
156 :
157 1 : subroutine check_TDC()
158 1 : real(dp) :: mixing_length_alpha, conv_vel_start
159 1 : real(dp) :: alpha_TDC_DAMP, alpha_TDC_DAMPR, alpha_TDC_PtdVdt, dt, cgrav, m, scale, max_conv_vel
160 : type(auto_diff_real_star_order1) :: &
161 : r, L, T, P, opacity, rho, dV, chiRho, chiT, Cp, gradr, grada, scale_height, gradL
162 : type(auto_diff_real_star_order1) :: gradT, Y_face, conv_vel, D
163 : integer :: mixing_type, ierr, tdc_num_iters
164 : logical :: report
165 : integer :: j
166 :
167 : include 'formats'
168 :
169 1 : call header('Test TDC')
170 :
171 : ! For limiting the conv_vel coming out of mlt/TDC with Csound.
172 1 : max_conv_vel = 1d99 ! we don't limit the conv_vel for testing.
173 :
174 1 : conv_vel_start = 52320587.415154047d0
175 :
176 1 : mixing_length_alpha = 2.0d0
177 1 : alpha_TDC_DAMP = 1.0d0
178 1 : alpha_TDC_DAMPR = 0.0d0
179 1 : alpha_TDC_PtdVdt = 0.0d0
180 1 : cgrav = 6.6743000000000004d-8
181 1 : m = 5.8707400456875664d34
182 1 : scale = 5.0386519362246294d45
183 1 : L = 1.0941528815883500015d0*(-5.0386519362246294d45)
184 1 : r = 10314294541.567163d0
185 1 : P = 5.0581587249808894d20
186 1 : T = 613193666.51783681d0
187 1 : rho = 5204.5732574745753d0
188 1 : dV = 3.8256494463482604d-7
189 1 : Cp = 6628075118.4606590d0
190 1 : opacity = 9.0750171231469945d-2
191 1 : scale_height = 2638686602.0063782d0
192 1 : gradL = 0.25207587267343501d0
193 1 : grada = 0.25204697256872738d0
194 1 : report = .false.
195 1 : chiT = 1d0
196 1 : chiRho = 1d0
197 1 : gradr = 3d0*P*opacity*L/(64*pi*boltz_sigma*pow4(T)*cgrav*m)
198 :
199 1 : write (*, *) "####################################"
200 1 : write (*, *) "Running dt test"
201 :
202 32 : do j = 0, 30
203 31 : dt = 500d0*pow(1.02d0, j)
204 : call set_TDC( &
205 : conv_vel_start, mixing_length_alpha, alpha_TDC_DAMP, alpha_TDC_DAMPR, alpha_TDC_PtdVdt, dt, cgrav, m, report, &
206 : mixing_type, scale, chiT, chiRho, gradr, r, P, T, rho, dV, Cp, opacity, &
207 31 : scale_height, gradL, grada, conv_vel, D, Y_face, gradT, tdc_num_iters, max_conv_vel, ierr)
208 :
209 31 : write (*, 1) 'dt, gradT, conv_vel_start, conv_vel', dt, gradT%val, conv_vel_start, conv_vel%val
210 63 : if (report) stop
211 : end do
212 :
213 1 : end subroutine check_TDC
214 :
215 : end program test_turb
|
