diff --git a/star/defaults/controls.defaults b/star/defaults/controls.defaults index d9e317ed1..17b96f9ed 100644 --- a/star/defaults/controls.defaults +++ b/star/defaults/controls.defaults @@ -2284,8 +2284,8 @@ ! mlt_Pturb_factor ! ~~~~~~~~~~~~~~~~ - ! include MLT turbulent pressure at face k = mlt_Pturb_factor*(1/3d0)*0.5*(rho(k) + rho(k-1))*mlt_vc(k)**2/3 - ! MLT turbulent pressure for cell k = avg of values at faces. For TDC, mlt_Pturb_factor = 3d0 + ! For MESA's momentume equation, include MLT turbulent pressure at face k = mlt_Pturb_factor*(1/3d0)*0.5*(rho(k) + rho(k-1))*mlt_vc(k)**2/3 + ! MLT turbulent pressure for cell k = avg of values at faces. For TDC, mlt_Pturb_factor = 2d0 ! corresponds to alpha_TDC_Ptdvdt = 1d0, where Pt = (2/3)*rho(k)*(mlt_vc(k)/sqrt(2/3))**2 ! this replaces conv_dP_term_factor. also see extra_pressure vector and other_pressure routine diff --git a/star/dev_cases_TDC_Pulsation/dev_TDC_Cepheid_6M/inlist_pgstar b/star/dev_cases_TDC_Pulsation/dev_TDC_Cepheid_6M/inlist_pgstar index 305f0b082..a001757aa 100644 --- a/star/dev_cases_TDC_Pulsation/dev_TDC_Cepheid_6M/inlist_pgstar +++ b/star/dev_cases_TDC_Pulsation/dev_TDC_Cepheid_6M/inlist_pgstar @@ -181,9 +181,9 @@ Profile_Panels1_title = '' Profile_Panels1_txt_scale = 0.7 Profile_Panels1_num_panels = 4 -Profile_Panels1_yaxis_name(1) = 'Lc_div_L' +Profile_Panels1_yaxis_name(1) = 'Eq' !Profile_Panels1_dymin(1) = 0.14 -Profile_Panels1_other_yaxis_name(1) = 'Frad_div_cUrad'!'lum_conv' +Profile_Panels1_other_yaxis_name(1) = 'Uq'!'lum_conv' !Profile_Panels1_other_dymin(1) = 0.14 Profile_Panels1_yaxis_name(2) = 'logT'!'logdq'!'radius'!'entropy' diff --git a/star/dev_cases_TDC_Pulsation/dev_TDC_Cepheid_6M/src/run_star_extras.f90 b/star/dev_cases_TDC_Pulsation/dev_TDC_Cepheid_6M/src/run_star_extras.f90 index 3fefedd79..258770b0d 100644 --- a/star/dev_cases_TDC_Pulsation/dev_TDC_Cepheid_6M/src/run_star_extras.f90 +++ b/star/dev_cases_TDC_Pulsation/dev_TDC_Cepheid_6M/src/run_star_extras.f90 @@ -490,6 +490,7 @@ integer function extras_start_step(id) ! if v= 0, turn on v so we can kick if (.not. s%v_flag .and. .not. s%u_flag) then call star_set_v_flag(id, .true., ierr) + write (*,*) 'turning on v_flag hydro for kick' end if call gyre_in_mesa_extras_set_velocities(s, .false., ierr) diff --git a/star/dev_cases_TDC_Pulsation/dev_TDC_Cepheid_6M/src/run_star_extras_TDC_pulsation.inc b/star/dev_cases_TDC_Pulsation/dev_TDC_Cepheid_6M/src/run_star_extras_TDC_pulsation.inc index e3254cb2a..e1d2ca8ca 100644 --- a/star/dev_cases_TDC_Pulsation/dev_TDC_Cepheid_6M/src/run_star_extras_TDC_pulsation.inc +++ b/star/dev_cases_TDC_Pulsation/dev_TDC_Cepheid_6M/src/run_star_extras_TDC_pulsation.inc @@ -158,8 +158,8 @@ v_surf = s% v(1) v_surf_start = s% v_start(1) else if (s% u_flag) then - v_surf = s% u_face_val(1) - v_surf_start = s% u_face_start(1) + v_surf = s% u(1) ! s% u_face_val(1) + v_surf_start = s% u_start(1) ! s% u_face_start(1) else ! v_surf = 0d0 v_surf_start = 0d0 diff --git a/star/private/auto_diff_support.f90 b/star/private/auto_diff_support.f90 index ac5e40539..e926d7728 100644 --- a/star/private/auto_diff_support.f90 +++ b/star/private/auto_diff_support.f90 @@ -1235,6 +1235,39 @@ function wrap_u_p1(s, k) result(v_p1) end if end function wrap_u_p1 + function wrap_opt_time_center_u_m1(s, k) result(v_tc) + type (star_info), pointer :: s + type(auto_diff_real_star_order1) :: v_tc + integer, intent(in) :: k + v_tc = 0d0 + if (k == 1) return + v_tc = wrap_u_m1(s,k) + if (s% using_velocity_time_centering) & + v_tc = 0.5d0*(v_tc + s% u_start(k-1)) + end function wrap_opt_time_center_u_m1 + + function wrap_opt_time_center_u_00(s, k) result(v_tc) + type (star_info), pointer :: s + type(auto_diff_real_star_order1) :: v_tc + integer, intent(in) :: k + v_tc = 0d0 + v_tc = wrap_u_00(s,k) + if (s% using_velocity_time_centering) & + v_tc = 0.5d0*(v_tc + s% u_start(k)) + end function wrap_opt_time_center_u_00 + + function wrap_opt_time_center_u_p1(s, k) result(v_tc) + type (star_info), pointer :: s + type(auto_diff_real_star_order1) :: v_tc + integer, intent(in) :: k + v_tc = 0d0 + if (k == s%nz) return + v_tc = wrap_u_p1(s,k) + if (s% using_velocity_time_centering) & + v_tc = 0.5d0*(v_tc + s% u_start(k+1)) + end function wrap_opt_time_center_u_p1 + + function wrap_Hp_m1(s, k) result(Hp_m1) type (star_info), pointer :: s type(auto_diff_real_star_order1) :: Hp_m1 diff --git a/star/private/hydro_energy.f90 b/star/private/hydro_energy.f90 index ef4cceb36..e24b15d7d 100644 --- a/star/private/hydro_energy.f90 +++ b/star/private/hydro_energy.f90 @@ -216,7 +216,6 @@ end subroutine setup_dL_dm subroutine setup_sources_and_others(ierr) ! sources_ad, others_ad use hydro_rsp2, only: compute_Eq_cell - use star_utils, only: get_face_weights use tdc_hydro, only: compute_tdc_Eq_div_w_face ! compute_Eq_cell real(dp) :: alfa, beta integer, intent(out) :: ierr @@ -267,14 +266,13 @@ subroutine setup_sources_and_others(ierr) ! sources_ad, others_ad Eq_ad = s% Eq_ad(k) ! compute_Eq_cell(s, k, ierr) if (ierr /= 0) return else if (s% alpha_TDC_DampM >0d0 .and. s% MLT_option == 'TDC' .and. & - s% TDC_include_eturb_in_energy_equation) then ! not checking for v or u flag. + s% TDC_include_eturb_in_energy_equation .and. (s% v_flag .or. s% u_flag)) then ! not checking for v or u flag. !Eq_ad = compute_tdc_Eq_cell(s, k, ierr) ! safe to just recompute - if (k == 1) then - Eq_ad = compute_tdc_Eq_div_w_face(s, k, ierr)*(s% mlt_vc_ad(k)/sqrt_2_div_3) - else - call get_face_weights(s, k, alfa, beta) - Eq_ad = alfa*compute_tdc_Eq_div_w_face(s, k, ierr)*(s% mlt_vc_ad(k)/sqrt_2_div_3) + & - beta*shift_m1(compute_tdc_Eq_div_w_face(s, k-1, ierr))*(shift_m1(s% mlt_vc_ad(k-1))/sqrt_2_div_3) + if (k < s% nz) then + Eq_ad = 0.5d0*(compute_tdc_Eq_div_w_face(s, k, ierr)*s% mlt_vc_ad(k) + & + shift_p1(compute_tdc_Eq_div_w_face(s, k+1, ierr))*shift_p1(s% mlt_vc_ad(k+1)))/sqrt_2_div_3 + else ! center cell is 0 at inner face + Eq_ad = 0.5d0*compute_tdc_Eq_div_w_face(s, k, ierr)*s% mlt_vc_ad(k)/sqrt_2_div_3 end if if (ierr /= 0) return end if @@ -355,11 +353,10 @@ subroutine setup_RTI_diffusion(diffusion_eps_ad) end subroutine setup_RTI_diffusion subroutine setup_d_turbulent_energy_dt(ierr) - use star_utils, only: get_face_weights use const_def, only: sqrt_2_div_3 integer, intent(out) :: ierr type(auto_diff_real_star_order1) :: TDC_eturb_cell - real (dp) :: TDC_eturb_cell_start, alfa, beta + real (dp) :: TDC_eturb_cell_start include 'formats' ierr = 0 if (s% RSP2_flag) then @@ -367,15 +364,14 @@ subroutine setup_d_turbulent_energy_dt(ierr) else if (s% mlt_vc_old(k) > 0d0 .and. s% MLT_option == 'TDC' .and. & s% TDC_include_eturb_in_energy_equation) then ! write a wrapper for this. - if (k == 1) then - TDC_eturb_cell_start = pow2(s% mlt_vc_old(k)/sqrt_2_div_3) - TDC_eturb_cell = pow2(s% mlt_vc(k)/sqrt_2_div_3) - else - call get_face_weights(s, k, alfa, beta) - TDC_eturb_cell_start = alfa*pow2(s% mlt_vc_old(k)/sqrt_2_div_3) + & - beta*pow2(s% mlt_vc_old(k-1)/sqrt_2_div_3) - TDC_eturb_cell = alfa*pow2(s% mlt_vc_ad(k)/sqrt_2_div_3) + & - beta*pow2(shift_m1(s% mlt_vc_ad(k-1))/sqrt_2_div_3) + if (k < s% nz) then + TDC_eturb_cell_start = 0.5d0*(pow2(s% mlt_vc_old(k)/sqrt_2_div_3) + & + pow2(s% mlt_vc_old(k+1)/sqrt_2_div_3)) + TDC_eturb_cell = 0.5d0*(pow2(s% mlt_vc_ad(k)/sqrt_2_div_3) + & + pow2(shift_p1(s% mlt_vc_ad(k+1))/sqrt_2_div_3)) + else ! center cell averaged with 0 for inner face + TDC_eturb_cell_start = 0.5d0*pow2(s% mlt_vc_old(k)/sqrt_2_div_3) + TDC_eturb_cell = 0.5d0*pow2(s% mlt_vc(k)/sqrt_2_div_3) end if d_turbulent_energy_dt_ad = (TDC_eturb_cell - TDC_eturb_cell_start) / dt else diff --git a/star/private/hydro_riemann.f90 b/star/private/hydro_riemann.f90 index 5a904ef37..49441a553 100644 --- a/star/private/hydro_riemann.f90 +++ b/star/private/hydro_riemann.f90 @@ -73,18 +73,19 @@ subroutine do1_dudt_eqn( & s, k, P_surf_ad, nvar, ierr) use accurate_sum_auto_diff_star_order1 use star_utils, only: get_area_info_opt_time_center, save_eqn_residual_info + use tdc_hydro, only: compute_tdc_Uq_dm_cell type (star_info), pointer :: s integer, intent(in) :: k type(auto_diff_real_star_order1), intent(in) :: P_surf_ad ! only for k=1 integer, intent(in) :: nvar integer, intent(out) :: ierr - integer :: nz, i_du_dt type(auto_diff_real_star_order1) :: & flux_in_ad, flux_out_ad, diffusion_source_ad, & geometry_source_ad, gravity_source_ad, & area_00, area_p1, inv_R2_00, inv_R2_p1, & - dudt_expected_ad, dudt_actual_ad, resid_ad + dudt_expected_ad, dudt_actual_ad, resid_ad, & + Uq_cell type(accurate_auto_diff_real_star_order1) :: sum_ad real(dp) :: dt, dm, ie_plus_ke, scal, residual logical :: dbg, do_diffusion, test_partials @@ -123,8 +124,15 @@ subroutine do1_dudt_eqn( & call setup_gravity_source call setup_diffusion_source + ! Add turbulent eddy viscous acceleration Uq for TDC as source + Uq_cell = 0d0 + if (s% MLT_option == 'TDC' .and. s%alpha_TDC_DampM > 0d0) then + Uq_cell = compute_tdc_Uq_dm_cell(s, k, ierr) ! Uq * dm + if (ierr /= 0) return + end if + sum_ad = flux_in_ad - flux_out_ad + & - geometry_source_ad + gravity_source_ad + diffusion_source_ad + geometry_source_ad + gravity_source_ad + diffusion_source_ad + Uq_cell dudt_expected_ad = sum_ad dudt_expected_ad = dudt_expected_ad/dm @@ -310,7 +318,6 @@ subroutine do1_uface_and_Pface(s, k, ierr) use eos_def, only: i_gamma1, i_lnfree_e, i_lnPgas use star_utils, only: calc_Ptot_ad_tw, get_face_weights use hydro_rsp2, only: compute_Uq_face - use tdc_hydro, only: compute_tdc_Uq_face type (star_info), pointer :: s integer, intent(in) :: k integer, intent(out) :: ierr @@ -427,16 +434,10 @@ subroutine do1_uface_and_Pface(s, k, ierr) end if - if (s% RSP2_flag) then ! include Uq in u_face + if (s% RSP2_flag) then ! include Uq in u_face, To do: implement in sources instead ~ EbF Uq_ad = compute_Uq_face(s, k, ierr) if (ierr /= 0) return s% u_face_ad(k) = s% u_face_ad(k) + Uq_ad - else if (s% alpha_TDC_DampM >0d0 .and. s% MLT_option == 'TDC') then ! include Uq in u_face - Uq_ad = compute_tdc_Uq_face(s, k, ierr) - if (ierr /= 0) return - s% u_face_ad(k) = s% u_face_ad(k) + Uq_ad - else - Uq_ad = 0d0 end if s% u_face_val(k) = s% u_face_ad(k)%val diff --git a/star/private/hydro_vars.f90 b/star/private/hydro_vars.f90 index bf8f1766d..a3261a742 100644 --- a/star/private/hydro_vars.f90 +++ b/star/private/hydro_vars.f90 @@ -599,7 +599,7 @@ subroutine set_hydro_vars( & s% gradr_factor(nzlo:nzhi) = 1d0 end if - if (s% alpha_TDC_DampM > 0) then + if (s% alpha_TDC_DampM > 0 .and. s% MLT_option == 'TDC') then call set_viscosity_vars_TDC(s,ierr) if (ierr /= 0) then if (len_trim(s% retry_message) == 0) s% retry_message = 'set_viscosity_vars_TDC failed' diff --git a/star/private/read_model.f90 b/star/private/read_model.f90 index 292edad9a..db5bd784f 100644 --- a/star/private/read_model.f90 +++ b/star/private/read_model.f90 @@ -161,7 +161,7 @@ subroutine finish_load_model(s, restart, ierr) s% doing_finish_load_model = .true. call set_vars(s, s% dt, ierr) if (ierr == 0 .and. s% RSP2_flag) call set_RSP2_vars(s,ierr) - if (ierr == 0 .and. s% alpha_TDC_DampM > 0 .and. s% mlt_option == 'TDC') call set_viscosity_vars_TDC(s,ierr) + if (ierr == 0 .and. s% alpha_TDC_DampM > 0 .and. s% MLT_option == 'TDC') call set_viscosity_vars_TDC(s,ierr) s% doing_finish_load_model = .false. if (ierr /= 0) then write(*,*) 'finish_load_model: failed in set_vars' diff --git a/star/private/tdc_hydro.f90 b/star/private/tdc_hydro.f90 index 34afd1b50..302037cee 100644 --- a/star/private/tdc_hydro.f90 +++ b/star/private/tdc_hydro.f90 @@ -32,7 +32,7 @@ module tdc_hydro private public :: & compute_tdc_Eq_cell, compute_tdc_Uq_face, compute_tdc_Eq_div_w_face, & - get_RSP2_alfa_beta_face_weights, set_viscosity_vars_TDC + get_TDC_alfa_beta_face_weights, set_viscosity_vars_TDC, compute_tdc_Uq_dm_cell contains @@ -46,6 +46,15 @@ subroutine set_viscosity_vars_TDC(s, ierr) ierr = 0 op_err = 0 + if (.not. (s%v_flag .or. s%u_flag)) then ! set values 0 if not using v_flag or u_flag. + do k = 1, s%nz + s%Eq(k) = 0d0; s%Eq_ad(k) = 0d0 + s%Chi(k) = 0d0; s%Chi_ad(k) = 0d0 + s%Uq(k) = 0d0 + end do + return + end if + !$OMP PARALLEL DO PRIVATE(k,op_err) SCHEDULE(dynamic,2) do k = 1, s%nz ! Hp_face(k) <= 0 means it needs to be set. e.g., after read file @@ -61,11 +70,19 @@ subroutine set_viscosity_vars_TDC(s, ierr) end if !$OMP PARALLEL DO PRIVATE(k,op_err) SCHEDULE(dynamic,2) do k = 1, s%nz - x = compute_Chi_cell(s, k, op_err) + if (s% v_flag) then + x = compute_Chi_cell(s, k, op_err) + else if (s% u_flag) then + x = compute_Chi_div_w_face(s, k, op_err) + end if if (op_err /= 0) ierr = op_err x = compute_tdc_Eq_cell(s, k, op_err) if (op_err /= 0) ierr = op_err - x = compute_tdc_Uq_face(s, k, op_err) + if (s% v_flag) then + x = compute_tdc_Uq_face(s, k, op_err) + else if (s% u_flag) then + x = compute_tdc_Uq_dm_cell(s, k, op_err) + end if if (op_err /= 0) ierr = op_err end do !$OMP END PARALLEL DO @@ -73,16 +90,9 @@ subroutine set_viscosity_vars_TDC(s, ierr) if (s%report_ierr) write (*, 2) 'failed in set_viscosity_vars_TDC loop 2', s%model_number return end if - if (.not. (s%v_flag .or. s%u_flag)) then ! set values 0 if not using v_flag or u_flag. - do k = 1, s%nz - s%Eq(k) = 0d0; s%Eq_ad(k) = 0d0 - s%Chi(k) = 0d0; s%Chi_ad(k) = 0d0 - s%Uq(k) = 0d0 - end do - end if end subroutine set_viscosity_vars_TDC - subroutine get_RSP2_alfa_beta_face_weights(s, k, alfa, beta) + subroutine get_TDC_alfa_beta_face_weights(s, k, alfa, beta) type(star_info), pointer :: s integer, intent(in) :: k real(dp), intent(out) :: alfa, beta @@ -95,7 +105,7 @@ subroutine get_RSP2_alfa_beta_face_weights(s, k, alfa, beta) alfa = 0.5d0 beta = 0.5d0 end if - end subroutine get_RSP2_alfa_beta_face_weights + end subroutine get_TDC_alfa_beta_face_weights function compute_d_v_div_r(s, k, ierr) result(d_v_div_r) ! s^-1 type(star_info), pointer :: s @@ -127,38 +137,56 @@ function compute_rho_form_of_d_v_div_r(s, k, ierr) result(d_v_div_r) type(star_info), pointer :: s integer, intent(in) :: k integer, intent(out) :: ierr - type(auto_diff_real_star_order1) :: d_v_div_r + type(auto_diff_real_star_order1) :: d_v_div_r, v_00, v_p1 type(auto_diff_real_star_order1) :: r_cell, rho_cell, v_cell, dlnrho_dt real(dp) :: dm_cell ierr = 0 r_cell = 0.5d0*(wrap_r_00(s, k) + wrap_r_p1(s, k)) rho_cell = wrap_d_00(s, k) - v_cell = wrap_v_00(s, k) ! cell-centred velocity (u_flag) + if (s% u_flag) then + v_cell = wrap_u_00(s,k) + else ! v flag + v_cell = 0.5d0*(wrap_v_00(s, k) + wrap_v_p1(s, k)) + end if + v_00 = wrap_opt_time_center_v_00(s, k) + v_p1 = wrap_opt_time_center_v_p1(s, k) dlnrho_dt = wrap_dxh_lnd(s, k)/s%dt ! (∂/∂t)lnρ dm_cell = s%dm(k) ! cell mass d_v_div_r = -dm_cell/(4d0*pi*rho_cell)*(dlnrho_dt/pow3(r_cell) + 3d0*v_cell/pow4(r_cell)) + +! d_v_div_r = ((v_00 - v_p1) - dm_cell*v_cell/(4d0*pi*rho_cell*pow3(r_cell)))/r_cell + end function compute_rho_form_of_d_v_div_r function compute_rho_form_of_d_v_div_r_opt_time_center(s, k, ierr) result(d_v_div_r) ! s^-1 type(star_info), pointer :: s integer, intent(in) :: k integer, intent(out) :: ierr - type(auto_diff_real_star_order1) :: d_v_div_r + type(auto_diff_real_star_order1) :: d_v_div_r, v_00, v_p1 type(auto_diff_real_star_order1) :: r_cell, rho_cell, v_cell, dlnrho_dt real(dp) :: dm_cell ierr = 0 r_cell = 0.5d0*(wrap_opt_time_center_r_00(s, k) + wrap_opt_time_center_r_p1(s, k)) rho_cell = wrap_d_00(s, k) - v_cell = wrap_opt_time_center_v_00(s, k) ! cell-centred velocity (u_flag) + if (s% u_flag) then + v_cell = wrap_opt_time_center_u_00(s,k) + else ! v flag + v_cell = 0.5d0*(wrap_opt_time_center_v_00(s, k) + wrap_opt_time_center_v_p1(s, k)) + end if + v_00 = wrap_opt_time_center_v_00(s, k) + v_p1 = wrap_opt_time_center_v_p1(s, k) + dlnrho_dt = wrap_dxh_lnd(s, k)/s%dt ! (∂/∂t)lnρ dm_cell = s%dm(k) ! cell mass - d_v_div_r = -dm_cell/(4d0*pi*rho_cell)* & - (dlnrho_dt/pow3(r_cell) & - + 3d0*v_cell/pow4(r_cell)) + d_v_div_r = -dm_cell/(4d0*pi*rho_cell)*(dlnrho_dt/pow3(r_cell) + 3d0*v_cell/pow4(r_cell)) + + ! want this on the cell now instead of face. +! d_v_div_r = ((v_00 - v_p1) - dm_cell*v_cell/(4d0*pi*rho_cell*pow3(r_cell)))/r_cell + end function compute_rho_form_of_d_v_div_r_opt_time_center @@ -168,25 +196,26 @@ function compute_rho_form_of_d_v_div_r_face(s, k, ierr) result(d_v_div_r) integer, intent(out) :: ierr type(auto_diff_real_star_order1) :: d_v_div_r type(auto_diff_real_star_order1) :: r_face, rho_face, v_face, dlnrho_dt - real(dp) :: dm_bar + real(dp) :: dm_bar, alfa, beta ierr = 0 r_face = wrap_r_00(s, k) rho_face = get_rho_face(s, k) - if (s% v_flag .and. .not. s% u_flag) then - v_face = wrap_v_00(s, k) ! face-centred velocity (v_flag) - else if (s% u_flag .and. .not. s% v_flag) then - v_face = s% u_face_ad(k) ! reconstructed face velocity (u_flag) - end if - dlnrho_dt = 0.5d0*(wrap_dxh_lnd(s, k) + wrap_dxh_lnd(s, k+1))/s%dt ! (∂/∂t)lnρ - + v_face = wrap_v_00(s, k) ! face-centered velocity if (k >= 2) then dm_bar = 0.5d0*(s% dm(k) + s% dm(k-1)) + call get_TDC_alfa_beta_face_weights(s, k, alfa, beta) + dlnrho_dt = (alfa*wrap_dxh_lnd(s, k) + beta*shift_m1(wrap_dxh_lnd(s, k-1)))/s%dt ! (∂/∂t)lnρ else dm_bar = 0.5d0*s% dm(k) + dlnrho_dt = 0.5d0*wrap_dxh_lnd(s, k)/s%dt ! (∂/∂t)lnρ end if - d_v_div_r = -dm_bar/(4d0*pi*rho_face)*(dlnrho_dt/pow3(r_face) + 3d0*v_face/pow4(r_face)) +! d_v_div_r = -dm_bar/(4d0*pi*rho_face)*(dlnrho_dt/pow3(r_face) + 3d0*v_face/pow4(r_face)) + + ! 1/r * du/dm - U/4/pi/rho/r^4 + d_v_div_r = ((wrap_u_m1(s,k) - wrap_u_00(s,k)) - dm_bar*v_face/(4d0*pi*rho_face*pow3(r_face)))/r_face + end function compute_rho_form_of_d_v_div_r_face function compute_rho_form_of_d_v_div_r_face_opt_time_center(s, k, ierr) result(d_v_div_r) ! s^-1 @@ -195,26 +224,26 @@ function compute_rho_form_of_d_v_div_r_face_opt_time_center(s, k, ierr) result(d integer, intent(out) :: ierr type(auto_diff_real_star_order1) :: d_v_div_r type(auto_diff_real_star_order1) :: r_face, rho_face, v_face, dlnrho_dt - real(dp) :: dm_bar + real(dp) :: dm_bar, alfa, beta ierr = 0 r_face = wrap_opt_time_center_r_00(s, k) rho_face = get_rho_face(s, k) - - if (s% v_flag .and. .not. s% u_flag) then - v_face = wrap_opt_time_center_v_00(s, k) ! face-centred velocity (v_flag) - else if (s% u_flag .and. .not. s% v_flag) then - v_face = s% u_face_ad(k) ! reconstructed face velocity (u_flag) - end if - dlnrho_dt = 0.5d0*(wrap_dxh_lnd(s, k) + wrap_dxh_lnd(s, k+1))/s%dt ! (∂/∂t)lnρ - + v_face = wrap_opt_time_center_v_00(s, k) ! face-centered velocity if (k >= 2) then dm_bar = 0.5d0*(s% dm(k) + s% dm(k-1)) + call get_TDC_alfa_beta_face_weights(s, k, alfa, beta) + dlnrho_dt = (alfa*wrap_dxh_lnd(s, k) + beta*shift_m1(wrap_dxh_lnd(s, k-1)))/s%dt ! (∂/∂t)lnρ else dm_bar = 0.5d0*s% dm(k) + dlnrho_dt = 0.5d0*wrap_dxh_lnd(s, k)/s%dt ! (∂/∂t)lnρ end if - d_v_div_r = -dm_bar/(4d0*pi*rho_face)*(dlnrho_dt/pow3(r_face) + 3d0*v_face/pow4(r_face)) +! d_v_div_r = -dm_bar/(4d0*pi*rho_face)*(dlnrho_dt/pow3(r_face) + 3d0*v_face/pow4(r_face)) + + ! dm_bar*(1/r * du/dm - U/4/pi/rho/r^4) + d_v_div_r = ((wrap_opt_time_center_u_m1(s,k) - wrap_opt_time_center_u_00(s,k)) - dm_bar*v_face/(4d0*pi*rho_face*pow3(r_face)))/r_face + end function compute_rho_form_of_d_v_div_r_face_opt_time_center @@ -241,7 +270,7 @@ function wrap_Hp_cell(s, k) result(Hp_cell) ! cm , different than rsp2 Hp0 = get_scale_height_face(s,k) Hp1 = 0d0 if (k+1 < s%nz) then - Hp1 = get_scale_height_face(s,k+1) + Hp1 = shift_p1(get_scale_height_face(s,k+1)) end if Hp_cell = 0.5d0*(Hp0 + Hp1) !0.5d0*(wrap_Hp_00(s, k) + wrap_Hp_p1(s, k)) @@ -317,10 +346,18 @@ function compute_Chi_cell(s, k, ierr) result(Chi_cell) if (ierr /= 0) return ! don't need to check if mlt_vc > 0 here. - if (s%MLT_option == 'TDC' .and. .not. s%RSP2_flag) then - w_00 = s% mlt_vc_ad(k)/sqrt_2_div_3! same as info%A0 from TDC - else ! normal RSP2 - w_00 = wrap_w_00(s, k) + if (k < s% nz) then + if (s% okay_to_set_mlt_vc .and. s% have_mlt_vc) then !add option for explicit mlt_vc, operator split in momentum eq. + w_00 = 0.5d0*(s% mlt_vc_old(k) + s% mlt_vc_old(k+1))/sqrt_2_div_3! same as info%A0 from TDC + else + w_00 = 0.5d0*(s% mlt_vc_ad(k) + shift_p1(s% mlt_vc_ad(k+1)))/sqrt_2_div_3! same as info%A0 from TDC + end if + else + if (s% okay_to_set_mlt_vc .and. s% have_mlt_vc) then !add option for explicit mlt_vc, operator split in momentum eq. + w_00 = 0.5d0*s% mlt_vc_old(k)/sqrt_2_div_3! same as info%A0 from TDC + else + w_00 = 0.5d0*s% mlt_vc_ad(k)/sqrt_2_div_3! same as info%A0 from TDC + end if end if d_00 = wrap_d_00(s, k) f = (16d0/3d0)*pi*ALFAM_ALFA/s%dm(k) @@ -444,10 +481,6 @@ function compute_Chi_div_w_face(s, k, ierr) result(Chi_face) if (ALFAM_ALFA == 0d0 .or. & k > s%nz - s% TDC_num_innermost_cells_forced_nonturbulent) then Chi_face = 0d0 - if (k >= 1 .and. k <= s%nz) then - s%Chi(k) = 0d0 - s%Chi_ad(k) = 0d0 - end if else Hp_face = get_scale_height_face(s,k) !Hp_cell_for_Chi(s, k, ierr) if (ierr /= 0) return @@ -476,8 +509,12 @@ function compute_Chi_div_w_face(s, k, ierr) result(Chi_face) ! = erg ! * (s / cm) - > [erg] * 1/w00 end if - !s%Chi(k) = Chi_face%val - !s%Chi_ad(k) = Chi_face + + ! since chi_cell is not called for u_flag when doing Uq_cell_dm + if (s% u_flag .and. .not. s% v_flag) then ! when v_flag, set inside chi_cell + s%Chi(k) = Chi_face%val + s%Chi_ad(k) = Chi_face + end if if (dbg .and. k == -100) then write (*, *) ' s% ALFAM_ALFA', ALFAM_ALFA @@ -500,7 +537,7 @@ function compute_tdc_Eq_div_w_face(s, k, ierr) result(Eq_face) ! erg g^-1 s^-1 integer, intent(in) :: k type(auto_diff_real_star_order1) :: Eq_face integer, intent(out) :: ierr - type(auto_diff_real_star_order1) :: d_v_div_r, Chi_face + type(auto_diff_real_star_order1) :: d_v_div_r, Chi_face, w_00 real(dp) :: dmbar include 'formats' ierr = 0 @@ -528,8 +565,16 @@ function compute_tdc_Eq_div_w_face(s, k, ierr) result(Eq_face) ! erg g^-1 s^-1 if (ierr /= 0) return Eq_face = 4d0*pi*Chi_face*d_v_div_r/dmbar ! erg s^-1 g^-1 * (cm^-1 s^1) end if - !s%Eq(k) = Eq_face%val - !s%Eq_ad(k) = Eq_face + + ! only for output, really only used for returning Eq to star pointers. + if (s% okay_to_set_mlt_vc .and. s% have_mlt_vc) then !add option for explicit mlt_vc, operator split in momentum eq. + w_00 = s% mlt_vc_old(k)/sqrt_2_div_3! same as info%A0 from TDC + else + w_00 = s% mlt_vc_ad(k)/sqrt_2_div_3! same as info%A0 from TDC + end if + + s%Eq(k) = Eq_face%val * w_00%val + s%Eq_ad(k) = Eq_face * w_00 end function compute_tdc_Eq_div_w_face function compute_d_v_div_r_face(s, k, ierr) result(d_v_div_r) ! s^-1 @@ -537,68 +582,125 @@ function compute_d_v_div_r_face(s, k, ierr) result(d_v_div_r) ! s^-1 integer, intent(in) :: k type(auto_diff_real_star_order1) :: d_v_div_r integer, intent(out) :: ierr - type(auto_diff_real_star_order1) :: v_00, v_p1, r_00, r_p1, term1, term2 + type(auto_diff_real_star_order1) :: v_00, v_m1, r_00, r_m1, term1, term2 logical :: dbg include 'formats' ierr = 0 dbg = .false. - if (k >1) then - v_00 = 0.5d0 * (wrap_v_00(s, k) + wrap_v_m1(s, k)) - else - v_00 = 0.5d0 * wrap_v_00(s, k) + if (s% v_flag) then + v_00 = 0.5d0 * (wrap_v_00(s, k) + wrap_v_p1(s, k)) + v_m1 = 0.5d0*(wrap_v_00(s, k) + wrap_v_m1(s, k)) + else if(s% u_flag) then + v_00 = wrap_u_00(s,k) + v_m1 = wrap_u_m1(s,k) end if - v_p1 = 0.5d0 * (wrap_v_00(s, k) + wrap_v_p1(s, k)) + r_00 = 0.5d0*(wrap_r_00(s, k) + wrap_r_p1(s, k)) + r_m1 = 0.5d0*(wrap_r_00(s, k) + wrap_r_m1(s, k)) - if (k >1) then - r_00 = 0.5d0 * (wrap_r_00(s, k) + wrap_r_m1(s, k)) - else - r_00 = 0.5d0 * wrap_r_00(s, k) - end if - r_p1 = 0.5d0*(wrap_r_00(s, k) + wrap_r_p1(s, k)) - if (r_p1%val == 0d0) r_p1 = 1d0 - d_v_div_r = v_00/r_00 - v_p1/r_p1 ! units s^-1 + if (r_00%val == 0d0) r_00 = 1d0 + if (r_m1%val == 0d0) r_m1 = 1d0 + d_v_div_r = v_m1/r_m1 - v_00/r_00 ! units s^-1 ! Debugging output to trace values if (dbg .and. k == -63) then write (*, *) 'test d_v_div_r, k:', k - write (*, *) 'v_00:', v_00%val, 'v_p1:', v_p1%val - write (*, *) 'r_00:', r_00%val, 'r_p1:', r_p1%val + write (*, *) 'v_00:', v_00%val, 'v_p1:', v_m1%val + write (*, *) 'r_00:', r_00%val, 'r_p1:', r_m1%val write (*, *) 'd_v_div_r:', d_v_div_r%val end if end function compute_d_v_div_r_face - function compute_d_v_div_r_opt_time_center_face(s, k, ierr) result(d_v_div_r) ! s^-1 type(star_info), pointer :: s integer, intent(in) :: k type(auto_diff_real_star_order1) :: d_v_div_r integer, intent(out) :: ierr - type(auto_diff_real_star_order1) :: v_00, v_p1, r_00, r_p1 + type(auto_diff_real_star_order1) :: v_00, v_m1, r_00, r_m1, term1, term2 + logical :: dbg include 'formats' ierr = 0 + dbg = .false. - if (k >1) then - v_00 = 0.5d0 * (wrap_opt_time_center_v_00(s, k) + wrap_opt_time_center_v_m1(s, k)) - else - v_00 = 0.5d0 * wrap_opt_time_center_v_00(s, k) + if (s% v_flag) then + v_00 = 0.5d0 *(wrap_opt_time_center_v_00(s, k) + wrap_opt_time_center_v_p1(s, k)) + v_m1 = 0.5d0*(wrap_opt_time_center_v_00(s, k) + wrap_opt_time_center_v_m1(s, k)) + else if(s% u_flag) then + v_00 = wrap_opt_time_center_u_00(s,k) + v_m1 = wrap_opt_time_center_u_m1(s,k) end if + + r_00 = 0.5d0*(wrap_opt_time_center_r_00(s, k) + wrap_opt_time_center_r_p1(s, k)) + r_m1 = 0.5d0*(wrap_opt_time_center_r_00(s, k) + wrap_opt_time_center_r_m1(s, k)) - v_p1 = 0.5d0 * (wrap_opt_time_center_v_00(s, k) + wrap_opt_time_center_v_p1(s, k)) + if (r_00%val == 0d0) r_00 = 1d0 + if (r_m1%val == 0d0) r_m1 = 1d0 + d_v_div_r = v_m1/r_m1 - v_00/r_00 ! units s^-1 - if (k >1) then - r_00 = 0.5d0 * (wrap_opt_time_center_r_00(s, k) + wrap_opt_time_center_r_m1(s, k)) - else - r_00 = 0.5d0 * wrap_opt_time_center_r_00(s, k) + ! Debugging output to trace values + if (dbg .and. k == -63) then + write (*, *) 'test d_v_div_r, k:', k + write (*, *) 'v_00:', v_00%val, 'v_p1:', v_m1%val + write (*, *) 'r_00:', r_00%val, 'r_p1:', r_m1%val + write (*, *) 'd_v_div_r:', d_v_div_r%val end if - r_p1 = 0.5d0*(wrap_opt_time_center_r_00(s, k) + wrap_opt_time_center_r_p1(s, k)) - if (r_p1%val == 0d0) r_p1 = 1d0 - d_v_div_r = v_00/r_00 - v_p1/r_p1 ! units s^-1 - - if (r_p1%val == 0d0) r_p1 = 1d0 - d_v_div_r = v_00/r_00 - v_p1/r_p1 ! units s^-1 end function compute_d_v_div_r_opt_time_center_face + + + ! for u_flag only. cell centered Uq as source for Reimann flux. + function compute_tdc_Uq_dm_cell(s, k, ierr) result(Uq_cell) ! cm s^-2, acceleration + type(star_info), pointer :: s + integer, intent(in) :: k + integer, intent(out) :: ierr + type(auto_diff_real_star_order1) :: Chi_00, Chi_p1, r_00, r_p1, w_00, w_p1, r_cell, Uq_cell + include 'formats' + ierr = 0 + if (s%mixing_length_alpha == 0d0 .or. & + k <= s%RSP2_num_outermost_cells_forced_nonturbulent .or. & + k > s%nz - int(s%nz/s%RSP2_nz_div_IBOTOM)) then + Uq_cell = 0d0 + else + r_00 = wrap_opt_time_center_r_00(s, k) + r_p1 = wrap_opt_time_center_r_p1(s, k) + r_cell = 0.5d0*(r_00+r_p1) + + if (s% okay_to_set_mlt_vc .and. s% have_mlt_vc) then !add option for explicit mlt_vc, operator split in momentum eq. + w_00 = s% mlt_vc_old(k)/sqrt_2_div_3! same as info%A0 from TDC + else + w_00 = s% mlt_vc_ad(k)/sqrt_2_div_3! same as info%A0 from TDC + end if + + Chi_00 = compute_Chi_div_w_face(s, k, ierr) * w_00 ! s% Chi_ad(k) XXX + !Chi_00 = s% Chi_ad(k) ! compute_Chi_cell(s,k,ierr) + + if (k < s% nz) then + if (s% okay_to_set_mlt_vc .and. s% have_mlt_vc) then !add option for explicit mlt_vc, operator split in momentum eq. + w_p1 = s% mlt_vc_old(k+1)/sqrt_2_div_3! same as info%A0 from TDC + else + w_p1 = shift_p1(s% mlt_vc_ad(k+1))/sqrt_2_div_3! same as info%A0 from TDC + end if + + Chi_p1 = shift_p1(compute_Chi_div_w_face(s, k + 1, ierr))*w_p1 + !Chi_m1 = shift_m1(s% Chi_ad(k-1)) XXX + if (ierr /= 0) return + else + Chi_p1 = 0d0 + w_p1 = 0d0 + end if + + Uq_cell = 4d0*pi*(Chi_00 - Chi_p1)/(r_cell) ! we have neglected the /dm here, because it is restored in the reimann flux calculation + + if (k == -56) then + write (*, 3) 'RSP2 Uq chi_m1 chi_00 r', k, s%solver_iter, & + Uq_cell%val, Chi_p1%val, Chi_00%val, r_00%val + end if + + end if + ! erg g^-1 cm^-1 = g cm^2 s^-2 g^-1 cm^-1 = cm s^-2, acceleration + s%Uq(k) = Uq_cell%val/ s% dm(k) + end function compute_tdc_Uq_dm_cell + end module tdc_hydro diff --git a/star/private/tdc_hydro_support.f90 b/star/private/tdc_hydro_support.f90 index 731d19c7b..161d1f679 100644 --- a/star/private/tdc_hydro_support.f90 +++ b/star/private/tdc_hydro_support.f90 @@ -475,7 +475,7 @@ subroutine revise_lnT_for_QHSE2(P_surf, ierr) end subroutine revise_lnT_for_QHSE2 subroutine set_Hp_face2(k) - use tdc_hydro, only: get_RSP2_alfa_beta_face_weights + use tdc_hydro, only: get_TDC_alfa_beta_face_weights integer, intent(in) :: k real(dp) :: r_00, d_00, Peos_00, Peos_div_rho, Hp_face, & d_m1, Peos_m1, alfa, beta @@ -488,7 +488,7 @@ subroutine set_Hp_face2(k) else d_m1 = s%rho(k - 1) Peos_m1 = s%Peos(k - 1) - call get_RSP2_alfa_beta_face_weights(s, k, alfa, beta) + call get_TDC_alfa_beta_face_weights(s, k, alfa, beta) Peos_div_rho = alfa*Peos_00/d_00 + beta*Peos_m1/d_m1 Hp_face = pow2(r_00)*Peos_div_rho/(s%cgrav(k)*s%m(k)) end if