diff --git a/star/private/hydro_energy.f90 b/star/private/hydro_energy.f90 index 41c1c8525..81f2c745b 100644 --- a/star/private/hydro_energy.f90 +++ b/star/private/hydro_energy.f90 @@ -221,7 +221,7 @@ subroutine setup_dL_dm(ierr) end subroutine setup_dL_dm subroutine setup_sources_and_others(ierr) ! sources_ad, others_ad - !use hydro_rsp2, only: compute_Eq_cell + use hydro_rsp2, only: compute_Eq_cell integer, intent(out) :: ierr type(auto_diff_real_star_order1) :: & eps_nuc_ad, non_nuc_neu_ad, extra_heat_ad, Eq_ad, RTI_diffusion_ad, & @@ -264,12 +264,9 @@ subroutine setup_sources_and_others(ierr) ! sources_ad, others_ad others_ad%val = others_ad%val + s% eps_pre_mix(k) if (s% do_phase_separation .and. s% do_phase_separation_heating) & others_ad%val = others_ad%val + s% eps_phase_separation(k) - - Eq_ad = 0d0 - if (s% RSP2_flag) then - Eq_ad = s% Eq_ad(k) ! compute_Eq_cell(s, k, ierr) - if (ierr /= 0) return - end if + + Eq_ad = compute_Eq_cell(s, k, ierr) ! s% Eq_ad(k) XXX + if (ierr /= 0) return call setup_RTI_diffusion(RTI_diffusion_ad) diff --git a/star/private/hydro_momentum.f90 b/star/private/hydro_momentum.f90 index 54b29fddc..6c753dbc8 100644 --- a/star/private/hydro_momentum.f90 +++ b/star/private/hydro_momentum.f90 @@ -434,11 +434,8 @@ subroutine expected_non_HSE_term( & end if ! v_flag - Uq_ad = 0d0 - if (s% RSP2_flag) then ! Uq(k) is turbulent viscosity drag at face k - Uq_ad = compute_Uq_face(s, k, ierr) - if (ierr /= 0) return - end if + Uq_ad = compute_Uq_face(s, k, ierr) + if (ierr /= 0) return other_ad = extra_ad - accel_ad + drag + Uq_ad other = other_ad%val diff --git a/star/private/hydro_riemann.f90 b/star/private/hydro_riemann.f90 index 3fc78617c..a7f869982 100644 --- a/star/private/hydro_riemann.f90 +++ b/star/private/hydro_riemann.f90 @@ -410,11 +410,10 @@ subroutine do1_uface_and_Pface(s, k, ierr) end if end if - if (s% RSP2_flag) then ! include Uq in u_face - Uq_ad = compute_Uq_face(s, k, ierr) - if (ierr /= 0) return - s% u_face_ad(k) = s% u_face_ad(k) + Uq_ad - end if + ! include Uq in u_face + Uq_ad = compute_Uq_face(s, k, ierr) + if (ierr /= 0) return + s% u_face_ad(k) = s% u_face_ad(k) + Uq_ad s% u_face_val(k) = s% u_face_ad(k)%val diff --git a/star/private/hydro_rsp2.f90 b/star/private/hydro_rsp2.f90 index 3f040ebea..fec0542cd 100644 --- a/star/private/hydro_rsp2.f90 +++ b/star/private/hydro_rsp2.f90 @@ -664,7 +664,7 @@ function compute_Chi_cell(s, k, ierr) result(Chi_cell) real(dp) :: f, ALFAM_ALFA include 'formats' ierr = 0 - ALFAM_ALFA = s% RSP2_alfam*s% mixing_length_alpha + ALFAM_ALFA = s% RSP2_alfam * s% mixing_length_alpha if (ALFAM_ALFA == 0d0 .or. & k <= s% RSP2_num_outermost_cells_forced_nonturbulent .or. & k > s% nz - int(s% nz/s% RSP2_nz_div_IBOTOM)) then @@ -735,10 +735,10 @@ function compute_Uq_face(s, k, ierr) result(Uq_face) ! cm s^-2, acceleration Uq_face = 0d0 else r_00 = wrap_opt_time_center_r_00(s,k) - Chi_00 = s% Chi_ad(k) ! compute_Chi_cell(s,k,ierr) + Chi_00 = compute_Chi_cell(s,k,ierr) ! s% Chi_ad(k) XXX if (k > 1) then - !Chi_m1 = shift_m1(compute_Chi_cell(s,k-1,ierr)) - Chi_m1 = shift_m1(s% Chi_ad(k-1)) + Chi_m1 = shift_m1(compute_Chi_cell(s,k-1,ierr)) + !Chi_m1 = shift_m1(s% Chi_ad(k-1)) XXX if (ierr /= 0) return else Chi_m1 = 0d0 diff --git a/star/private/turb_support.f90 b/star/private/turb_support.f90 index c99508b69..06bd96634 100644 --- a/star/private/turb_support.f90 +++ b/star/private/turb_support.f90 @@ -170,6 +170,7 @@ subroutine Get_results(s, k, MLT_option, & ! NOTE: k=0 is a valid arg alpha_semiconvection, thermohaline_coeff, & mixing_type, gradT, Y_face, conv_vel, D, Gamma, ierr) use star_utils + use hydro_rsp2, only: compute_Eq_cell type (star_info), pointer :: s integer, intent(in) :: k character (len=*), intent(in) :: MLT_option @@ -189,13 +190,27 @@ subroutine Get_results(s, k, MLT_option, & ! NOTE: k=0 is a valid arg ! these are used by use_superad_reduction real(dp) :: Gamma_limit, scale_value1, scale_value2, diff_grads_limit, reduction_limit, lambda_limit type(auto_diff_real_star_order1) :: Lrad_div_Ledd, Gamma_inv_threshold, Gamma_factor, alfa0, & - diff_grads_factor, Gamma_term, exp_limit, grad_scale, gradr_scaled + diff_grads_factor, Gamma_term, exp_limit, grad_scale, gradr_scaled, Eq_div_w, check_Eq logical :: test_partials, using_TDC logical, parameter :: report = .false. include 'formats' + ! check if this particular k can be done with TDC + using_TDC = .false. + if (s% MLT_option == 'TDC') using_TDC = .true. + if (.not. s% have_mlt_vc) using_TDC = .false. + if (k <= 0 .or. s%dt <= 0d0) using_TDC = .false. + if (using_TDC) using_TDC = .not. check_if_must_fall_back_to_MLT(s, k) + ! Pre-calculate some things. + Eq_div_w = 0d0 + if (using_TDC) then + if (s% mlt_vc_old(k) > 0d0) then + check_Eq = compute_Eq_cell(s, k, ierr) + Eq_div_w = check_Eq/(s% mlt_vc_old(k)/sqrt_2_div_3) ! maybe should be using conv_vel??? + end if + end if Pr = crad*pow4(T)/3d0 Pg = P - Pr beta = Pg / P @@ -234,14 +249,6 @@ subroutine Get_results(s, k, MLT_option, & ! NOTE: k=0 is a valid arg k, s% solver_iter, s% model_number, gradr%val, grada%val, scale_height%val end if - - ! check if this particular k can be done with TDC - using_TDC = .false. - if (s% MLT_option == 'TDC') using_TDC = .true. - if (.not. s% have_mlt_vc) using_TDC = .false. - if (k <= 0 .or. s%dt <= 0d0) using_TDC = .false. - if (using_TDC) using_TDC = .not. check_if_must_fall_back_to_MLT(s, k) - if (k >= 1) then s% dvc_dt_TDC(k) = 0d0 end if @@ -264,7 +271,7 @@ subroutine Get_results(s, k, MLT_option, & ! NOTE: k=0 is a valid arg call set_TDC(& conv_vel_start, mixing_length_alpha, s% alpha_TDC_DAMP, s%alpha_TDC_DAMPR, s%alpha_TDC_PtdVdt, s%dt, cgrav, m, report, & mixing_type, scale, chiT, chiRho, gradr, r, P, T, rho, dV, Cp, opacity, & - scale_height, gradL, grada, conv_vel, D, Y_face, gradT, s%tdc_num_iters(k), ierr) + scale_height, gradL, grada, conv_vel, D, Y_face, gradT, s%tdc_num_iters(k), Eq_div_w, ierr) s% dvc_dt_TDC(k) = (conv_vel%val - conv_vel_start) / s%dt if (ierr /= 0) then @@ -282,7 +289,7 @@ subroutine Get_results(s, k, MLT_option, & ! NOTE: k=0 is a valid arg call set_TDC(& conv_vel_start, mixing_length_alpha, s% alpha_TDC_DAMP, s%alpha_TDC_DAMPR, s%alpha_TDC_PtdVdt, s%dt, cgrav, m, report, & mixing_type, scale, chiT, chiRho, gradr_scaled, r, P, T, rho, dV, Cp, opacity, & - scale_height, gradL, grada, conv_vel, D, Y_face, gradT, s%tdc_num_iters(k), ierr) + scale_height, gradL, grada, conv_vel, D, Y_face, gradT, s%tdc_num_iters(k), Eq_div_w, ierr) s% dvc_dt_TDC(k) = (conv_vel%val - conv_vel_start) / s%dt if (ierr /= 0) then if (s% report_ierr) write(*,*) 'ierr from set_TDC when using superad_reduction' diff --git a/turb/private/tdc_support.f90 b/turb/private/tdc_support.f90 index 25795ebc1..24aeee5a4 100644 --- a/turb/private/tdc_support.f90 +++ b/turb/private/tdc_support.f90 @@ -65,7 +65,7 @@ module tdc_support logical :: report real(dp) :: mixing_length_alpha, alpha_TDC_DAMP, alpha_TDC_DAMPR, alpha_TDC_PtdVdt, dt type(auto_diff_real_tdc) :: A0, c0, L, L0, gradL, grada - type(auto_diff_real_star_order1) :: T, rho, dV, Cp, kap, Hp, Gamma + type(auto_diff_real_star_order1) :: T, rho, dV, Cp, kap, Hp, Gamma, Eq_div_w end type tdc_info contains @@ -522,7 +522,7 @@ subroutine eval_xis(info, Y, xi0, xi1, xi2) real(dp), parameter :: x_GAMMAR = 2.d0*sqrt(3.d0) S0 = convert(x_ALFAS*info%mixing_length_alpha*info%Cp*info%T/info%Hp)*info%grada - S0 = S0*Y + S0 = S0*Y + convert(info%Eq_div_w) D0 = convert(info%alpha_TDC_DAMP*x_CEDE/(info%mixing_length_alpha*info%Hp)) gammar_div_alfa = info%alpha_TDC_DAMPR*x_GAMMAR/(info%mixing_length_alpha*info%Hp) DR0 = convert(4d0*boltz_sigma*pow2(gammar_div_alfa)*pow3(info%T)/(pow2(info%rho)*info%Cp*info%kap)) diff --git a/turb/public/turb.f90 b/turb/public/turb.f90 index 4c4ab2dff..08c22aab8 100644 --- a/turb/public/turb.f90 +++ b/turb/public/turb.f90 @@ -92,12 +92,12 @@ end subroutine set_thermohaline subroutine set_TDC( & conv_vel_start, mixing_length_alpha, alpha_TDC_DAMP, alpha_TDC_DAMPR, alpha_TDC_PtdVdt, dt, cgrav, m, report, & mixing_type, scale, chiT, chiRho, gradr, r, P, T, rho, dV, Cp, opacity, & - scale_height, gradL, grada, conv_vel, D, Y_face, gradT, tdc_num_iters, ierr) + scale_height, gradL, grada, conv_vel, D, Y_face, gradT, tdc_num_iters, Eq_div_w, ierr) use tdc use tdc_support real(dp), intent(in) :: conv_vel_start, mixing_length_alpha, alpha_TDC_DAMP, alpha_TDC_DAMPR, alpha_TDC_PtdVdt, dt, cgrav, m, scale type(auto_diff_real_star_order1), intent(in) :: & - chiT, chiRho, gradr, r, P, T, rho, dV, Cp, opacity, scale_height, gradL, grada + chiT, chiRho, gradr, r, P, T, rho, dV, Cp, opacity, scale_height, gradL, grada, Eq_div_w logical, intent(in) :: report type(auto_diff_real_star_order1),intent(out) :: conv_vel, Y_face, gradT, D integer, intent(out) :: tdc_num_iters, mixing_type, ierr @@ -139,6 +139,7 @@ subroutine set_TDC( & info%kap = opacity info%Hp = scale_height info%Gamma = Gamma + info%Eq_div_w = Eq_div_w ! Get solution Zub = upper_bound_Z diff --git a/turb/test/src/test_turb.f90 b/turb/test/src/test_turb.f90 index f6a0f5a42..6bb6fd11f 100644 --- a/turb/test/src/test_turb.f90 +++ b/turb/test/src/test_turb.f90 @@ -79,7 +79,7 @@ subroutine compare_TDC_and_Cox_MLT() real(dp) :: mixing_length_alpha, conv_vel_start, alpha_TDC_DAMP, alpha_TDC_DAMPR, alpha_TDC_PtdVdt, dt, cgrav, m, scale type(auto_diff_real_star_order1) :: & r, L, T, P, opacity, rho, dV, chiRho, chiT, Cp, gradr, grada, scale_height, gradL, grav, Lambda - type(auto_diff_real_star_order1) :: gradT, Y_face, conv_vel, D, Gamma + type(auto_diff_real_star_order1) :: gradT, Y_face, conv_vel, D, Gamma, Eq_div_w real(dp) :: Henyey_MLT_nu_param, Henyey_MLT_y_param character(len=3) :: MLT_option integer :: mixing_type, ierr, tdc_num_iters @@ -125,6 +125,7 @@ subroutine compare_TDC_and_Cox_MLT() scale = L%val*1d-3 report = .false. dt = 1d40 ! Long time-step so we get into equilibrium + Eq_div_w = 0d0 ! MLT MLT_option = 'Cox' @@ -136,7 +137,7 @@ subroutine compare_TDC_and_Cox_MLT() call set_TDC(& conv_vel_start, mixing_length_alpha, alpha_TDC_DAMP, alpha_TDC_DAMPR, alpha_TDC_PtdVdt, dt, cgrav, m, report, & mixing_type, scale, chiT, chiRho, gradr, r, P, T, rho, dV, Cp, opacity, & - scale_height, gradL, grada, conv_vel, D, Y_face, gradT, tdc_num_iters, ierr) + scale_height, gradL, grada, conv_vel, D, Y_face, gradT, tdc_num_iters, Eq_div_w, ierr) write(*,1) 'TDC: Y, conv_vel_start, conv_vel, dt ', Y_face%val, conv_vel_start, conv_vel% val, dt @@ -154,7 +155,7 @@ subroutine check_TDC() real(dp) :: mixing_length_alpha, conv_vel_start, alpha_TDC_DAMP, alpha_TDC_DAMPR, alpha_TDC_PtdVdt, dt, cgrav, m, scale type(auto_diff_real_star_order1) :: & r, L, T, P, opacity, rho, dV, chiRho, chiT, Cp, gradr, grada, scale_height, gradL - type(auto_diff_real_star_order1) :: gradT, Y_face, conv_vel, D + type(auto_diff_real_star_order1) :: gradT, Y_face, conv_vel, D, Eq_div_w integer :: mixing_type, ierr, tdc_num_iters logical :: report integer :: j @@ -187,6 +188,7 @@ subroutine check_TDC() chiT = 1d0 chiRho = 1d0 gradr = 3d0 * P * opacity * L / (64 * pi * boltz_sigma * pow4(T) * cgrav * m) + Eq_div_w = 0d0 write(*,*) "####################################" write(*,*) "Running dt test" @@ -196,7 +198,7 @@ subroutine check_TDC() call set_TDC(& conv_vel_start, mixing_length_alpha, alpha_TDC_DAMP, alpha_TDC_DAMPR, alpha_TDC_PtdVdt, dt, cgrav, m, report, & mixing_type, scale, chiT, chiRho, gradr, r, P, T, rho, dV, Cp, opacity, & - scale_height, gradL, grada, conv_vel, D, Y_face, gradT, tdc_num_iters, ierr) + scale_height, gradL, grada, conv_vel, D, Y_face, gradT, tdc_num_iters, Eq_div_w, ierr) write(*,1) 'dt, gradT, conv_vel_start, conv_vel', dt, gradT%val, conv_vel_start, conv_vel% val if (report) stop