diff --git a/docs/source/changelog.rst b/docs/source/changelog.rst index 8ef03e1a5..676cdb988 100644 --- a/docs/source/changelog.rst +++ b/docs/source/changelog.rst @@ -26,6 +26,13 @@ Removed ``file_extension`` option because it is redundant with ``file_device``. Renamed pgstar ``pause``` option to ``pause_flag``` because pause is a reserved Fortran 77 keyword. +For greater consistency and clarity between ``TDC`` and ``RSP``, the controls for ``mlt_option = TDC`` have been renamed to match their respective ``RSP`` +counter part. See below : + +- ``alpha_TDC_DAMP`` has been renamed to ``TDC_alpha_D`` and is analogous to ``RSP_alfad`` +- ``alpha_TDC_DAMPR`` has been renamed to ``TDC_alpha_R`` and is analogous to ``RSP_gammar`` +- ``alpha_TDC_PtdVdt`` has been renamed to ``TDC_alpha_Pt`` and is analogous to ``RSP_alfap`` + .. _New Features main: New Features diff --git a/star/defaults/controls.defaults b/star/defaults/controls.defaults index d9e317ed1..3f3617464 100644 --- a/star/defaults/controls.defaults +++ b/star/defaults/controls.defaults @@ -2284,9 +2284,9 @@ ! 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 - ! corresponds to alpha_TDC_Ptdvdt = 1d0, where Pt = (2/3)*rho(k)*(mlt_vc(k)/sqrt(2/3))**2 + ! 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 TDC_alpha_Pt = 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 @@ -2328,16 +2328,16 @@ ! TDC options ! ~~~~~~~~~~~ - ! + ``alpha_TDC_DAMP`` : The turbulent viscous damping parameter which determines the saturation of TDC. Increasing this decreases convection speeds. - ! + ``alpha_TDC_DAMPR`` : The radiative damping parameter which determines the saturation of TDC. Increasing this decreases convection speeds. - ! + ``alpha_TDC_PtdVdt`` : The prefactor on the term accounting for work done against turbulent pressure (P_turb * dV/dt). Physically this should be unity. + ! + ``TDC_alpha_D`` : The turbulent viscous damping parameter which determines the saturation of TDC. Increasing this decreases convection speeds. This control is analogous to ``RSP_alfad``, and previously called ``alpha_TDC_DAMP`` + ! + ``TDC_alpha_R`` : The radiative damping parameter which determines the saturation of TDC. Increasing this decreases convection speeds. This control is analogous to ``RSP_gammar``, and previously called ``alpha_TDC_DAMPR`` + ! + ``TDC_alpha_Pt`` : The prefactor on the term accounting for work done against turbulent pressure (P_turb * dV/dt). Physically this should be unity. This control is analogous to ``RSP_alfap``, and previously called ``alpha_TDC_PtdVdt`` ! + ``steps_before_use_TDC`` : TDC often struggles with models on the pre-main-sequence. Set this option to pick MLT_option='Cox' for the first several steps to get past the pre-MS. Note that if this option is positive then only either TDC or Cox will be used (depending on model number). THIS OVERRIDES MLT_option! ! :: - alpha_TDC_DAMP = 1d0 - alpha_TDC_DAMPR = 0d0 - alpha_TDC_PtdVdt = 0d0 + TDC_alpha_D = 1d0 + TDC_alpha_R = 0d0 + TDC_alpha_Pt = 0d0 steps_before_use_TDC = 0 diff --git a/star/defaults/controls_dev.defaults b/star/defaults/controls_dev.defaults index a7ab75c1e..3a8719025 100644 --- a/star/defaults/controls_dev.defaults +++ b/star/defaults/controls_dev.defaults @@ -180,40 +180,48 @@ ! === - ! alpha_TDC_DAMPM - ! ~~~~~~~~~~~~~~~ - ! alpha_TDC_C + ! TDC_alpha_M + ! ~~~~~~~~~~~ + ! TDC_alpha_C ! ~~~~~~~~~~~ - ! alpha_TDC_S + ! TDC_alpha_S ! ~~~~~~~~~~~ ! TDC_use_density_form_for_eddy_viscosity ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + ! TDC_alpha_M_use_explicit_mlt_vc_in_momentum_equation + ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ - - ! If ``alpha_TDC_DAMPM>0``, then include eddy viscous damping in TDC ``alpha_TDC_DAMPM`` + ! If ``TDC_alpha_M>0``, then include eddy viscous damping in TDC ``TDC_alpha_M`` ! This control is analogous to ``RSP_alfam``, where the default is ``RSP_alfam = 0.25d0``. ! If hydrostatic (``v_flag``, ``u_flag = .false.``, ``v = 0`` ) there are no velocity gradients, ! and thus no shear to drive turbulence. Without shear, the eddy viscosity term becomes zero. - ! ``alpha_TDC_C`` and ``alpha_TDC_S`` are pre-factors to scale the flux and source terms. + ! ``TDC_alpha_C`` and ``TDC_alpha_S`` are pre-factors to scale the flux and source terms. ! If ``u_flag = .true.`` or ``TDC_use_density_form_for_eddy_viscosity = .true.``, use density ! derivative from newton solver to form d(v/r)/dr, used to compute Eq and Uq. + ! if ``TDC_alpha_M_use_explicit_mlt_vc_in_momentum_equation = .true.``, ``mlt_vc_old`` (the start of step value) + ! is used in place of the current solver iterate ``mlt_vc_ad`` inside the hydro momentum equation. This option can + ! provide greater numerical stability in extremely dynamic situations, with very little loss in accuracy. + ! :: - alpha_TDC_DAMPM = 0d0 - alpha_TDC_C = 1.0d0 - alpha_TDC_S = 1.0d0 + TDC_alpha_M = 0d0 + TDC_alpha_C = 1.0d0 + TDC_alpha_S = 1.0d0 TDC_use_density_form_for_eddy_viscosity = .false. - + TDC_alpha_M_use_explicit_mlt_vc_in_momentum_equation = .false. ! TDC_include_eturb_in_energy_equation ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ - ! calculates and incorporates det/dt from TDC as a negative source term in - ! the total energy equation, along with the eddy viscous heating Eq + ! While ``TDC`` is already implicitly coupled to the total energy equation, this + ! controls feedsback the turbulent energy into the total energy equation by + ! calculating and incorporating det/dt from TDC as a negative source term in + ! the total energy equation, along with the eddy viscous heating Eq. This comes + ! at the cost of numerical stability on long timescales, be warned. ! :: @@ -247,14 +255,18 @@ ! TDC_num_innermost_cells_forced_nonturbulent ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + ! TDC_num_outermost_cells_forced_nonturbulent + ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ - ! Forces innermost TDC_num_innermost_cells_forced_nonturbulent cells + ! Forces innermost TDC_num_innermost_cells_forced_nonturbulent or + ! TDC_num_outermost_cells_forced_nonturbulent cells ! to be nonturbulent, and sets gradT = gradr in these cells. ! Useful for pulsation models in TDC. ! :: TDC_num_innermost_cells_forced_nonturbulent = 0 + TDC_num_outermost_cells_forced_nonturbulent = 0 ! include_mlt_corr_to_TDC @@ -268,6 +280,16 @@ include_mlt_corr_to_TDC = .true. + ! use_TDC_enthalpy_flux_limiter + ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + + ! If use_TDC_enthalpy_flux_limiter = .true. , apply enthalpy flux limiter to TDC + ! similar to the form presented in Wuchterl & Feuchtinger 1998, and Smolec 2008. + + ! :: + + use_TDC_enthalpy_flux_limiter = .false. + ! include_mlt_Pturb_in_thermodynamic_gradients ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ @@ -293,6 +315,9 @@ include_mlt_in_velocity_time_centering = .false. +! Remeshing +! ========= + ! remesh_for_TDC_pulsations_log_core_zoning ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ @@ -302,9 +327,6 @@ ! :: - ! Remeshing parameters - ! ~~~~~~~~~~~~~~~~~~~~ - remesh_for_TDC_pulsations_log_core_zoning = .false. TDC_hydro_use_mass_interp_face_values = .false. @@ -312,3 +334,15 @@ TDC_hydro_nz_outer = 40 TDC_hydro_T_anchor = 11d3 TDC_hydro_dq_1_factor = 2d0 + + + ! use_hydro_merge_limits_in_mesh_plan + ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + + ! Allows use of ``merge_amr_max_abs_du_div_cs``, ``merge_amr_du_div_cs_limit_only_for_compression``, + ! and ``merge_amr_inhibit_at_jumps`` in default mesh_plan when ``use_split_merge_amr = .false.``. + ! Only supports ``v_flag``. + + ! :: + + use_hydro_merge_limits_in_mesh_plan = .false. diff --git a/star/defaults/profile_columns.list b/star/defaults/profile_columns.list index eb34d84fa..89e4d6d5c 100644 --- a/star/defaults/profile_columns.list +++ b/star/defaults/profile_columns.list @@ -420,6 +420,7 @@ !cell_internal_energy_fraction !cell_internal_energy_fraction_start !cell_specific_PE + !dwork_dm ! cell specific work per unit time, Work = dwork_dm*dm*dt !log_cell_ie_div_star_ie !log_cell_specific_IE diff --git a/star/dev_cases_TDC_Pulsation/dev_TDC_Cepheid_6M/he_dep.mod b/star/dev_cases_TDC_Pulsation/dev_TDC_Cepheid_6M/he_dep.mod index 5ca3ecb05..42abcf322 100644 --- a/star/dev_cases_TDC_Pulsation/dev_TDC_Cepheid_6M/he_dep.mod +++ b/star/dev_cases_TDC_Pulsation/dev_TDC_Cepheid_6M/he_dep.mod @@ -1,3 +1,3 @@ version https://git-lfs.github.com/spec/v1 -oid sha256:c978df6f7409a68f159431f320e6be00ee8a6c6fe0bb7ccf86ca878532c6576a -size 771380 +oid sha256:a84daed587eca4e68a8dc2bb2217eaf2cb9eb19c1259f299edcf39542e6bd4d0 +size 777148 diff --git a/star/dev_cases_TDC_Pulsation/dev_TDC_Cepheid_6M/inlist_common b/star/dev_cases_TDC_Pulsation/dev_TDC_Cepheid_6M/inlist_common index a6b0e7690..1485f1f47 100644 --- a/star/dev_cases_TDC_Pulsation/dev_TDC_Cepheid_6M/inlist_common +++ b/star/dev_cases_TDC_Pulsation/dev_TDC_Cepheid_6M/inlist_common @@ -124,11 +124,14 @@ use_Type2_opacities = .true. ! convection controls MLT_option = 'TDC' - alpha_TDC_DAMP = 1d0 - alpha_TDC_DAMPR = 0d0 - alpha_TDC_PtdVdt = 0d0 - alpha_TDC_DAMPM = 0.25d0 + TDC_alpha_D = 1d0 + TDC_alpha_R = 0d0 + TDC_alpha_Pt = 0d0 + TDC_alpha_M = 0.25d0 + TDC_alpha_M_use_explicit_mlt_vc_in_momentum_equation = .false. steps_before_use_TDC = 0 + include_mlt_corr_to_TDC = .false. + TDC_include_eturb_in_energy_equation = .false. ! unstable for evol use_ledoux_criterion = .true. alpha_semiconvection = 1d0 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..4a101208a 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 @@ -3,7 +3,7 @@ !pause = .true. -pgstar_interval = 50 ! making this too small slows the model down. +pgstar_interval = 200 ! making this too small slows the model down. pgstar_show_age_in_years = .true. pgstar_show_age_in_seconds = .false. pgstar_sleep = 0.0 @@ -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/inlist_pulses b/star/dev_cases_TDC_Pulsation/dev_TDC_Cepheid_6M/inlist_pulses index 2353fdec1..02b1836f8 100644 --- a/star/dev_cases_TDC_Pulsation/dev_TDC_Cepheid_6M/inlist_pulses +++ b/star/dev_cases_TDC_Pulsation/dev_TDC_Cepheid_6M/inlist_pulses @@ -52,7 +52,7 @@ cubic_interpolation_in_Z = .false. &controls ! limit max model number as part of test suite. - max_model_number = 3000 + !max_model_number = 3000 ! controls for analyzing pulsations using GYRE in MESA x_integer_ctrl(1) = 1000 ! gyre interval to check @@ -66,7 +66,7 @@ cubic_interpolation_in_Z = .false. ! new rsp style meshing, thanks to Bill P. - x_logical_ctrl(23) = .true. ! .true. = Remesh for TDC + x_logical_ctrl(23) = .true. ! .true. = Remesh for TDC TDC_hydro_use_mass_interp_face_values = .true. TDC_hydro_nz = 150 TDC_hydro_nz_outer = 40 @@ -75,9 +75,13 @@ cubic_interpolation_in_Z = .false. remesh_for_TDC_pulsations_log_core_zoning = .false. ! .false means do rsp style core okay_to_remesh = .false. ! freeze mesh after initial remesh. - x_logical_ctrl(24) = .true. ! if true turn off remesh at the following model number + x_logical_ctrl(24) = .false. ! if true turn off remesh at the following model number x_ctrl(12) = 200 ! model number to turn off remesh ( only if if okay_to_remesh = .true.) + !max_num_merge_surface_cells = 5 + !max_surface_cell_dq= 1d-4 + !min_surface_cell_dq = 1d-6 + ! GYRE set starting velocities, kick! ! kick when true and not restarting. x_logical_ctrl(5) = .true. ! to turn on gyre kick @@ -151,16 +155,16 @@ cubic_interpolation_in_Z = .false. include_mlt_in_velocity_time_centering = .false. include_mlt_Pturb_in_thermodynamic_gradients = .false. mixing_length_alpha = 1.5d0 - alpha_TDC_DAMP = 1d0 - alpha_TDC_DAMPR = 0d0 - alpha_TDC_PtdVdt = 0d0 - alpha_TDC_DAMPM = 0.25d0 + TDC_alpha_D = 1d0 + TDC_alpha_R = 0d0 + TDC_alpha_Pt = 0d0 + TDC_alpha_M = 0.25d0 + use_TDC_enthalpy_flux_limiter = .false. steps_before_use_TDC = 0 mlt_Pturb_factor = 0d0 ! use 2d0 = 1d0 in TDC. time lagged. alt_scale_height_flag = .false. ! ignore eggleton in the core. TDC_num_innermost_cells_forced_nonturbulent = 2 ! for envelope models only. - TDC_use_density_form_for_eddy_viscosity = .false. ! always true for u_flag. mlt_make_surface_no_mixing = .false. diff --git a/star/dev_cases_TDC_Pulsation/dev_TDC_Cepheid_6M/profile_columns.list b/star/dev_cases_TDC_Pulsation/dev_TDC_Cepheid_6M/profile_columns.list index 1078b36b2..4db75a670 100644 --- a/star/dev_cases_TDC_Pulsation/dev_TDC_Cepheid_6M/profile_columns.list +++ b/star/dev_cases_TDC_Pulsation/dev_TDC_Cepheid_6M/profile_columns.list @@ -418,6 +418,7 @@ !cell_internal_energy_fraction !cell_internal_energy_fraction_start !cell_specific_PE + dwork_dm !log_cell_ie_div_star_ie !log_cell_specific_IE 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/other/other_mlt_results.f90 b/star/other/other_mlt_results.f90 index 358d9f32b..28a23107b 100644 --- a/star/other/other_mlt_results.f90 +++ b/star/other/other_mlt_results.f90 @@ -31,7 +31,7 @@ subroutine null_other_mlt_results(id, k, MLT_option, & ! NOTE: k=0 is a valid a r, L, T, P, opacity, rho, chiRho, chiT, Cp, gradr, grada, scale_height, & iso, XH1, cgrav, m, gradL_composition_term, mixing_length_alpha, & alpha_semiconvection, thermohaline_coeff, & - mixing_type, gradT, Y_face, conv_vel, D, Gamma, ierr) + mixing_type, gradT, Y_face, conv_vel, D, Gamma, energy, ierr) use const_def, only: dp use auto_diff use star_def @@ -39,7 +39,7 @@ subroutine null_other_mlt_results(id, k, MLT_option, & ! NOTE: k=0 is a valid a integer, intent(in) :: k character(len=*), intent(in) :: MLT_option type(auto_diff_real_star_order1), intent(in) :: & - r, L, T, P, opacity, rho, chiRho, chiT, Cp, gradr, grada, scale_height + r, L, T, P, opacity, rho, chiRho, chiT, Cp, gradr, grada, scale_height, energy integer, intent(in) :: iso real(dp), intent(in) :: & XH1, cgrav, m, gradL_composition_term, & diff --git a/star/private/auto_diff_support.f90 b/star/private/auto_diff_support.f90 index ac5e40539..b5665322f 100644 --- a/star/private/auto_diff_support.f90 +++ b/star/private/auto_diff_support.f90 @@ -1235,6 +1235,38 @@ 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 = 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/ctrls_io.f90 b/star/private/ctrls_io.f90 index 97832f197..5f9e6e525 100644 --- a/star/private/ctrls_io.f90 +++ b/star/private/ctrls_io.f90 @@ -105,10 +105,13 @@ module ctrls_io mixing_D_limit_for_log, trace_mass_location, min_tau_for_max_abs_v_location, & min_q_for_inner_mach1_location, max_q_for_outer_mach1_location, & conv_core_gap_dq_limit, & - alpha_TDC_DAMP, alpha_TDC_DAMPR, alpha_TDC_PtdVdt, alpha_TDC_DAMPM, & - alpha_TDC_C, alpha_TDC_S, TDC_use_density_form_for_eddy_viscosity, & - TDC_num_innermost_cells_forced_nonturbulent, include_mlt_Pturb_in_thermodynamic_gradients, & - include_mlt_corr_to_TDC, TDC_include_eturb_in_energy_equation, & + TDC_alpha_D, TDC_alpha_R, TDC_alpha_Pt, TDC_alpha_M, & + TDC_alpha_C, TDC_alpha_S, & + TDC_alpha_M_use_explicit_mlt_vc_in_momentum_equation, & + TDC_use_density_form_for_eddy_viscosity, & + TDC_num_innermost_cells_forced_nonturbulent, TDC_num_outermost_cells_forced_nonturbulent, & + include_mlt_Pturb_in_thermodynamic_gradients, & + include_mlt_corr_to_TDC, use_TDC_enthalpy_flux_limiter, TDC_include_eturb_in_energy_equation, & use_rsp_form_of_scale_height, include_mlt_in_velocity_time_centering, & TDC_hydro_use_mass_interp_face_values, TDC_hydro_nz, TDC_hydro_nz_outer, TDC_hydro_T_anchor, TDC_hydro_dq_1_factor, & @@ -259,7 +262,7 @@ module ctrls_io merge_amr_ignore_core_cells, merge_amr_logT_for_ignore_core_cells, & split_amr_ignore_core_cells, split_amr_logT_for_ignore_core_cells, & merge_amr_du_div_cs_limit_only_for_compression, split_merge_amr_avoid_repeated_remesh, split_merge_amr_r_core_cm, & - split_merge_amr_dq_min, split_merge_amr_dq_max, split_merge_amr_max_iters, trace_split_merge_amr, equal_split_density_amr, & + split_merge_amr_dq_min, split_merge_amr_dq_max, split_merge_amr_max_iters, trace_split_merge_amr, equal_split_density_amr, use_hydro_merge_limits_in_mesh_plan, & ! nuclear reaction parameters screening_mode, default_net_name, net_logTcut_lo, net_logTcut_lim, & @@ -1611,6 +1614,7 @@ subroutine store_controls(s, ierr) s% split_merge_amr_max_iters = split_merge_amr_max_iters s% trace_split_merge_amr = trace_split_merge_amr s% equal_split_density_amr = equal_split_density_amr + s% use_hydro_merge_limits_in_mesh_plan = use_hydro_merge_limits_in_mesh_plan ! nuclear reaction parameters s% screening_mode = screening_mode @@ -2079,16 +2083,19 @@ subroutine store_controls(s, ierr) s% max_safe_logT_for_rates = max_safe_logT_for_rates s% eps_mdot_leak_frac_factor = eps_mdot_leak_frac_factor - s% alpha_TDC_DAMP = alpha_TDC_DAMP - s% alpha_TDC_DAMPR = alpha_TDC_DAMPR - s% alpha_TDC_PtdVdt = alpha_TDC_PtdVdt - s% alpha_TDC_DAMPM = alpha_TDC_DAMPM - s% alpha_TDC_C = alpha_TDC_C - s% alpha_TDC_S = alpha_TDC_S + s% TDC_alpha_D = TDC_alpha_D + s% TDC_alpha_R = TDC_alpha_R + s% TDC_alpha_Pt = TDC_alpha_Pt + s% TDC_alpha_M = TDC_alpha_M + s% TDC_alpha_C = TDC_alpha_C + s% TDC_alpha_S = TDC_alpha_S + s% TDC_alpha_M_use_explicit_mlt_vc_in_momentum_equation = TDC_alpha_M_use_explicit_mlt_vc_in_momentum_equation s% TDC_use_density_form_for_eddy_viscosity = TDC_use_density_form_for_eddy_viscosity s% TDC_num_innermost_cells_forced_nonturbulent = TDC_num_innermost_cells_forced_nonturbulent + s% TDC_num_outermost_cells_forced_nonturbulent = TDC_num_outermost_cells_forced_nonturbulent s% include_mlt_Pturb_in_thermodynamic_gradients = include_mlt_Pturb_in_thermodynamic_gradients s% include_mlt_corr_to_TDC = include_mlt_corr_to_TDC + s% use_TDC_enthalpy_flux_limiter = use_TDC_enthalpy_flux_limiter s% TDC_include_eturb_in_energy_equation = TDC_include_eturb_in_energy_equation s% use_rsp_form_of_scale_height = use_rsp_form_of_scale_height s% include_mlt_in_velocity_time_centering = include_mlt_in_velocity_time_centering @@ -3322,7 +3329,7 @@ subroutine set_controls_for_writing(s, ierr) split_merge_amr_max_iters = s% split_merge_amr_max_iters trace_split_merge_amr = s% trace_split_merge_amr equal_split_density_amr = s% equal_split_density_amr - + use_hydro_merge_limits_in_mesh_plan = s% use_hydro_merge_limits_in_mesh_plan ! nuclear reaction parameters screening_mode = s% screening_mode default_net_name = s% default_net_name @@ -3787,16 +3794,19 @@ subroutine set_controls_for_writing(s, ierr) max_safe_logT_for_rates = s% max_safe_logT_for_rates eps_mdot_leak_frac_factor = s% eps_mdot_leak_frac_factor - alpha_TDC_DAMP = s% alpha_TDC_DAMP - alpha_TDC_DAMPR = s% alpha_TDC_DAMPR - alpha_TDC_PtdVdt = s% alpha_TDC_PtdVdt - alpha_TDC_DAMPM = s% alpha_TDC_DAMPM - alpha_TDC_C = s% alpha_TDC_C - alpha_TDC_S = s% alpha_TDC_S + TDC_alpha_D = s% TDC_alpha_D + TDC_alpha_R = s% TDC_alpha_R + TDC_alpha_Pt = s% TDC_alpha_Pt + TDC_alpha_M = s% TDC_alpha_M + TDC_alpha_C = s% TDC_alpha_C + TDC_alpha_S = s% TDC_alpha_S + TDC_alpha_M_use_explicit_mlt_vc_in_momentum_equation = s% TDC_alpha_M_use_explicit_mlt_vc_in_momentum_equation TDC_use_density_form_for_eddy_viscosity = s% TDC_use_density_form_for_eddy_viscosity TDC_num_innermost_cells_forced_nonturbulent = s% TDC_num_innermost_cells_forced_nonturbulent + TDC_num_outermost_cells_forced_nonturbulent = s% TDC_num_outermost_cells_forced_nonturbulent include_mlt_Pturb_in_thermodynamic_gradients = s% include_mlt_Pturb_in_thermodynamic_gradients include_mlt_corr_to_TDC = s% include_mlt_corr_to_TDC + use_TDC_enthalpy_flux_limiter = s% use_TDC_enthalpy_flux_limiter TDC_include_eturb_in_energy_equation = s% TDC_include_eturb_in_energy_equation use_rsp_form_of_scale_height = s% use_rsp_form_of_scale_height include_mlt_in_velocity_time_centering = s% include_mlt_in_velocity_time_centering diff --git a/star/private/hydro_energy.f90 b/star/private/hydro_energy.f90 index ef4cceb36..b22d3c801 100644 --- a/star/private/hydro_energy.f90 +++ b/star/private/hydro_energy.f90 @@ -216,8 +216,7 @@ 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 + use tdc_hydro, only: compute_tdc_Eq_div_w_face real(dp) :: alfa, beta integer, intent(out) :: ierr type(auto_diff_real_star_order1) :: & @@ -266,15 +265,13 @@ subroutine setup_sources_and_others(ierr) ! sources_ad, others_ad if (s% RSP2_flag) then 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. - !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) + else if (s% TDC_alpha_M >0d0 .and. s% MLT_option == 'TDC' .and. & + s% TDC_include_eturb_in_energy_equation .and. (s% v_flag .or. s% u_flag)) then + 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,29 +352,34 @@ 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 d_turbulent_energy_dt_ad = (wrap_etrb_00(s,k) - get_etrb_start(s,k))/dt - else if (s% mlt_vc_old(k) > 0d0 .and. s% MLT_option == 'TDC' .and. & - s% TDC_include_eturb_in_energy_equation) then + else if (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) - end if - d_turbulent_energy_dt_ad = (TDC_eturb_cell - TDC_eturb_cell_start) / dt + if (k < s% nz) then + if (s% okay_to_set_mlt_vc) then ! have mlt_vc_old + TDC_eturb_cell_start = 0.75d0*(pow2(s% mlt_vc_old(k)) + & + pow2(s% mlt_vc_old(k+1))) + else + TDC_eturb_cell_start = 0d0 + end if + TDC_eturb_cell = 0.75d0*(pow2(s% mlt_vc_ad(k)) + & + pow2(shift_p1(s% mlt_vc_ad(k+1)))) + else ! center cell averaged with 0 for inner face + if (s% okay_to_set_mlt_vc) then ! have mlt_vc_old + TDC_eturb_cell_start = 0.75d0*pow2(s% mlt_vc_old(k)) + else + TDC_eturb_cell_start = 0d0 + end if + TDC_eturb_cell = 0.75d0*pow2(s% mlt_vc_ad(k)) + end if + d_turbulent_energy_dt_ad = (TDC_eturb_cell - TDC_eturb_cell_start)/dt else d_turbulent_energy_dt_ad = 0d0 end if diff --git a/star/private/hydro_momentum.f90 b/star/private/hydro_momentum.f90 index 9ff009e58..ef8e620a3 100644 --- a/star/private/hydro_momentum.f90 +++ b/star/private/hydro_momentum.f90 @@ -430,7 +430,7 @@ subroutine expected_non_HSE_term( & 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 - else if (s% alpha_TDC_DampM > 0 .and. s% MLT_option == 'TDC') then ! Uq(k) is turbulent viscosity drag at face k + else if (s% TDC_alpha_M > 0 .and. s% MLT_option == 'TDC') then ! Uq(k) is turbulent viscosity drag at face k Uq_ad = compute_tdc_Uq_face(s, k, ierr) if (ierr /= 0) return end if diff --git a/star/private/hydro_riemann.f90 b/star/private/hydro_riemann.f90 index 5a904ef37..ca7e83578 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%TDC_alpha_M > 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..32fc189ac 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% TDC_alpha_M > 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/mesh_plan.f90 b/star/private/mesh_plan.f90 index 17a40bf78..a55247825 100644 --- a/star/private/mesh_plan.f90 +++ b/star/private/mesh_plan.f90 @@ -602,9 +602,9 @@ subroutine pick_new_points(s, ierr) type (star_info), pointer :: s integer, intent(out) :: ierr - logical :: dbg, force_merge_with_one_more + logical :: dbg, force_merge_with_one_more, du_div_cs_limit_flag real(dp) :: maxval_delta_xa, next_dq_max, beta_limit, & - remaining_dq_old, min_dr + remaining_dq_old, min_dr, abs_du_div_cs integer :: kk, k_old_init, k_old_next, k_old_next_max, j00, jm1, i, max_merge include 'formats' @@ -616,6 +616,7 @@ subroutine pick_new_points(s, ierr) k_new = 1 xq_new(1) = 0 min_dr = s% mesh_min_dr_div_dRstar*(s% r(1) - s% R_center) + abs_du_div_cs = 0d0 do ! pick next point location @@ -738,6 +739,53 @@ subroutine pick_new_points(s, ierr) k_old_next_max = min(nz_old, k_old + max_merge) k_old_next = k_old_next_max ! will cut this back as necessary do kk=k_old+1,k_old_next_max + + if (s% use_hydro_merge_limits_in_mesh_plan) then ! limit merges over steep velocity gradients + ! begin hydro check_merge_limits section + du_div_cs_limit_flag = .false. + if (.not. s% merge_amr_du_div_cs_limit_only_for_compression) then + du_div_cs_limit_flag = .true. + else if (associated(s% v)) then + ! Only set flag for compressive flow across interface (kk-1, kk) + if (kk <= nz_old .and. s% v(kk)*pow2(s% r(kk)) > s% v(kk-1)*pow2(s% r(kk-1))) du_div_cs_limit_flag = .true. + if (.not. du_div_cs_limit_flag .and. kk-1 > 1) then + if (s% v(kk-1)*pow2(s% r(kk-1)) > s% v(kk-2)*pow2(s% r(kk-2))) du_div_cs_limit_flag = .true. + end if + end if + + if (du_div_cs_limit_flag .and. associated(s% v)) then + if (kk-1 == 1) then + abs_du_div_cs = abs(s% v(1) - s% v(2)) / s% csound(1) + else if (kk == nz_old) then + abs_du_div_cs = abs(s% v(nz_old-1) - s% v(nz_old)) / s% csound(nz_old) + else + abs_du_div_cs = abs(s% v(kk-1) - s% v(kk)) / s% csound(kk-1) + end if + else + abs_du_div_cs = 0.0_dp + end if + + if (du_div_cs_limit_flag) then + if (.not. s% merge_amr_inhibit_at_jumps) then + if (abs_du_div_cs > s% merge_amr_max_abs_du_div_cs) then + ! Shock jump too large, so block merge at this interface + k_old_next = kk-1 + exit + end if + else ! inhibit_at_jumps = .true. + if (abs_du_div_cs > s% merge_amr_max_abs_du_div_cs) then + if (dq_old(k_old) >= min_dq) then + ! Shock is large and zone is not extremely small, so inhibit merge + k_old_next = kk-1 + exit + end if + ! else: if zone is very small, allow merge despite the shock + end if + end if + end if + end if + ! end hydro check_merge_limits section + maxval_delta_xa = maxval(abs(s% xa(:,kk)-s% xa(:,kk-1))) j00 = maxloc(s% xa(:,kk),dim=1) jm1 = maxloc(s% xa(:,kk-1),dim=1) diff --git a/star/private/profile_getval.f90 b/star/private/profile_getval.f90 index 946d17439..ebc2b1f2d 100644 --- a/star/private/profile_getval.f90 +++ b/star/private/profile_getval.f90 @@ -2026,6 +2026,10 @@ subroutine getval_for_profile(s, c, k, val, int_flag, int_val) end if val = dble(int_val) int_flag = .true. + case (p_dwork_dm) + ! differential work per unit mass per unit mass*time dW/dm + ! W = dwork_dm*dm*dt + val = s% dwork_dm(k) ! returns (dw/dt)/dm case (p_cell_specific_IE) val = s% energy(k) diff --git a/star/private/read_model.f90 b/star/private/read_model.f90 index 292edad9a..7a1014cbf 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% TDC_alpha_M > 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/star_profile_def.f90 b/star/private/star_profile_def.f90 index 1489ffd2d..7da296819 100644 --- a/star/private/star_profile_def.f90 +++ b/star/private/star_profile_def.f90 @@ -572,8 +572,9 @@ module star_profile_def integer, parameter :: p_total_energy_sign = p_cs_at_cell_bdy + 1 integer, parameter :: p_total_energy = p_total_energy_sign + 1 + integer, parameter :: p_dwork_dm = p_total_energy + 1 - integer, parameter :: p_Ptrb = p_total_energy + 1 + integer, parameter :: p_Ptrb = p_dwork_dm + 1 integer, parameter :: p_log_Ptrb = p_Ptrb + 1 integer, parameter :: p_log_w = p_log_Ptrb + 1 integer, parameter :: p_w = p_log_w + 1 @@ -1249,6 +1250,7 @@ subroutine profile_column_names_init(ierr) profile_column_name(p_total_energy_sign) = 'total_energy_sign' profile_column_name(p_total_energy) = 'total_energy' + profile_column_name(p_dwork_dm) = 'dwork_dm' profile_column_name(p_Ptrb) = 'Ptrb' profile_column_name(p_log_Ptrb) = 'log_Ptrb' diff --git a/star/private/star_utils.f90 b/star/private/star_utils.f90 index 774340c5c..dc3132179 100644 --- a/star/private/star_utils.f90 +++ b/star/private/star_utils.f90 @@ -86,9 +86,10 @@ module star_utils public :: get_kap_face public :: get_rho_face public :: get_rho_start_face - public :: get_chirho_face - public :: get_chit_face - public :: get_t_face + public :: get_e_face + public :: get_ChiRho_face + public :: get_ChiT_face + public :: get_T_face public :: get_Peos_face public :: get_Peos_face_val public :: get_cp_face @@ -2490,15 +2491,13 @@ subroutine eval_deltaM_total_energy_integrals( & cell1 = dm*pow2(s% w(k)) cell_total = cell_total + cell1 total_turbulent_energy = total_turbulent_energy + cell1 - else if (.not. s% RSP2_flag .and. s% mlt_vc(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 = pow2(s% mlt_vc(k)/sqrt_2_div_3) - else - call get_face_weights(s, k, alfa, beta) - TDC_eturb_cell = alfa*pow2(s% mlt_vc(k)/sqrt_2_div_3) + & - beta*pow2(s% mlt_vc(k-1)/sqrt_2_div_3) + else if ( s% MLT_option == 'TDC' .and. & + s% TDC_include_eturb_in_energy_equation) then ! needs corrected s% mlt_vc(k) >= 0 causes failures. + if (k < s% nz) then + TDC_eturb_cell = 0.75d0*(pow2(s% mlt_vc(k)) + & + pow2(s% mlt_vc(k+1))) + else ! k == s% nz + TDC_eturb_cell = 0.75d0*pow2(s% mlt_vc(k)) end if cell1 = dm*TDC_eturb_cell cell_total = cell_total + cell1 @@ -3870,6 +3869,18 @@ function get_Peos_face_val(s,k) result(Peos_face) Peos_face = alfa*s% Peos(k) + beta*s% Peos(k-1) end function get_Peos_face_val + function get_e_face(s,k) result(e_face) ! specific energy on face + type (star_info), pointer :: s + integer, intent(in) :: k + type(auto_diff_real_star_order1) :: e_face + real(dp) :: alfa, beta + if (k == 1) then + e_face = wrap_e_00(s,k) + return + end if + call get_face_weights(s, k, alfa, beta) + e_face = alfa*wrap_e_00(s,k) + beta*wrap_e_m1(s,k) + end function get_e_face function get_Cp_face(s,k) result(Cp_face) type (star_info), pointer :: s diff --git a/star/private/tdc_hydro.f90 b/star/private/tdc_hydro.f90 index 34afd1b50..f0b3145d9 100644 --- a/star/private/tdc_hydro.f90 +++ b/star/private/tdc_hydro.f90 @@ -31,8 +31,8 @@ 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 + compute_tdc_Uq_face, compute_tdc_Eq_div_w_face, & + 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,15 @@ 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) + x = compute_Chi_div_w_face(s, k, op_err) ! Sets Chi_face if (op_err /= 0) ierr = op_err - x = compute_tdc_Eq_cell(s, k, op_err) + x = compute_tdc_Eq_div_w_face(s, k, op_err) ! Sets Eq_face 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,21 +86,14 @@ 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 ! face_value = alfa*cell_value(k) + beta*cell_value(k-1) - if (k == 1) call mesa_error(__FILE__, __LINE__, 'bad k==1 for get_RSP2_alfa_beta_face_weights') + if (k == 1) call mesa_error(__FILE__, __LINE__, 'bad k==1 for get_TDC_alfa_beta_face_weights') if (s%TDC_hydro_use_mass_interp_face_values) then alfa = s%dq(k - 1)/(s%dq(k - 1) + s%dq(k)) beta = 1d0 - alfa @@ -95,145 +101,9 @@ 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 - - function compute_d_v_div_r(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, term1, term2 - logical :: dbg - include 'formats' - ierr = 0 - dbg = .false. - v_00 = wrap_v_00(s, k) - v_p1 = wrap_v_p1(s, k) - r_00 = wrap_r_00(s, k) - r_p1 = 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 - - ! 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 (*, *) 'd_v_div_r:', d_v_div_r%val - end if - end function compute_d_v_div_r - - 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) :: 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) - 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)) - 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) :: 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) - 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)) - end function compute_rho_form_of_d_v_div_r_opt_time_center - - - function compute_rho_form_of_d_v_div_r_face(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) :: r_face, rho_face, v_face, dlnrho_dt - real(dp) :: dm_bar - 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ρ - - if (k >= 2) then - dm_bar = 0.5d0*(s% dm(k) + s% dm(k-1)) - else - dm_bar = 0.5d0*s% dm(k) - end if - - d_v_div_r = -dm_bar/(4d0*pi*rho_face)*(dlnrho_dt/pow3(r_face) + 3d0*v_face/pow4(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 - 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) :: r_face, rho_face, v_face, dlnrho_dt - real(dp) :: dm_bar - 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ρ - - if (k >= 2) then - dm_bar = 0.5d0*(s% dm(k) + s% dm(k-1)) - else - dm_bar = 0.5d0*s% dm(k) - end if - - d_v_div_r = -dm_bar/(4d0*pi*rho_face)*(dlnrho_dt/pow3(r_face) + 3d0*v_face/pow4(r_face)) - end function compute_rho_form_of_d_v_div_r_face_opt_time_center + end subroutine get_TDC_alfa_beta_face_weights - function compute_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 - 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 - include 'formats' - ierr = 0 - v_00 = wrap_opt_time_center_v_00(s, k) - v_p1 = wrap_opt_time_center_v_p1(s, k) - r_00 = wrap_opt_time_center_r_00(s, k) - r_p1 = 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 - end function compute_d_v_div_r_opt_time_center - function wrap_Hp_cell(s, k) result(Hp_cell) ! cm , different than rsp2 type(star_info), pointer :: s integer, intent(in) :: k @@ -241,7 +111,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)) @@ -277,7 +147,8 @@ function Hp_cell_for_Chi(s, k, ierr) result(Hp_cell) ! cm end if end function Hp_cell_for_Chi - function compute_Chi_cell(s, k, ierr) result(Chi_cell) + ! this function is only called internally in TDC_Uq_face, and for v_flag only. + function compute_Chi_cell(s, k, ierr) result(Chi_cell) ! does not update s% Chi or Chi_ad ! eddy viscosity energy (Kuhfuss 1986) [erg] type(star_info), pointer :: s integer, intent(in) :: k @@ -293,22 +164,19 @@ function compute_Chi_cell(s, k, ierr) result(Chi_cell) ! check where we are getting alfam from. if (s%MLT_option == 'TDC' .and. .not. s%RSP2_flag) then - ALFAM_ALFA = s%alpha_TDC_DampM*s%mixing_length_alpha + ALFAM_ALFA = s%TDC_alpha_M*s%mixing_length_alpha else ! this is for safety, but probably is never called. ALFAM_ALFA = 0d0 end if if (ALFAM_ALFA == 0d0 .or. & + k <= s% TDC_num_outermost_cells_forced_nonturbulent .or. & k > s% nz - s% TDC_num_innermost_cells_forced_nonturbulent) then Chi_cell = 0d0 - if (k >= 1 .and. k <= s%nz) then - s%Chi(k) = 0d0 - s%Chi_ad(k) = 0d0 - end if else Hp_cell = Hp_cell_for_Chi(s, k, ierr) if (ierr /= 0) return - if (s%u_flag .or. s%TDC_use_density_form_for_eddy_viscosity) then + if (s%TDC_use_density_form_for_eddy_viscosity) then ! new density derivative term d_v_div_r = compute_rho_form_of_d_v_div_r(s, k, ierr) else @@ -317,10 +185,20 @@ 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% TDC_alpha_M_use_explicit_mlt_vc_in_momentum_equation) 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% TDC_alpha_M_use_explicit_mlt_vc_in_momentum_equation) 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) @@ -334,8 +212,9 @@ function compute_Chi_cell(s, k, ierr) result(Chi_cell) ! = erg end if - s%Chi(k) = Chi_cell%val - s%Chi_ad(k) = Chi_cell + ! this is set in Chi_div_w_face + !s%Chi(k) = Chi_cell%val + !s%Chi_ad(k) = Chi_cell if (dbg .and. k == -100) then write (*, *) ' s% ALFAM_ALFA', ALFAM_ALFA @@ -351,74 +230,6 @@ function compute_Chi_cell(s, k, ierr) result(Chi_cell) end if end function compute_Chi_cell - function compute_tdc_Eq_cell(s, k, ierr) result(Eq_cell) ! erg g^-1 s^-1 - type(star_info), pointer :: s - integer, intent(in) :: k - type(auto_diff_real_star_order1) :: Eq_cell - integer, intent(out) :: ierr - type(auto_diff_real_star_order1) :: d_v_div_r, Chi_cell - include 'formats' - ierr = 0 - if (s%mixing_length_alpha == 0d0 .or. & - k > s% nz - s% TDC_num_innermost_cells_forced_nonturbulent) then - Eq_cell = 0d0 - if (k >= 1 .and. k <= s%nz) s%Eq_ad(k) = 0d0 - else - Chi_cell = s%Chi_ad(k) ! compute_Chi_cell(s,k,ierr) - if (ierr /= 0) return - - if (s%u_flag .or. s%TDC_use_density_form_for_eddy_viscosity) then - ! new density derivative term - d_v_div_r = compute_rho_form_of_d_v_div_r_opt_time_center(s, k, ierr) - else - d_v_div_r = compute_d_v_div_r_opt_time_center(s, k, ierr) - end if - - if (ierr /= 0) return - Eq_cell = 4d0*pi*Chi_cell*d_v_div_r/s%dm(k) ! erg s^-1 g^-1 - end if - s%Eq(k) = Eq_cell%val - s%Eq_ad(k) = Eq_cell - end function compute_tdc_Eq_cell - - function compute_tdc_Uq_face(s, k, ierr) result(Uq_face) ! cm s^-2, acceleration - type(star_info), pointer :: s - integer, intent(in) :: k - type(auto_diff_real_star_order1) :: Uq_face - integer, intent(out) :: ierr - type(auto_diff_real_star_order1) :: Chi_00, Chi_m1, r_00 - 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_face = 0d0 - else - r_00 = wrap_opt_time_center_r_00(s, k) - - ! which do we adopt? - Chi_00 = compute_Chi_cell(s, k, ierr) ! s% Chi_ad(k) XXX - !Chi_00 = s% Chi_ad(k) ! compute_Chi_cell(s,k,ierr) - - if (k > 1) then - 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 - end if - Uq_face = 4d0*pi*(Chi_m1 - Chi_00)/(r_00*s%dm_bar(k)) - - if (k == -56) then - write (*, 3) 'RSP2 Uq chi_m1 chi_00 r', k, s%solver_iter, & - Uq_face%val, Chi_m1%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_face%val - end function compute_tdc_Uq_face - ! face centered variables for tdc update below function compute_Chi_div_w_face(s, k, ierr) result(Chi_face) ! eddy viscosity energy (Kuhfuss 1986) [erg] @@ -436,7 +247,7 @@ function compute_Chi_div_w_face(s, k, ierr) result(Chi_face) ! check where we are getting alfam from. if (s%MLT_option == 'TDC' .and. .not. s%RSP2_flag) then - ALFAM_ALFA = s%alpha_TDC_DampM*s%mixing_length_alpha + ALFAM_ALFA = s%TDC_alpha_M*s%mixing_length_alpha else ! this is for safety, but probably is never called. ALFAM_ALFA = 0d0 end if @@ -444,15 +255,11 @@ 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 if (s%u_flag .or. s%TDC_use_density_form_for_eddy_viscosity) then - ! new density derivative term + ! new density derivative form d_v_div_r = compute_rho_form_of_d_v_div_r_face(s, k, ierr) else d_v_div_r = compute_d_v_div_r_face(s, k, ierr) @@ -471,36 +278,42 @@ function compute_Chi_div_w_face(s, k, ierr) result(Chi_face) !r_p1 = wrap_r_p1(s, k) r6_face = pow6(r_00) !0.5d0*(pow6(r_00) + pow6(r_p1)) Chi_face = f*rho2*r6_face*d_v_div_r*Hp_face!*w_00 - ! units = g^-1 cm s^-1 g^2 cm^-6 cm^6 s^-1 cm - ! = g cm^2 s^-2 - ! = erg ! * (s / cm) - > [erg] * 1/w00 + ! units = g^-1 cm s^-1 g^2 cm^-6 cm^6 s^-1 cm * [s/cm] ! [1/w_00] = [s/cm] + ! = g cm^2 s^-2 * [s/cm] + ! = erg ! * [s / cm] - > [erg] * [s/cm] end if - !s%Chi(k) = Chi_face%val - !s%Chi_ad(k) = Chi_face - - if (dbg .and. k == -100) then - write (*, *) ' s% ALFAM_ALFA', ALFAM_ALFA - write (*, *) 'Hp_face', Hp_face%val - write (*, *) 'd_v_div_r', d_v_div_r%val - write (*, *) ' f', f - write (*, *) 'w_00', w_00%val - write (*, *) 'd_00 ', d_00%val - write (*, *) 'rho2 ', rho2%val - write (*, *) 'r_00', r_00%val - write (*, *) 'r_p1 ', r_p1%val - write (*, *) 'r6_cell', r6_face%val - end if - end function compute_Chi_div_w_face + ! Chi_cell does not set Chi, we store Chi_face in s% Chi and s% Chi_ad + if (s% okay_to_set_mlt_vc .and. & + s% TDC_alpha_M_use_explicit_mlt_vc_in_momentum_equation) 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%Chi(k) = Chi_face%val*w_00%val + s%Chi_ad(k) = Chi_face*w_00 + if (dbg .and. k == -100) then + write (*, *) ' s% ALFAM_ALFA', ALFAM_ALFA + write (*, *) 'Hp_face', Hp_face%val + write (*, *) 'd_v_div_r', d_v_div_r%val + write (*, *) ' f', f + write (*, *) 'w_00', w_00%val + write (*, *) 'd_00 ', d_00%val + write (*, *) 'rho2 ', rho2%val + write (*, *) 'r_00', r_00%val + write (*, *) 'r_p1 ', r_p1%val + write (*, *) 'r6_cell', r6_face%val + end if + end function compute_Chi_div_w_face function compute_tdc_Eq_div_w_face(s, k, ierr) result(Eq_face) ! erg g^-1 s^-1 * (cm^-1 s^1) type(star_info), pointer :: s 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,35 +341,125 @@ 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% TDC_alpha_M_use_explicit_mlt_vc_in_momentum_equation) 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 + ! for v_flag only. face centered Uq for hydro_momentum + function compute_tdc_Uq_face(s, k, ierr) result(Uq_face) !(v_flag only) ! cm s^-2, acceleration type(star_info), pointer :: s integer, intent(in) :: k - type(auto_diff_real_star_order1) :: d_v_div_r + type(auto_diff_real_star_order1) :: Uq_face integer, intent(out) :: ierr - type(auto_diff_real_star_order1) :: v_00, v_p1, r_00, r_p1, term1, term2 - logical :: dbg + type(auto_diff_real_star_order1) :: Chi_00, Chi_m1, r_00 include 'formats' ierr = 0 - dbg = .false. - - if (k >1) then - v_00 = 0.5d0 * (wrap_v_00(s, k) + wrap_v_m1(s, k)) + if (s%mixing_length_alpha == 0d0 .or. & + k <= s% TDC_num_outermost_cells_forced_nonturbulent .or. & + k > s%nz - s% TDC_num_innermost_cells_forced_nonturbulent) then + Uq_face = 0d0 else - v_00 = 0.5d0 * wrap_v_00(s, k) - end if + r_00 = wrap_opt_time_center_r_00(s, k) - v_p1 = 0.5d0 * (wrap_v_00(s, k) + wrap_v_p1(s, k)) + ! which do we adopt? + Chi_00 = compute_Chi_cell(s, k, ierr) ! s% Chi_ad(k) XXX - if (k >1) then - r_00 = 0.5d0 * (wrap_r_00(s, k) + wrap_r_m1(s, k)) + if (k > 1) then + Chi_m1 = shift_m1(compute_Chi_cell(s, k-1, ierr)) + if (ierr /= 0) return + else + Chi_m1 = 0d0 + end if + Uq_face = 4d0*pi*(Chi_m1 - Chi_00)/(r_00*s%dm_bar(k)) + + if (k == -56) then + write (*, 3) 'TDC Uq chi_m1 chi_00 r', k, s%solver_iter, & + Uq_face%val, Chi_m1%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_face%val + end function compute_tdc_Uq_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% TDC_num_outermost_cells_forced_nonturbulent .or. & + k > s%nz - s% TDC_num_innermost_cells_forced_nonturbulent) then + Uq_cell = 0d0 else - r_00 = 0.5d0 * wrap_r_00(s, k) + 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% TDC_alpha_M_use_explicit_mlt_vc_in_momentum_equation) then + w_00 = s% mlt_vc_old(k)/sqrt_2_div_3 + else + w_00 = s% mlt_vc_ad(k)/sqrt_2_div_3 + end if + + Chi_00 = compute_Chi_div_w_face(s, k, ierr) * w_00 + + if (k < s% nz) then + if (s% okay_to_set_mlt_vc .and. & + s% TDC_alpha_M_use_explicit_mlt_vc_in_momentum_equation) then + w_p1 = s% mlt_vc_old(k+1)/sqrt_2_div_3 + else + w_p1 = shift_p1(s% mlt_vc_ad(k+1))/sqrt_2_div_3 + end if + + Chi_p1 = shift_p1(compute_Chi_div_w_face(s, k+1, ierr))*w_p1 + 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 + ! erg g^-1 cm^-1 = g cm^2 s^-2 g^-1 cm^-1 = cm s^-2 [g], acceleration*mass = Force + + if (k == -56) then + write (*, 3) 'TDC 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 - r_p1 = 0.5d0*(wrap_r_00(s, k) + wrap_r_p1(s, k)) + s%Uq(k) = Uq_cell%val/ s% dm(k) + end function compute_tdc_Uq_dm_cell + + +! all the forms of d(v/r)/dr, below + function compute_d_v_div_r(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, term1, term2 + logical :: dbg + include 'formats' + ierr = 0 + dbg = .false. + v_00 = wrap_v_00(s, k) + v_p1 = wrap_v_p1(s, k) + r_00 = wrap_r_00(s, k) + r_p1 = 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 @@ -567,11 +470,9 @@ function compute_d_v_div_r_face(s, k, ierr) result(d_v_div_r) ! s^-1 write (*, *) 'r_00:', r_00%val, 'r_p1:', r_p1%val write (*, *) 'd_v_div_r:', d_v_div_r%val end if - end function compute_d_v_div_r_face - - + end function compute_d_v_div_r - function compute_d_v_div_r_opt_time_center_face(s, k, ierr) result(d_v_div_r) ! s^-1 + function compute_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 type(auto_diff_real_star_order1) :: d_v_div_r @@ -579,26 +480,169 @@ function compute_d_v_div_r_opt_time_center_face(s, k, ierr) result(d_v_div_r) ! type(auto_diff_real_star_order1) :: v_00, v_p1, r_00, r_p1 include 'formats' ierr = 0 + v_00 = wrap_opt_time_center_v_00(s, k) + v_p1 = wrap_opt_time_center_v_p1(s, k) + r_00 = wrap_opt_time_center_r_00(s, k) + r_p1 = 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 + end function compute_d_v_div_r_opt_time_center - if (k >1) then - v_00 = 0.5d0 * (wrap_opt_time_center_v_00(s, k) + wrap_opt_time_center_v_m1(s, k)) + function compute_rho_form_of_d_v_div_r(s, k, ierr) result(d_v_div_r) ! used in Chi_cell + type(star_info), pointer :: s + integer, intent(in) :: k + integer, intent(out) :: ierr + 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) + 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 + + ! density form + d_v_div_r = -dm_cell/(4d0*pi*rho_cell)*(dlnrho_dt/pow3(r_cell) + 3d0*v_cell/pow4(r_cell)) + + ! dm_cell*(1/r * du/dm - U/4/pi/rho/r^4), more sensitive to geometry + !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_face(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) :: r_face, rho_face, v_face, dlnrho_dt + real(dp) :: dm_bar, alfa, beta + ierr = 0 + + r_face = wrap_r_00(s, k) + rho_face = get_rho_face(s, k) + 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 - v_00 = 0.5d0 * wrap_opt_time_center_v_00(s, k) + dm_bar = 0.5d0*s% dm(k) + dlnrho_dt = 0.5d0*wrap_dxh_lnd(s, k)/s%dt ! (∂/∂t)lnρ end if - v_p1 = 0.5d0 * (wrap_opt_time_center_v_00(s, k) + wrap_opt_time_center_v_p1(s, k)) + ! density form + 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), more sensitive to geometry + !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 - if (k >1) then - r_00 = 0.5d0 * (wrap_opt_time_center_r_00(s, k) + wrap_opt_time_center_r_m1(s, k)) + 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 + 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) :: r_face, rho_face, v_face, dlnrho_dt + 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) + 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 - r_00 = 0.5d0 * wrap_opt_time_center_r_00(s, k) + dm_bar = 0.5d0*s% dm(k) + dlnrho_dt = 0.5d0*wrap_dxh_lnd(s, k)/s%dt ! (∂/∂t)lnρ 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 + ! density form + 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), more sensitive to geometry + !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 + + function compute_d_v_div_r_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_m1, r_00, r_m1, term1, term2 + logical :: dbg + include 'formats' + ierr = 0 + dbg = .false. + + 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 + + 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 (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_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_m1, r_00, r_m1, term1, term2 + logical :: dbg + include 'formats' + ierr = 0 + dbg = .false. + + 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)) + + 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_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_opt_time_center_face 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 diff --git a/star/private/turb_support.f90 b/star/private/turb_support.f90 index d35f66e3a..53d5242b0 100644 --- a/star/private/turb_support.f90 +++ b/star/private/turb_support.f90 @@ -69,7 +69,7 @@ subroutine get_gradT(s, MLT_option, & ! used to create models integer, intent(out) :: mixing_type, ierr type(auto_diff_real_star_order1) :: & gradr_ad, grada_ad, scale_height_ad, gradT_ad, Y_face_ad, mlt_vc_ad, D_ad, & - Gamma_ad, r_ad, L_ad, T_ad, P_ad, opacity_ad, rho_ad, dV_ad, chiRho_ad, chiT_ad, Cp_ad + Gamma_ad, r_ad, L_ad, T_ad, P_ad, opacity_ad, rho_ad, dV_ad, chiRho_ad, chiT_ad, Cp_ad, energy_ad ierr = 0 r_ad = r L_ad = L @@ -83,13 +83,14 @@ subroutine get_gradT(s, MLT_option, & ! used to create models Cp_ad = Cp gradr_ad = gradr grada_ad = grada + energy_ad = 0d0 ! correct to a value scale_height_ad = scale_height call Get_results(s, 0, MLT_option, & r_ad, L_ad, T_ad, P_ad, opacity_ad, rho_ad, dV_ad, chiRho_ad, chiT_ad, Cp_ad, & gradr_ad, grada_ad, scale_height_ad, & iso, XH1, cgrav, m, gradL_composition_term, mixing_length_alpha, & s% alpha_semiconvection, s% thermohaline_coeff, & - mixing_type, gradT_ad, Y_face_ad, mlt_vc_ad, D_ad, Gamma_ad, ierr) + mixing_type, gradT_ad, Y_face_ad, mlt_vc_ad, D_ad, Gamma_ad, energy_ad, ierr) gradT = gradT_ad%val Y_face = Y_face_ad%val conv_vel = mlt_vc_ad%val @@ -117,7 +118,8 @@ subroutine do1_mlt_eval( & real(dp) :: cgrav, m, XH1, P_theta, L_theta integer :: iso - type(auto_diff_real_star_order1) :: gradr, r, L, T, P, opacity, rho, dV, chiRho, chiT, Cp, rho_start + type(auto_diff_real_star_order1) :: gradr, r, L, T, P, opacity, rho, dV, & + chiRho, chiT, Cp, rho_start, energy include 'formats' ierr = 0 @@ -158,6 +160,7 @@ subroutine do1_mlt_eval( & chiRho = get_ChiRho_face(s,k) chiT = get_ChiT_face(s,k) Cp = get_Cp_face(s,k) + energy = get_e_face(s,k) iso = s% dominant_iso_for_thermohaline(k) XH1 = s% xa(s% net_iso(ih1),k) @@ -166,7 +169,7 @@ subroutine do1_mlt_eval( & r, L, T, P, opacity, rho, chiRho, chiT, Cp, gradr, grada, scale_height, & iso, XH1, cgrav, m, gradL_composition_term, mixing_length_alpha, & s% alpha_semiconvection, s% thermohaline_coeff, & - mixing_type, gradT, Y_face, mlt_vc, D, Gamma, ierr) + mixing_type, gradT, Y_face, mlt_vc, D, Gamma, energy, ierr) else ! starspot YREC routine if (s% do_starspots) then @@ -177,7 +180,7 @@ subroutine do1_mlt_eval( & r, L, T, P, opacity, rho, dV, chiRho, chiT, Cp, gradr, grada, scale_height, & iso, XH1, cgrav, m, gradL_composition_term, mixing_length_alpha, & s% alpha_semiconvection, s% thermohaline_coeff, & - mixing_type, gradT, Y_face, mlt_vc, D, Gamma, ierr) + mixing_type, gradT, Y_face, mlt_vc, D, Gamma, energy, ierr) end if end subroutine do1_mlt_eval @@ -187,14 +190,14 @@ subroutine Get_results(s, k, MLT_option, & ! NOTE: k=0 is a valid arg r, L, T, P, opacity, rho, dV, chiRho, chiT, Cp, gradr, grada, scale_height, & iso, XH1, cgrav, m, gradL_composition_term, mixing_length_alpha, & alpha_semiconvection, thermohaline_coeff, & - mixing_type, gradT, Y_face, conv_vel, D, Gamma, ierr) + mixing_type, gradT, Y_face, conv_vel, D, Gamma, energy, ierr) use star_utils - use tdc_hydro, only: compute_tdc_Eq_cell, compute_tdc_Eq_div_w_face + use tdc_hydro, only: compute_tdc_Eq_div_w_face type (star_info), pointer :: s integer, intent(in) :: k character (len=*), intent(in) :: MLT_option type(auto_diff_real_star_order1), intent(in) :: & - r, L, T, P, opacity, rho, dV, chiRho, chiT, Cp, gradr, grada, scale_height + r, L, T, P, opacity, rho, dV, chiRho, chiT, Cp, gradr, grada, scale_height, energy integer, intent(in) :: iso real(dp), intent(in) :: & XH1, cgrav, m, gradL_composition_term, & @@ -224,11 +227,9 @@ subroutine Get_results(s, k, MLT_option, & ! NOTE: k=0 is a valid arg ! Pre-calculate some things. Eq_div_w = 0d0 if ((s% v_flag .or. s% u_flag) .and. k > 0 ) then ! only include Eq_div_w if v_flag or u_flag is true. - if (using_TDC .and. s% alpha_TDC_DampM > 0) then - if (s% mlt_vc(k) > 0) then ! calculate using mlt_vc from current timestep. - check_Eq = compute_tdc_Eq_div_w_face(s, k, ierr) - Eq_div_w = check_Eq - end if + if (using_TDC .and. s% TDC_alpha_M > 0) then + check_Eq = compute_tdc_Eq_div_w_face(s, k, ierr) + Eq_div_w = check_Eq end if end if @@ -316,11 +317,11 @@ subroutine Get_results(s, k, MLT_option, & ! NOTE: k=0 is a valid arg end if call set_TDC(& - conv_vel_start, mixing_length_alpha, s% alpha_TDC_DAMP, s%alpha_TDC_DAMPR, s%alpha_TDC_PtdVdt, & + conv_vel_start, mixing_length_alpha, s%TDC_alpha_D, s%TDC_alpha_R, s%TDC_alpha_Pt, & s%dt, cgrav, m, report, & mixing_type, scale, chiT, chiRho, gradr, r, Ptot, T, rho, dV, Cp, opacity, & scale_height, gradL, grada, conv_vel, D, Y_face, gradT, s%tdc_num_iters(k), max_conv_vel, & - Eq_div_w, grav, s% include_mlt_corr_to_TDC, s% alpha_TDC_C, s% alpha_TDC_S, ierr) + Eq_div_w, grav, s% include_mlt_corr_to_TDC, s% TDC_alpha_C, s% TDC_alpha_S, s% use_TDC_enthalpy_flux_limiter, energy, ierr) s% dvc_dt_TDC(k) = (conv_vel%val - conv_vel_start) / s%dt if (ierr /= 0) then @@ -336,11 +337,11 @@ subroutine Get_results(s, k, MLT_option, & ! NOTE: k=0 is a valid arg call set_superad_reduction if (Gamma_factor > 1d0) then call set_TDC(& - conv_vel_start, mixing_length_alpha, s% alpha_TDC_DAMP, s%alpha_TDC_DAMPR, s%alpha_TDC_PtdVdt, & + conv_vel_start, mixing_length_alpha, s%TDC_alpha_D, s%TDC_alpha_R, s%TDC_alpha_Pt, & s%dt, cgrav, m, report, & mixing_type, scale, chiT, chiRho, gradr_scaled, r, Ptot, T, rho, dV, Cp, opacity, & scale_height, gradL, grada, conv_vel, D, Y_face, gradT, s%tdc_num_iters(k), max_conv_vel, & - Eq_div_w, grav, s% include_mlt_corr_to_TDC, s% alpha_TDC_C, s% alpha_TDC_S, ierr) + Eq_div_w, grav, s% include_mlt_corr_to_TDC, s% TDC_alpha_C, s% TDC_alpha_S, s% use_TDC_enthalpy_flux_limiter, energy, 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' @@ -419,8 +420,8 @@ subroutine Get_results(s, k, MLT_option, & ! NOTE: k=0 is a valid arg ! Prevent convection near center of model for MLT or TDC pulsations ! We don't check for the using_TDC flag, because mlt is sometimes called when using TDC - if ( s% TDC_num_innermost_cells_forced_nonturbulent > 0 .and. & - k > s% nz - s% TDC_num_innermost_cells_forced_nonturbulent) then + if (k > s% nz - s% TDC_num_innermost_cells_forced_nonturbulent .or. & + k < s% TDC_num_outermost_cells_forced_nonturbulent) then if (report) write(*,2) 'make TDC center cells non-turbulent', k mixing_type = no_mixing gradT = gradr diff --git a/star/public/star_lib.f90 b/star/public/star_lib.f90 index b3d5dac28..020538322 100644 --- a/star/public/star_lib.f90 +++ b/star/public/star_lib.f90 @@ -3119,16 +3119,17 @@ subroutine star_mlt_results(id, k, MLT_option, & ! NOTE: k=0 is a valid arg type(auto_diff_real_star_order1), intent(out) :: & gradT, Y_face, conv_vel, D, Gamma integer, intent(out) :: ierr - type(auto_diff_real_star_order1) :: dV + type(auto_diff_real_star_order1) :: dV, energy type (star_info), pointer :: s call star_ptr(id, s, ierr) if (ierr /= 0) return dV = 0d0 ! dV = 1/rho - 1/rho_start and we assume rho = rho_start. + energy = 0d0 ! for flux limiting TDC only call Get_results(s, k, MLT_option, & r, L, T, P, opacity, rho, dV, chiRho, chiT, Cp, gradr, grada, scale_height, & iso, XH1, cgrav, m, gradL_composition_term, mixing_length_alpha, & alpha_semiconvection, thermohaline_coeff, & - mixing_type, gradT, Y_face, conv_vel, D, Gamma, ierr) + mixing_type, gradT, Y_face, conv_vel, D, Gamma, energy, ierr) end subroutine star_mlt_results diff --git a/star/test_suite/12M_pre_ms_to_core_collapse/inlist_common b/star/test_suite/12M_pre_ms_to_core_collapse/inlist_common index cf0da6303..e49489993 100644 --- a/star/test_suite/12M_pre_ms_to_core_collapse/inlist_common +++ b/star/test_suite/12M_pre_ms_to_core_collapse/inlist_common @@ -99,7 +99,7 @@ ! If you would like to be bold, try radiative damping, although it ! is a departure from the mlt limit of 'TDC' - alpha_TDC_DAMPR = 0d0 ! 0d0 is default for mlt limit + TDC_alpha_R = 0d0 ! 0d0 is default for mlt limit use_other_alpha_mlt = .true. ! implemented in run_star_extras diff --git a/star/test_suite/20M_pre_ms_to_core_collapse/inlist_common b/star/test_suite/20M_pre_ms_to_core_collapse/inlist_common index f73ea365c..7bc607728 100644 --- a/star/test_suite/20M_pre_ms_to_core_collapse/inlist_common +++ b/star/test_suite/20M_pre_ms_to_core_collapse/inlist_common @@ -99,7 +99,7 @@ ! If you would like to be bold, try radiative damping, although it ! is a departure from the mlt limit of 'TDC' - alpha_TDC_DAMPR = 0d0 ! 0d0 is default for mlt limit + TDC_alpha_R = 0d0 ! 0d0 is default for mlt limit use_other_alpha_mlt = .true. ! implemented in run_star_extras diff --git a/star/test_suite/zams_to_cc_80/inlist_common b/star/test_suite/zams_to_cc_80/inlist_common index ed4ecd563..d2056e06f 100644 --- a/star/test_suite/zams_to_cc_80/inlist_common +++ b/star/test_suite/zams_to_cc_80/inlist_common @@ -99,7 +99,7 @@ ! If you would like to be bold, try radiative damping, although it ! is a departure from the mlt limit of 'TDC' - alpha_TDC_DAMPR = 0d0 ! 0d0 is default for mlt limit + TDC_alpha_R = 0d0 ! 0d0 is default for mlt limit use_other_alpha_mlt = .true. ! implemented in run_star_extras diff --git a/star_data/private/star_controls.inc b/star_data/private/star_controls.inc index d79ef7bfe..cff5b944e 100644 --- a/star_data/private/star_controls.inc +++ b/star_data/private/star_controls.inc @@ -255,9 +255,9 @@ logical :: use_Ledoux_criterion integer :: num_cells_for_smooth_gradL_composition_term integer :: steps_before_use_TDC - real(dp) :: alpha_TDC_DAMP ! ~ RSP2_alfad - real(dp) :: alpha_TDC_DAMPR ! ~ RSP2_alfar - real(dp) :: alpha_TDC_PtdVdt ! RSP2_alfap + real(dp) :: TDC_alpha_D ! ~ RSP2_alfad + real(dp) :: TDC_alpha_R ! ~ RSP2_alfar + real(dp) :: TDC_alpha_Pt ! RSP2_alfap real(dp) :: threshold_for_smooth_gradL_composition_term real(dp) :: clip_D_limit logical :: alt_scale_height_flag diff --git a/star_data/private/star_controls_dev.inc b/star_data/private/star_controls_dev.inc index 4a4a326c9..cffb14c1f 100644 --- a/star_data/private/star_controls_dev.inc +++ b/star_data/private/star_controls_dev.inc @@ -1,18 +1,21 @@ logical :: compare_TDC_to_MLT - real(dp) :: alpha_TDC_DAMPM - real(dp) :: alpha_TDC_C - real(dp) :: alpha_TDC_S + real(dp) :: TDC_alpha_M + real(dp) :: TDC_alpha_C + real(dp) :: TDC_alpha_S + logical :: TDC_alpha_M_use_explicit_mlt_vc_in_momentum_equation logical :: remesh_for_TDC_pulsations_log_core_zoning logical :: TDC_use_density_form_for_eddy_viscosity logical :: use_RSP_L_eqn_outer_BC real(dp) :: TDC_num_innermost_cells_forced_nonturbulent + real(dp) :: TDC_num_outermost_cells_forced_nonturbulent logical :: include_mlt_Pturb_in_thermodynamic_gradients - logical :: include_mlt_corr_to_TDC + logical :: include_mlt_corr_to_TDC, use_TDC_enthalpy_flux_limiter logical :: TDC_include_eturb_in_energy_equation, use_rsp_form_of_scale_height logical :: include_mlt_in_velocity_time_centering logical :: TDC_hydro_use_mass_interp_face_values + logical :: use_hydro_merge_limits_in_mesh_plan integer :: TDC_hydro_nz, TDC_hydro_nz_outer real(dp) :: TDC_hydro_T_anchor, TDC_hydro_dq_1_factor diff --git a/star_data/public/star_data_def.inc b/star_data/public/star_data_def.inc index 044d893a8..d265c52a6 100644 --- a/star_data/public/star_data_def.inc +++ b/star_data/public/star_data_def.inc @@ -799,14 +799,14 @@ r, L, T, P, opacity, rho, chiRho, chiT, Cp, gradr, grada, scale_height, & iso, XH1, cgrav, m, gradL_composition_term, mixing_length_alpha, & alpha_semiconvection, thermohaline_coeff, & - mixing_type, gradT, Y_face, conv_vel, D, Gamma, ierr) + mixing_type, gradT, Y_face, conv_vel, D, Gamma, energy, ierr) use const_def, only: dp use auto_diff integer, intent(in) :: id integer, intent(in) :: k character (len=*), intent(in) :: MLT_option type(auto_diff_real_star_order1), intent(in) :: & - r, L, T, P, opacity, rho, chiRho, chiT, Cp, gradr, grada, scale_height + r, L, T, P, opacity, rho, chiRho, chiT, Cp, gradr, grada, scale_height, energy integer, intent(in) :: iso real(dp), intent(in) :: & XH1, cgrav, m, gradL_composition_term, & diff --git a/turb/private/tdc_support.f90 b/turb/private/tdc_support.f90 index 161662135..99d02e5b7 100644 --- a/turb/private/tdc_support.f90 +++ b/turb/private/tdc_support.f90 @@ -44,9 +44,9 @@ module tdc_support !! !! @param report Write debug output if true, not if false. !! @param mixing_length_alpha Mixing length parameter - !! @param alpha_TDC_DAMP TDC turbulent damping parameter - !! @param alpha_TDC_DAMPR TDC radiative damping parameter - !! @param alpha_TDC_PtdVdt TDC coefficient on P_turb*dV/dt. Physically should probably be 1. + !! @param TDC_alpha_D TDC turbulent damping parameter + !! @param TDC_alpha_R TDC radiative damping parameter + !! @param TDC_alpha_Pt TDC coefficient on P_turb*dV/dt. Physically should probably be 1. !! @param dt Time-step !! @param c0 A proportionality factor for the convective luminosity !! @param L luminosity @@ -62,10 +62,10 @@ module tdc_support !! @param grada grada is the adiabatic dlnT/dlnP, !! @param Gamma Gamma is the MLT Gamma efficiency parameter, which we evaluate in steady state from MLT. type tdc_info - logical :: report, include_mlt_corr_to_TDC - real(dp) :: mixing_length_alpha, alpha_TDC_c, alpha_TDC_s, alpha_TDC_DAMP, alpha_TDC_DAMPR, alpha_TDC_PtdVdt, dt, e + logical :: report, include_mlt_corr_to_TDC, use_TDC_enthalpy_flux_limiter + real(dp) :: mixing_length_alpha, TDC_alpha_C, TDC_alpha_S, TDC_alpha_D, TDC_alpha_R, TDC_alpha_Pt, dt, e type(auto_diff_real_tdc) :: A0, c0, L, L0, gradL, grada - type(auto_diff_real_star_order1) :: T, rho, dV, Cp, kap, Hp, Gamma, Eq_div_w, P + type(auto_diff_real_star_order1) :: T, rho, dV, Cp, kap, Hp, Gamma, Eq_div_w, P, h end type tdc_info contains @@ -518,17 +518,43 @@ subroutine eval_xis(info, Y, xi0, xi1, xi2) type(auto_diff_real_tdc), intent(out) :: xi0, xi1, xi2 type(auto_diff_real_tdc) :: S0, D0, DR0 type(auto_diff_real_star_order1) :: gammar_div_alfa, Pt0, dVdt + type(auto_diff_real_tdc) :: X, FL, scale real(dp), parameter :: x_ALFAS = (1.d0/2.d0)*sqrt_2_div_3 real(dp), parameter :: x_CEDE = (8.d0/3.d0)*sqrt_2_div_3 real(dp), parameter :: x_ALFAP = 2.d0/3.d0 real(dp), parameter :: x_GAMMAR = 2.d0*sqrt(3.d0) - S0 = convert(info%alpha_TDC_s * x_ALFAS * info%mixing_length_alpha*info%Cp*info%T/info%Hp)*info%grada - 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) + ! flux-limiter scale factor from Wuchterl & Feuchtinger form + ! X ~ G/F, with G ∝ α α_s c_p Y_env and F = √(2/3)·w/T, w = E + P/rho + ! Here Y is already Y_env (from compute_Q), so: + ! X = α α_s x_ALFAS * Y_env / sqrt(2/3) + ! scale = 1 for small X (linear regime) + ! scale ~ FL(X)/X for large X (saturating regime) + + scale = 1d0 + if (Y > 0d0 .and. info%use_TDC_enthalpy_flux_limiter) then + ! X = G/F + X = convert(info%Cp*info%T/info%h)*info%mixing_length_alpha * info%TDC_alpha_S * x_ALFAS * Y / sqrt_2_div_3 + FL = flux_limiter_function(X) + ! Avoid 0/0 or tiny/tiny; for X ≈ 0, FL ≈ X so scale ~ 1 anyway. + if (abs(X%val) >= 0.95d0) then + scale = FL / X + else + scale = 1d0 + end if + end if + + S0 = convert(info%TDC_alpha_S * x_ALFAS * info%mixing_length_alpha*info%Cp*info%T/info%Hp)*info%grada + if (info%use_TDC_enthalpy_flux_limiter) then + S0 = S0 * (scale * Y) + convert(info%Eq_div_w) + else ! no flux limiting + S0 = S0*Y + convert(info%Eq_div_w) + end if + + D0 = convert(info%TDC_alpha_D*x_CEDE/(info%mixing_length_alpha*info%Hp)) + gammar_div_alfa = info%TDC_alpha_R*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)) - Pt0 = info%alpha_TDC_PtdVdt*x_ALFAP*info%rho + Pt0 = info%TDC_alpha_Pt*x_ALFAP*info%rho dVdt = info%dV/info%dt xi0 = S0 @@ -536,6 +562,42 @@ subroutine eval_xis(info, Y, xi0, xi1, xi2) xi2 = -D0 end subroutine eval_xis + type(auto_diff_real_tdc) function flux_limiter_function(X) result(FL) ! should be c2 continuous + type(auto_diff_real_tdc), intent(in) :: X + real(dp), parameter :: X0 = 0.95_dp ! start of transition + real(dp), parameter :: delta = 0.05_dp ! width of transition + real(dp), parameter :: X1 = 1d0 !X0 + delta ! end of transition + + type(auto_diff_real_tdc) :: s, p + + ! Region 1: purely linear, FL = X + if (X%val < X0) then ! should not be encountered + FL = X + + ! Region 3: saturated, FL = 1 + else if (X%val >= X1) then + FL = 1.0_dp + + ! Region 2: smooth C² transition between the two + else + ! Normalized coordinate in [0,1] + s = (X - X0) / (X1 - X0) + + ! Quintic "smootherstep" polynomial: + ! p(s) = 10 s^3 - 15 s^4 + 6 s^5 + ! p(0)=0, p(1)=1, p'(0)=p'(1)=0, p''(0)=p''(1)=0 + p = pow3(s) * (10.0_dp + s * (-15.0_dp + 6.0_dp * s)) + + ! Blend between line FL=X and flat FL=1 + ! At s=0: FL = X + ! At s=1: FL = 1 + ! Because p', p'' vanish at 0 and 1, FL, FL', FL'' all match. + FL = X + (1.0_dp - X) * p + end if + end function flux_limiter_function + + + !! Calculates the solution to the TDC velocity equation. !! The velocity equation is !! diff --git a/turb/public/turb.f90 b/turb/public/turb.f90 index 6b70fa2aa..8407bcb88 100644 --- a/turb/public/turb.f90 +++ b/turb/public/turb.f90 @@ -83,9 +83,9 @@ end subroutine set_thermohaline !! !! @param conv_vel_start The convection speed at the start of the step. !! @param mixing_length_alpha The mixing length parameter. - !! @param alpha_TDC_DAMP TDC turbulent damping parameter - !! @param alpha_TDC_DAMPR TDC radiative damping parameter - !! @param alpha_TDC_PtdVdt TDC coefficient on P_turb*dV/dt. Physically should probably be 1. + !! @param TDC_alpha_D TDC turbulent damping parameter + !! @param TDC_alpha_R TDC radiative damping parameter + !! @param TDC_alpha_Pt TDC coefficient on P_turb*dV/dt. Physically should probably be 1. !! @param The time-step (s). !! @param cgrav gravitational constant (erg*cm/g^2). !! @param m Mass inside the face (g). @@ -112,21 +112,22 @@ end subroutine set_thermohaline !! @param tdc_num_iters Number of iterations taken in the TDC solver. !! @param ierr Tracks errors (output). subroutine set_TDC( & - conv_vel_start, mixing_length_alpha, alpha_TDC_DAMP, alpha_TDC_DAMPR, alpha_TDC_PtdVdt, dt, cgrav, m, report, & + conv_vel_start, mixing_length_alpha, TDC_alpha_D, TDC_alpha_R, TDC_alpha_Pt, 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, max_conv_vel, & - Eq_div_w, grav, include_mlt_corr_to_TDC, alpha_TDC_C, alpha_TDC_S, ierr) + scale_height, gradL, grada, conv_vel, D, Y_face, gradT, tdc_num_iters, & + max_conv_vel, Eq_div_w, grav, include_mlt_corr_to_TDC, TDC_alpha_C, & + TDC_alpha_S, use_TDC_enthalpy_flux_limiter, energy, 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 - real(dp), intent(in) :: dt, cgrav, m, scale, max_conv_vel, alpha_TDC_c, alpha_TDC_s + real(dp), intent(in) :: conv_vel_start, mixing_length_alpha, TDC_alpha_D, TDC_alpha_R, TDC_alpha_Pt + real(dp), intent(in) :: dt, cgrav, m, scale, max_conv_vel, TDC_alpha_C, TDC_alpha_S type(auto_diff_real_star_order1), intent(in) :: & - chiT, chiRho, gradr, r, P, T, rho, dV, Cp, opacity, scale_height, gradL, grada, Eq_div_w, grav - logical, intent(in) :: report, include_mlt_corr_to_TDC + chiT, chiRho, gradr, r, P, T, rho, dV, Cp, opacity, scale_height, gradL, grada, Eq_div_w, grav, energy + logical, intent(in) :: report, include_mlt_corr_to_TDC, use_TDC_enthalpy_flux_limiter type(auto_diff_real_star_order1),intent(out) :: conv_vel, Y_face, gradT, D integer, intent(out) :: tdc_num_iters, mixing_type, ierr type(tdc_info) :: info - type(auto_diff_real_star_order1) :: L, Lambda, Gamma + type(auto_diff_real_star_order1) :: L, Lambda, Gamma, h real(dp), parameter :: alpha_c = (1d0/2d0)*sqrt_2_div_3 real(dp), parameter :: lower_bound_Z = -1d2 real(dp), parameter :: upper_bound_Z = 1d2 @@ -147,17 +148,19 @@ subroutine set_TDC( & ! Pack TDC info info%report = report info%include_mlt_corr_to_TDC = include_mlt_corr_to_TDC + info%use_TDC_enthalpy_flux_limiter = use_TDC_enthalpy_flux_limiter info%mixing_length_alpha = mixing_length_alpha - info%alpha_TDC_DAMP = alpha_TDC_DAMP - info%alpha_TDC_DAMPR = alpha_TDC_DAMPR - info%alpha_TDC_PtdVdt = alpha_TDC_PtdVdt + info%TDC_alpha_D = TDC_alpha_D + info%TDC_alpha_R = TDC_alpha_R + info%TDC_alpha_Pt = TDC_alpha_Pt info%dt = dt info%L = convert(L) info%gradL = convert(gradL) info%grada = convert(grada) - info%c0 = convert(alpha_TDC_C * mixing_length_alpha * alpha_c * rho * T * Cp * 4d0 * pi * pow2(r)) - info%L0 = convert((16d0*pi*crad*clight/3d0)*cgrav*m*pow4(T)/(P*opacity)) ! assumes QHSE for dP/dm, needs correction for if s% make_mlt_hydrodynamic = .false. + info%c0 = convert(TDC_alpha_C * mixing_length_alpha * alpha_c * rho * T * Cp * 4d0 * pi * pow2(r)) + info%L0 = convert((16d0*pi*crad*clight/3d0)*cgrav*m*pow4(T)/(P*opacity)) ! assumes QHSE for dP/dm info%A0 = conv_vel_start/sqrt_2_div_3 + info%h = energy + P/rho ! actual enthalpy info%T = T info%rho = rho info%dV = dV @@ -166,8 +169,8 @@ subroutine set_TDC( & info%Hp = scale_height info%Gamma = Gamma info%Eq_div_w = Eq_div_w - info%alpha_TDC_c = alpha_TDC_C - info%alpha_TDC_s = alpha_TDC_S + info%TDC_alpha_C = TDC_alpha_C + info%TDC_alpha_S = TDC_alpha_S ! Get solution Zub = upper_bound_Z diff --git a/turb/test/src/test_turb.f90 b/turb/test/src/test_turb.f90 index 8c4d91d85..2e71a34b3 100644 --- a/turb/test/src/test_turb.f90 +++ b/turb/test/src/test_turb.f90 @@ -81,15 +81,15 @@ end subroutine check_efficient_MLT_scaling 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, L_start, alpha_TDC_C, alpha_TDC_S + TDC_alpha_D, TDC_alpha_R, TDC_alpha_Pt, dt, cgrav, m, scale, L_start, TDC_alpha_C, TDC_alpha_S 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, Eq_div_w + type(auto_diff_real_star_order1) :: gradT, Y_face, conv_vel, D, Gamma, Eq_div_w, energy real(dp) :: Henyey_MLT_nu_param, Henyey_MLT_y_param, max_conv_vel character(len=3) :: MLT_option integer :: mixing_type, ierr, tdc_num_iters - logical :: report, include_mlt_corr_to_TDC + logical :: report, include_mlt_corr_to_TDC, use_TDC_enthalpy_flux_limiter include 'formats' @@ -124,18 +124,20 @@ subroutine compare_TDC_and_Cox_MLT() gradr = 3d0*P*opacity*L/(64*pi*boltz_sigma*pow4(T)*grav*pow2(r)) ! TDC - alpha_TDC_DAMP = 1.0d0 - alpha_TDC_DAMPR = 0.0d0 - alpha_TDC_PtdVdt = 0.0d0 - alpha_TDC_C = 1.0d0 - alpha_TDC_S = 1.0d0 + TDC_alpha_D = 1.0d0 + TDC_alpha_R = 0.0d0 + TDC_alpha_Pt = 0.0d0 + TDC_alpha_C = 1.0d0 + TDC_alpha_S = 1.0d0 dV = 0d0 + energy = 0d0 conv_vel_start = 0d0 !1d10 scale = L%val*1d-3 report = .false. dt = 1d40 ! Long time-step so we get into equilibrium - Eq_div_w = 0d0 + Eq_div_w = 0d0 ! TDC_alpha_M is implicit in this term include_mlt_corr_to_TDC = .true. + use_TDC_enthalpy_flux_limiter = .false. ! MLT MLT_option = 'Cox' @@ -145,10 +147,11 @@ subroutine compare_TDC_and_Cox_MLT() write (*, 1) 'gradR - gradA', gradr%val - grada%val call set_TDC( & - conv_vel_start, mixing_length_alpha, alpha_TDC_DAMP, alpha_TDC_DAMPR, alpha_TDC_PtdVdt, dt, cgrav, m, report, & + conv_vel_start, mixing_length_alpha, TDC_alpha_D, TDC_alpha_R, TDC_alpha_Pt, 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, max_conv_vel, & - Eq_div_w, grav, include_mlt_corr_to_TDC, alpha_TDC_C, alpha_TDC_S, ierr) + Eq_div_w, grav, include_mlt_corr_to_TDC, TDC_alpha_C, TDC_alpha_S, use_TDC_enthalpy_flux_limiter, & + energy, ierr) write (*, 1) 'TDC: Y, conv_vel_start, conv_vel, dt ', Y_face%val, conv_vel_start, conv_vel%val, dt @@ -164,13 +167,12 @@ end subroutine compare_TDC_and_Cox_MLT subroutine check_TDC() real(dp) :: mixing_length_alpha, conv_vel_start - real(dp) :: alpha_TDC_DAMP, alpha_TDC_DAMPR, alpha_TDC_PtdVdt, dt, cgrav, m, scale, max_conv_vel, L_start, & - alpha_TDC_C, alpha_TDC_S + 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 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, Eq_div_w, grav + type(auto_diff_real_star_order1) :: gradT, Y_face, conv_vel, D, Eq_div_w, grav, energy integer :: mixing_type, ierr, tdc_num_iters - logical :: report, include_mlt_corr_to_TDC + logical :: report, include_mlt_corr_to_TDC, use_TDC_enthalpy_flux_limiter integer :: j include 'formats' @@ -183,11 +185,11 @@ subroutine check_TDC() conv_vel_start = 52320587.415154047d0 mixing_length_alpha = 2.0d0 - alpha_TDC_DAMP = 1.0d0 - alpha_TDC_DAMPR = 0.0d0 - alpha_TDC_PtdVdt = 0.0d0 - alpha_TDC_C = 1.0d0 - alpha_TDC_S = 1.0d0 + TDC_alpha_D = 1.0d0 + TDC_alpha_R = 0.0d0 + TDC_alpha_Pt = 0.0d0 + TDC_alpha_C = 1.0d0 + TDC_alpha_S = 1.0d0 cgrav = 6.6743000000000004d-8 m = 5.8707400456875664d34 scale = 5.0386519362246294d45 @@ -208,8 +210,10 @@ subroutine check_TDC() gradr = 3d0 * P * opacity * L / (64 * pi * boltz_sigma * pow4(T) * cgrav * m) grav = m * cgrav / pow2(r) - Eq_div_w = 0d0 + Eq_div_w = 0d0 ! TDC_alpha_M is implicit in this term + energy = 0d0 include_mlt_corr_to_TDC = .true. + use_TDC_enthalpy_flux_limiter = .false. write (*, *) "####################################" write (*, *) "Running dt test" @@ -217,10 +221,11 @@ subroutine check_TDC() do j = 0, 30 dt = 500d0*pow(1.02d0, j) call set_TDC( & - conv_vel_start, mixing_length_alpha, alpha_TDC_DAMP, alpha_TDC_DAMPR, alpha_TDC_PtdVdt, dt, cgrav, m, report, & + conv_vel_start, mixing_length_alpha, TDC_alpha_D, TDC_alpha_R, TDC_alpha_Pt, 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, max_conv_vel, & - Eq_div_w, grav, include_mlt_corr_to_TDC, alpha_TDC_C, alpha_TDC_S, ierr) + Eq_div_w, grav, include_mlt_corr_to_TDC, TDC_alpha_C, TDC_alpha_S, use_TDC_enthalpy_flux_limiter, & + energy, ierr) write (*, 1) 'dt, gradT, conv_vel_start, conv_vel', dt, gradT%val, conv_vel_start, conv_vel%val