From 4facd2108bcbd64b83fb49f3e64ce021e553cf8a Mon Sep 17 00:00:00 2001 From: Loren Date: Fri, 23 Dec 2022 18:29:23 -0800 Subject: [PATCH 1/4] cleaned up the initialization block for PDE_Coefficients (no functional change) --- src/Physics/PDE_Coefficients.F90 | 100 ++++++++++++++++--------------- 1 file changed, 52 insertions(+), 48 deletions(-) diff --git a/src/Physics/PDE_Coefficients.F90 b/src/Physics/PDE_Coefficients.F90 index ed546c29..6acd01d1 100755 --- a/src/Physics/PDE_Coefficients.F90 +++ b/src/Physics/PDE_Coefficients.F90 @@ -39,10 +39,10 @@ Module PDE_Coefficients Use General_MPI, Only : BCAST2D Implicit None - Integer, Parameter :: n_scalar_max = 50 !/////////////////////////////////////////////////////////// ! I. Variables describing the background reference state + ! The following object (handle for which is "ref") is used by the code to access the PDE coefficients Type ReferenceInfo Real*8, Allocatable :: Density(:) Real*8, Allocatable :: dlnrho(:) @@ -57,20 +57,52 @@ Module PDE_Coefficients Real*8 :: Coriolis_Coeff ! Multiplies z_hat x u in momentum eq. Real*8 :: Lorentz_Coeff ! Multiplies (Del X B) X B in momentum eq. - Real*8, Allocatable :: Buoyancy_Coeff(:) ! Multiplies {S,T} in momentum eq. ..typically = gravity/cp - Real*8, Allocatable :: chi_buoyancy_coeff(:,:) ! Multiplies Chis in momentum eq. - Real*8, Allocatable :: dpdr_w_term(:) ! multiplies d_by_dr{P/rho} in momentum eq. - Real*8, Allocatable :: pressure_dwdr_term(:) !multiplies l(l+1)/r^2 (P/rho) in Div dot momentum eq. + Real*8, Allocatable :: Buoyancy_Coeff(:) ! Multiplies {S,T} in momentum eq. ..typically = gravity/cp + Real*8, Allocatable :: chi_buoyancy_coeff(:,:) ! Multiplies Chis in momentum eq. + Real*8, Allocatable :: dpdr_w_term(:) ! Multiplies d_by_dr{P/rho} in momentum eq. + Real*8, Allocatable :: pressure_dwdr_term(:) ! Multiplies l(l+1)/r^2 (P/rho) in Div dot momentum eq. ! The following two terms are used to compute the ohmic and viscous heating - Real*8, Allocatable :: ohmic_amp(:) !multiplied by {eta(r),H(r)}J^2 in dSdt eq. - Real*8, Allocatable :: viscous_amp(:) !multiplied by {nu(r),N(r)}{e_ij terms) in dSdt eq. + Real*8, Allocatable :: ohmic_amp(:) ! Multiplied by {eta(r),H(r)}J^2 in dSdt eq. + Real*8, Allocatable :: viscous_amp(:) ! Multiplied by {nu(r),N(r)}{e_ij terms) in dSdt eq. End Type ReferenceInfo + Type(ReferenceInfo) :: ref + ! Allow up to 50 active/passive scalar fields + Integer, Parameter :: n_scalar_max = 50 + + ! Version number for the "equation_coefficients" container that is output to the simulation directory + ! (i.e., the human-obtainable version of "ref") Integer, Parameter :: eqn_coeff_version = 1 - ! Custom reference state variables + ! Which background state to use; default 1 (non-dimensional Boussinesq) + Integer :: reference_type = 1 + + ! Nondimensional variables (reference_type = 1,3) + Real*8 :: Rayleigh_Number = 1.0d0 + Real*8 :: Ekman_Number = 1.0d0 + Real*8 :: Prandtl_Number = 1.0d0 + Real*8 :: Magnetic_Prandtl_Number = 1.0d0 + Real*8 :: gravity_power = 0.0d0 + Real*8 :: Dissipation_Number = 0.0d0 + Real*8 :: Modified_Rayleigh_Number = 0.0d0 + + ! Nondimensional variables for the active/passive scalar fields + Real*8 :: chi_a_rayleigh_number(1:n_scalar_max) = 0.0d0 + Real*8 :: chi_a_prandtl_number(1:n_scalar_max) = 1.0d0 + Real*8 :: chi_a_modified_rayleigh_number(1:n_scalar_max) = 0.0d0 + Real*8 :: chi_p_prandtl_number(1:n_scalar_max) = 1.0d0 + + ! Dimensional anelastic variables (reference_type = 2) + Real*8 :: pressure_specific_heat = 1.0d0 ! CP (not CV) + Real*8 :: poly_n = 0.0d0 ! Polytropic index + Real*8 :: poly_Nrho = 0.0d0 ! Number of density scale heights across domain + Real*8 :: poly_mass = 0.0d0 ! Stellar mass; g(r) = G*poly_mass/r^2 + Real*8 :: poly_rho_i = 0.0d0 ! Density (g/cm^3) at the inner radius rmin + Real*8 :: Angular_Velocity = -1.0d0 ! Frame rotation rate (sets Coriolis force) + + ! Custom reference-state variables (reference_type = 4) Integer, Parameter :: max_ra_constants = 10 + 2*n_scalar_max Integer, Parameter :: max_ra_functions = 14 + 2*n_scalar_max Integer :: n_ra_constants @@ -89,57 +121,29 @@ Module PDE_Coefficients Logical :: custom_reference_read = .false. Character*120 :: custom_reference_file ='nothing' - Real*8, Allocatable :: s_conductive(:) - - Integer :: reference_type =1 + ! Internal heating variables Integer :: heating_type = 0 ! 0 means no reference heating. > 0 selects optional reference heating - Real*8 :: Luminosity = 0.0d0 ! specifies the integral of the heating function - Real*8 :: Heating_Integral = 0.0d0 !same as luminosity (for non-star watchers) - Real*8 :: Heating_EPS = 1.0D-12 !Small number to test whether luminosity specified - Logical :: adjust_reference_heating = .false. ! Flag used to decide if luminosity determined via boundary conditions - - Type(ReferenceInfo) :: ref - - Real*8 :: pressure_specific_heat = 1.0d0 ! CP (not CV) - Real*8 :: poly_n = 0.0d0 !polytropic index - Real*8 :: poly_Nrho = 0.0d0 - Real*8 :: poly_mass = 0.0d0 - Real*8 :: poly_rho_i =0.0d0 - Real*8 :: Angular_Velocity = -1.0d0 - - !///////////////////////////////////////////////////////////////////////////////////// - ! Nondimensional Parameters - Real*8 :: Rayleigh_Number = 1.0d0 - Real*8 :: Ekman_Number = 1.0d0 - Real*8 :: Prandtl_Number = 1.0d0 - Real*8 :: Magnetic_Prandtl_Number = 1.0d0 - Real*8 :: gravity_power = 0.0d0 - Real*8 :: Dissipation_Number = 0.0d0 - Real*8 :: Modified_Rayleigh_Number = 0.0d0 - - Real*8 :: chi_a_rayleigh_number(1:n_scalar_max) = 0.0d0 - Real*8 :: chi_a_prandtl_number(1:n_scalar_max) = 1.0d0 - Real*8 :: chi_a_modified_rayleigh_number(1:n_scalar_max) = 0.0d0 - Real*8 :: chi_p_prandtl_number(1:n_scalar_max) = 1.0d0 + Real*8 :: Luminosity = 0.0d0 ! Specifies the integral of the heating function + Real*8 :: Heating_Integral = 0.0d0 ! Same as luminosity (for non-star watchers) + Real*8 :: Heating_EPS = 1.0d-12 ! Small number to test whether luminosity was specified + Logical :: adjust_reference_heating = .false. ! Flag used to decide if luminosity determined via boundary conditions + Real*8, Allocatable :: s_conductive(:) - !/////////////////////////////////////////////// ! Minimum time step based on rotation rate ! (determined by the reference state) Real*8 :: max_dt_rotation = 0.0d0 - - - Namelist /Reference_Namelist/ reference_type,poly_n, poly_Nrho, poly_mass,poly_rho_i, & - & pressure_specific_heat, heating_type, luminosity, Angular_Velocity, & + ! Alter some of the above (I) by reading the main_input file + Namelist /Reference_Namelist/ reference_type,poly_n, poly_Nrho, poly_mass, poly_rho_i, & + & pressure_specific_heat, heating_type, luminosity, Angular_Velocity, & & Rayleigh_Number, Ekman_Number, Prandtl_Number, Magnetic_Prandtl_Number, & - & gravity_power, custom_reference_file, & + & gravity_power, custom_reference_file, & & Dissipation_Number, Modified_Rayleigh_Number, & & Heating_Integral, override_constants, override_constant, ra_constants, with_custom_constants, & & with_custom_functions, with_custom_reference, & & chi_a_rayleigh_number, chi_a_prandtl_number, & & chi_a_modified_rayleigh_number, chi_p_prandtl_number - !/////////////////////////////////////////////////////////////////////////////////////// ! II. Variables Related to the Transport Coefficients @@ -157,7 +161,7 @@ Module PDE_Coefficients Real*8, Allocatable :: A_Diffusion_Coefs_1(:) real*8, allocatable :: chi_a_diffusion_coefs_1(:,:), chi_p_diffusion_coefs_1(:,:) - Integer :: kappa_type =1, nu_type = 1, eta_type = 1 + Integer :: kappa_type = 1, nu_type = 1, eta_type = 1 Real*8 :: nu_top = -1.0d0, kappa_top = -1.0d0, eta_top = -1.0d0 Real*8 :: nu_power = 0, eta_power = 0, kappa_power = 0 Integer :: kappa_chi_a_type(1:n_scalar_max) = 1 @@ -171,6 +175,7 @@ Module PDE_Coefficients Real*8 :: hyperdiffusion_beta = 0.0d0 Real*8 :: hyperdiffusion_alpha = 1.0d0 + ! Alter some of the above (II) by reading the main_input file Namelist /Transport_Namelist/ nu_type, kappa_type, eta_type, & & nu_power, kappa_power, eta_power, & & nu_top, kappa_top, eta_top, & @@ -178,7 +183,6 @@ Module PDE_Coefficients & kappa_chi_a_type, kappa_chi_a_top, kappa_chi_a_power, & & kappa_chi_p_type, kappa_chi_p_top, kappa_chi_p_power - Contains !///////////////////////////////////////////////////////////////// From 9422e3333623d171477e517a5de0a13d5f9e7f23 Mon Sep 17 00:00:00 2001 From: Loren Date: Fri, 23 Dec 2022 18:36:18 -0800 Subject: [PATCH 2/4] fixed errors in subroutine Restore_Reference_Defaults() --- src/Physics/PDE_Coefficients.F90 | 128 +++++++++++++++++++++++++++---- 1 file changed, 113 insertions(+), 15 deletions(-) diff --git a/src/Physics/PDE_Coefficients.F90 b/src/Physics/PDE_Coefficients.F90 index 6acd01d1..8eaf8ebb 100755 --- a/src/Physics/PDE_Coefficients.F90 +++ b/src/Physics/PDE_Coefficients.F90 @@ -1142,38 +1142,136 @@ End Subroutine log_deriv Subroutine Restore_Reference_Defaults Implicit None - !Restore all values in this module to their default state. - !Deallocates all allocatable module variables. - reference_type =1 - heating_type = 0 - Luminosity =0.0d0 + ! Restore all variables in this module to their default state. + ! Deallocates all allocatable module variables. - pressure_specific_heat = 0.0d0 ! CP (not CV) - poly_n = 0.0d0 - poly_Nrho = 0.0d0 - poly_mass = 0.0d0 - poly_rho_i = 0.0d0 - - Angular_Velocity = -1.0d0 + ! Which background state to use; default 1 (non-dimensional Boussinesq) + reference_type = 1 + ! Nondimensional variables (reference_type = 1,3) Rayleigh_Number = 1.0d0 Ekman_Number = 1.0d0 Prandtl_Number = 1.0d0 Magnetic_Prandtl_Number = 1.0d0 gravity_power = 0.0d0 + Dissipation_Number = 0.0d0 + Modified_Rayleigh_Number = 0.0d0 + + ! Nondimensional variables for the active/passive scalar fields + chi_a_rayleigh_number(1:n_scalar_max) = 0.0d0 + chi_a_prandtl_number(1:n_scalar_max) = 1.0d0 + chi_a_modified_rayleigh_number(1:n_scalar_max) = 0.0d0 + chi_p_prandtl_number(1:n_scalar_max) = 1.0d0 + + ! Dimensional anelastic variables (reference_type = 2) + pressure_specific_heat = 1.0d0 ! CP (not CV) + poly_n = 0.0d0 ! Polytropic index + poly_Nrho = 0.0d0 ! Number of density scale heights across domain + poly_mass = 0.0d0 ! Stellar mass; g(r) = G*poly_mass/r^2 + poly_rho_i = 0.0d0 ! Density (g/cm^3) at the inner radius rmin + Angular_Velocity = -1.0d0 ! Frame rotation rate (sets Coriolis force) + + ! Custom reference-state variables (reference_type = 4) + ! NOTE: n_ra_constants / n_ra_functions do not have default values (but maybe do, if you consider the Allocate_Reference_State() routine) + with_custom_reference = .false. + override_constants = .false. + override_constant(1:max_ra_constants) = .false. ! in namelist + with_custom_constants(1:max_ra_constants) = 0 ! in namelist + with_custom_functions(1:max_ra_functions) = 0 ! in namelist + ra_constants(1:max_ra_constants) = 0.0d0 ! in namelist + custom_reference_read = .false. + custom_reference_file ='nothing' + + ! Internal heating variables + heating_type = 0 ! 0 means no reference heating. > 0 selects optional reference heating + Luminosity = 0.0d0 ! Specifies the integral of the heating function + Heating_Integral = 0.0d0 ! Same as luminosity (for non-star watchers) + Heating_EPS = 1.0d-12 ! Small number to test whether luminosity was specified + adjust_reference_heating = .false. ! Flag used to decide if luminosity determined via boundary conditions + + ! Minimum time step based on rotation rate + ! (determined by the reference state) + max_dt_rotation = 0.0d0 + + ! Diffusion-coefficient (transport-namelist) variables + kappa_type = 1 + nu_type = 1 + eta_type = 1 + nu_top = -1.0d0 + kappa_top = -1.0d0 + eta_top = -1.0d0 + nu_power = 0 + eta_power = 0 + kappa_power = 0 - custom_reference_file ='nothing' + kappa_chi_a_type(1:n_scalar_max) = 1 + kappa_chi_a_top(1:n_scalar_max) = -1.0d0 + kappa_chi_a_power(1:n_scalar_max) = 0 + kappa_chi_p_type(1:n_scalar_max) = 1 + kappa_chi_p_top(1:n_scalar_max) = -1.0d0 + kappa_chi_p_power(1:n_scalar_max) = 0 - If (allocated(s_conductive)) DeAllocate(s_conductive) + hyperdiffusion = .false. + hyperdiffusion_beta = 0.0d0 + hyperdiffusion_alpha = 1.0d0 + + ! Deallocate "ref" object If (allocated(ref%Density)) DeAllocate(ref%density) If (allocated(ref%dlnrho)) DeAllocate(ref%dlnrho) If (allocated(ref%d2lnrho)) DeAllocate(ref%d2lnrho) + If (allocated(ref%Temperature)) DeAllocate(ref%Temperature) If (allocated(ref%dlnT)) DeAllocate(ref%dlnT) + If (allocated(ref%dsdr)) DeAllocate(ref%dsdr) + + If (allocated(ref%heating)) DeAllocate(ref%heating) + + ! NOTE: ref%Coriolis_Coeff and ref%Lorentz_Coeff have no default values (but maybe do, if you consider the Allocate_Reference_State() routine) If (allocated(ref%Buoyancy_Coeff)) DeAllocate(ref%Buoyancy_Coeff) If (allocated(ref%chi_buoyancy_coeff)) DeAllocate(ref%chi_buoyancy_coeff) - If (allocated(ref%Heating)) DeAllocate(ref%Heating) + If (allocated(ref%dpdr_w_term)) DeAllocate(ref%dpdr_w_term) + If (allocated(ref%pressure_dwdr_term)) DeAllocate(ref%pressure_dwdr_term) + + If (allocated(ref%ohmic_amp)) DeAllocate(ref%ohmic_amp) + If (allocated(ref%viscous_amp)) DeAllocate(ref%viscous_amp) + + ! Deallocate s_conductive + If (allocated(s_conductive)) DeAllocate(s_conductive) + + ! Deallocate custom-reference stuff + If (allocated(ra_constant_set)) DeAllocate(ra_constant_set) + If (allocated(ra_function_set)) DeAllocate(ra_function_set) + If (allocated(use_custom_constant)) DeAllocate(use_custom_constant) + If (allocated(use_custom_function)) DeAllocate(use_custom_function) + If (allocated(ra_functions)) DeAllocate(ra_functions) + + ! Deallocate transport-coefficient stuff + If (allocated(nu)) DeAllocate(nu) + If (allocated(kappa)) DeAllocate(kappa) + If (allocated(eta)) DeAllocate(eta) + If (allocated(dlnu)) DeAllocate(dlnu) + If (allocated(dlnkappa)) DeAllocate(dlnkappa) + If (allocated(dlneta)) DeAllocate(dlneta) + If (allocated(kappa_chi_a)) DeAllocate(kappa_chi_a) + If (allocated(kappa_chi_p)) DeAllocate(kappa_chi_p) + If (allocated(dlnkappa_chi_a)) DeAllocate(dlnkappa_chi_a) + If (allocated(dlnkappa_chi_p)) DeAllocate(dlnkappa_chi_p) + + If (allocated(ohmic_heating_coeff)) DeAllocate(ohmic_heating_coeff) + If (allocated(viscous_heating_coeff)) DeAllocate(viscous_heating_coeff) + + If (allocated(W_Diffusion_Coefs_0)) DeAllocate(W_Diffusion_Coefs_0) + If (allocated(W_Diffusion_Coefs_1)) DeAllocate(W_Diffusion_Coefs_1) + If (allocated(dW_Diffusion_Coefs_0)) DeAllocate(dW_Diffusion_Coefs_0) + If (allocated(dW_Diffusion_Coefs_1)) DeAllocate(dW_Diffusion_Coefs_1) + If (allocated(dW_Diffusion_Coefs_2)) DeAllocate(dW_Diffusion_Coefs_2) + If (allocated(S_Diffusion_Coefs_1)) DeAllocate(S_Diffusion_Coefs_1) + If (allocated(Z_Diffusion_Coefs_0)) DeAllocate(Z_Diffusion_Coefs_0) + If (allocated(Z_Diffusion_Coefs_1)) DeAllocate(Z_Diffusion_Coefs_1) + If (allocated(A_Diffusion_Coefs_1)) DeAllocate(A_Diffusion_Coefs_1) + If (allocated(chi_a_diffusion_coefs_1)) DeAllocate(chi_a_diffusion_coefs_1) + If (allocated(chi_p_diffusion_coefs_1)) DeAllocate(chi_p_diffusion_coefs_1) End Subroutine Restore_Reference_Defaults From 1625f627928726b74ea76c373a50c887f6c2ad3e Mon Sep 17 00:00:00 2001 From: Loren Date: Tue, 3 Jan 2023 16:05:26 -0800 Subject: [PATCH 3/4] tried removing new 'custom' stuff in Restore_Reference_Defaults() --- src/Physics/PDE_Coefficients.F90 | 24 ++++++++++++------------ 1 file changed, 12 insertions(+), 12 deletions(-) diff --git a/src/Physics/PDE_Coefficients.F90 b/src/Physics/PDE_Coefficients.F90 index 8eaf8ebb..570d2087 100755 --- a/src/Physics/PDE_Coefficients.F90 +++ b/src/Physics/PDE_Coefficients.F90 @@ -1173,13 +1173,13 @@ Subroutine Restore_Reference_Defaults ! Custom reference-state variables (reference_type = 4) ! NOTE: n_ra_constants / n_ra_functions do not have default values (but maybe do, if you consider the Allocate_Reference_State() routine) - with_custom_reference = .false. - override_constants = .false. - override_constant(1:max_ra_constants) = .false. ! in namelist - with_custom_constants(1:max_ra_constants) = 0 ! in namelist - with_custom_functions(1:max_ra_functions) = 0 ! in namelist - ra_constants(1:max_ra_constants) = 0.0d0 ! in namelist - custom_reference_read = .false. + !with_custom_reference = .false. + !override_constants = .false. + !override_constant(1:max_ra_constants) = .false. ! in namelist + !with_custom_constants(1:max_ra_constants) = 0 ! in namelist + !with_custom_functions(1:max_ra_functions) = 0 ! in namelist + !ra_constants(1:max_ra_constants) = 0.0d0 ! in namelist + !custom_reference_read = .false. custom_reference_file ='nothing' ! Internal heating variables @@ -1240,11 +1240,11 @@ Subroutine Restore_Reference_Defaults If (allocated(s_conductive)) DeAllocate(s_conductive) ! Deallocate custom-reference stuff - If (allocated(ra_constant_set)) DeAllocate(ra_constant_set) - If (allocated(ra_function_set)) DeAllocate(ra_function_set) - If (allocated(use_custom_constant)) DeAllocate(use_custom_constant) - If (allocated(use_custom_function)) DeAllocate(use_custom_function) - If (allocated(ra_functions)) DeAllocate(ra_functions) + !If (allocated(ra_constant_set)) DeAllocate(ra_constant_set) + !If (allocated(ra_function_set)) DeAllocate(ra_function_set) + !If (allocated(use_custom_constant)) DeAllocate(use_custom_constant) + !If (allocated(use_custom_function)) DeAllocate(use_custom_function) + !If (allocated(ra_functions)) DeAllocate(ra_functions) ! Deallocate transport-coefficient stuff If (allocated(nu)) DeAllocate(nu) From b6c514510d138415d763728ce0c4f5eaac76a5ee Mon Sep 17 00:00:00 2001 From: Loren Matilsky Date: Thu, 5 Jan 2023 12:32:08 -0800 Subject: [PATCH 4/4] Remove commented-out reset of custom-reference variables ...in Restore_Reference_Defaults(). As per Nick's request. --- src/Physics/PDE_Coefficients.F90 | 14 -------------- 1 file changed, 14 deletions(-) diff --git a/src/Physics/PDE_Coefficients.F90 b/src/Physics/PDE_Coefficients.F90 index 570d2087..cdcd2dba 100755 --- a/src/Physics/PDE_Coefficients.F90 +++ b/src/Physics/PDE_Coefficients.F90 @@ -1173,13 +1173,6 @@ Subroutine Restore_Reference_Defaults ! Custom reference-state variables (reference_type = 4) ! NOTE: n_ra_constants / n_ra_functions do not have default values (but maybe do, if you consider the Allocate_Reference_State() routine) - !with_custom_reference = .false. - !override_constants = .false. - !override_constant(1:max_ra_constants) = .false. ! in namelist - !with_custom_constants(1:max_ra_constants) = 0 ! in namelist - !with_custom_functions(1:max_ra_functions) = 0 ! in namelist - !ra_constants(1:max_ra_constants) = 0.0d0 ! in namelist - !custom_reference_read = .false. custom_reference_file ='nothing' ! Internal heating variables @@ -1238,13 +1231,6 @@ Subroutine Restore_Reference_Defaults ! Deallocate s_conductive If (allocated(s_conductive)) DeAllocate(s_conductive) - - ! Deallocate custom-reference stuff - !If (allocated(ra_constant_set)) DeAllocate(ra_constant_set) - !If (allocated(ra_function_set)) DeAllocate(ra_function_set) - !If (allocated(use_custom_constant)) DeAllocate(use_custom_constant) - !If (allocated(use_custom_function)) DeAllocate(use_custom_function) - !If (allocated(ra_functions)) DeAllocate(ra_functions) ! Deallocate transport-coefficient stuff If (allocated(nu)) DeAllocate(nu)