diff --git a/include/aspect/boundary_temperature/dynamic_core.h b/include/aspect/boundary_temperature/dynamic_core.h index 1f7752a38c1..e9d4dc25a54 100644 --- a/include/aspect/boundary_temperature/dynamic_core.h +++ b/include/aspect/boundary_temperature/dynamic_core.h @@ -103,6 +103,27 @@ namespace aspect */ DynamicCore(); + /** + * This function update the core-mantle boundary (CMB) temperature by + * the core energy balance solver using the core-mantle boundary heat flux. + */ + void + update() override; + + /** + * Pass core data to other modules + */ + const internal::CoreData & + get_core_data() const; + + /** + * Check if other energy source in the core is in use. The 'other energy source' is used for external core energy source. + * For example if someone want to test the early lunar core powered by precession + * (Dwyer, C. A., et al. (2011). "A long-lived lunar dynamo driven by continuous mechanical stirring." Nature 479(7372): 212-214.) + */ + bool + is_OES_used() const; + /** * Return the temperature that is to hold at a particular location on the * boundary of the domain. This function returns the temperatures @@ -113,8 +134,9 @@ namespace aspect * temperature. * @param location The location of the point at which we ask for the temperature. */ - double boundary_temperature (const types::boundary_id boundary_indicator, - const Point &location) const override; + double + boundary_temperature (const types::boundary_id boundary_indicator, + const Point &location) const override; /** * Return the minimal temperature on that part of the boundary @@ -123,7 +145,8 @@ namespace aspect * This value is used in computing dimensionless numbers such as the * Nusselt number indicating heat flux. */ - double minimal_temperature (const std::set &fixed_boundary_ids) const override; + double + minimal_temperature (const std::set &fixed_boundary_ids) const override; /** * Return the maximal temperature on that part of the boundary @@ -132,7 +155,8 @@ namespace aspect * This value is used in computing dimensionless numbers such as the * Nusselt number indicating heat flux. */ - double maximal_temperature (const std::set &fixed_boundary_ids) const override; + double + maximal_temperature (const std::set &fixed_boundary_ids) const override; /** * Declare the parameters this class takes through input files. @@ -149,27 +173,6 @@ namespace aspect void parse_parameters (ParameterHandler &prm) override; - /** - * This function update the core-mantle boundary (CMB) temperature by - * the core energy balance solver using the core-mantle boundary heat flux. - */ - void update() override; - - /** - * Pass core data to other modules - */ - const internal::CoreData & - get_core_data() const; - - /** - * Check if other energy source in the core is in use. The 'other energy source' is used for external core energy source. - * For example if someone want to test the early lunar core powered by precession - * (Dwyer, C. A., et al. (2011). "A long-lived lunar dynamo driven by continuous mechanical stirring." Nature 479(7372): 212-214.) - */ - bool - is_OES_used() const; - - private: /** @@ -361,14 +364,12 @@ namespace aspect void read_data_OES(); double get_OES(double t) const; - - /** * Solve core energy balance for each time step. - * When solving the change in core-mantle boundary temperature T, inner core radius R, and - * light component (e.g. S, O, Si) composition X, the following relations has to be respected: + * When solving the change in core-mantle boundary temperature @p T, inner core radius @p R, and + * light element (e.g. S, O, Si) composition @p X, the following relations have to be respected: * 1. At the inner core boundary the adiabatic temperature should be equal to solidus temperature - * 2. The following energy production rate should be balanced in core: + * 2. The following energy production rates should be balanced in the core: * Heat flux at core-mantle boundary Q * Specific heat Qs*dT/dt * Radioactive heating Qr @@ -378,134 +379,137 @@ namespace aspect * 3. The light component composition X depends on inner core radius (See function get_X() ), * and core solidus may dependent on X as well. * This becomes a small nonlinear problem. Directly iterate through the above three equations doesn't - * converge well. Alternatively we solve the inner core radius by bisection method. - * A single solution between fully liquid and fully solid core is expected. Otherwise this function will throw exception and terminate. + * converge well. Alternatively we solve the inner core radius using the bisection method. * - * At Earth core condition, a inner core is forming at the center of the Earth and surrounded by a liquid outer core. - * However, the core solidus is influenced by light components (e.g. S) and its slope is very closed to core adiabatic. So there is an alternative + * At the conditions of the Earth's core, an inner core is forming at the center of the Earth and surrounded by a liquid outer core. + * However, the core solidus is influenced by light components (e.g. S) and its slope is very close to an adiabat. So there is an alternative * scenario that the crystallization happens first at the core mantle boundary instead of at the center, which is called a 'snowing core' * (Stewart, A. J., et al. (2007). "Mars: a new core-crystallization regime." Science 316(5829): 1323-1325.). This also - * provides a valid solution for the solver. So the returning bool is set to true for normal core, and false for 'snowing core'. - * TODO: The current code is only able to treat normal core scenario, treating 'snowing core' scenario may be possible and could be added. + * provides a valid solution for the solver. The return value of the function is true for a 'normal core', and false for 'snowing core'. + * TODO: The current code is only able to treat the normal core scenario, treating 'snowing core' scenario may be possible and could be added. */ - bool solve_time_step(double &X, double &T, double &R); + bool solve_time_step(double &X, double &T, double &R) const; /** - * Compute the difference between solidus and adiabatic at inner - * core boundary for a given inner core radius. + * Compute the difference between solidus and adiabatic temperature at inner + * core boundary for a given inner core radius @p r. */ - double get_dT(double r) const; + double get_dT(const double r) const; /** * Use energy balance to calculate core mantle boundary temperature - * with a given inner core radius. + * with a given inner core radius @p r. */ - double get_Tc(double r) const; + double get_Tc(const double r) const; /** * Get the solidus temperature at inner core boundary - * with a given inner core radius. + * with a given inner core radius @p r. */ - double get_Ts(double r) const; + double get_Ts(const double r) const; /** - * Compute the core solidus at certain pressure + * Compute the core solidus at certain pressure @p pressure + * and light element concentration @p X (in wt.%). */ - double get_solidus(double X,double p) const; + double get_solidus(const double X, const double pressure) const; /** - * Get initial inner core radius with given initial core mantle temperature. + * Get initial inner core radius with given initial core mantle temperature + * @p T. */ - double get_initial_Ri(double T); + double get_initial_Ri(const double T) const; /** - * Get the light composition concentration in the outer core from given - * inner core radius r + * Get the light element concentration (in wt.%) in the outer core from given + * inner core radius @p r. */ - double get_X(double r) const; + double get_X(const double r) const; /** - * Compute the mass inside certain radius within the core_data. + * Compute the core mass inside a certain radius @p r. */ - double get_Mass(double r) const; + double get_mass(const double r) const; /** * Calculate Sn(B,R), referring to \cite NPB+04 . */ - double fun_Sn(double B,double R,double n) const; + double fun_Sn(const double B, const double R, const double n) const; /** - * Calculate density at given r + * Calculate density at given radius @p r. */ - double get_Rho(double r) const; + double get_rho(const double r) const; /** - * Calculate gravitational acceleration at given r + * Calculate gravitational acceleration at given radius @p r. */ - double get_g(double r) const; + double get_g(const double r) const; /** - * Calculate the core temperature at given r - * Tc is the temperature at CMB + * Calculate the core temperature at given radius @p r and + * temperature at CMB @p Tc. */ - double get_T(double Tc, double r) const; + double get_T(const double Tc, const double r) const; /** - * Calculate pressure at given r + * Calculate pressure at given radius @p r */ - double get_Pressure(double r) const; + double get_pressure(const double r) const; /** - * Calculate the gravitational potential at given r + * Calculate the gravitational potential at given radius @p r */ - double get_gravity_potential(double r) const; + double get_gravity_potential(const double r) const; /** - * Calculate energy and entropy change rate factor (regarding the core - * cooling rated Tc/dt) Qs and Es with given core-mantle boundary (CMB) - * temperature Tc + * Calculate energy (@p Qs) and entropy (@p Es) change rate factor + * (regarding the core cooling rated Tc/dt) for a given core-mantle boundary (CMB) + * temperature @p Tc */ - void get_specific_heating(double Tc, double &Qs,double &Es); + void get_specific_heating(const double Tc, double &Qs, double &Es) const; /** - * Calculate energy and entropy change rate factor (regarding the - * radioactive heating rate H) Qr and Er with given CMB temperature Tc + * Calculate energy (@p Qr) and entropy (@p Er) change rate factor (regarding the + * radioactive heating rate H) for a given CMB temperature @p Tc */ - void get_radio_heating(double Tc, double &Qr, double &Er); + void get_radio_heating(const double Tc, double &Qr, double &Er) const; /** - * Calculate energy and entropy change rate factor (regarding the inner core - * growth rate dR/dt) Qg and Eg with given Tc(CMB temperature), r(inner core - * radius), X(light composition concentration) + * Calculate energy (@p Qg) and entropy (@p Eg) change rate factor + * (regarding the inner core growth rate dR/dt) for a given + * @p Tc (CMB temperature), @p r (inner core radius), and @p X + * (light element concentration) */ - void get_gravity_heating(double Tc, double r,double X,double &Qg,double &Eg); + void get_gravity_heating(const double Tc, const double r, const double X, double &Qg, double &Eg) const; /** - * Calculate energy and entropy change rate factor (regarding the core - * cooling rate Tc/dt) Qk and Ek with given Tc(CMB temperature) + * Calculate energy (@p Qk) and entropy (@p Ek) change rate factor + * (regarding the core cooling rate Tc/dt) for a given @p Tc (CMB temperature) */ - void get_adiabatic_heating(double Tc, double &Ek, double &Qk); + void get_adiabatic_heating(const double Tc, double &Ek, double &Qk) const; /** - * Calculate energy and entropy change rate factor (regarding the inner core - * growth rate dR/dt) Ql and El with given Tc(CMB temperature), r(inner core - * radius) + * Calculate energy (@p Ql) and entropy (@p El) change rate factor + * (regarding the inner core growth rate dR/dt) for a given @p Tc (CMB temperature) + * and @p r (inner core radius) */ - void get_latent_heating(double Tc, double r, double &El, double &Ql); + void get_latent_heating(const double Tc, const double r, double &El, double &Ql) const; /** - * Calculate entropy of heat of solution Eh + * Calculate entropy of heat of solution @p Eh for a given @p Tc (CMB temperature), + * @p r (inner core radius), and @p X (light element concentration) */ - void get_heat_solution(double Tc, double r, double X, double &Eh); + void get_heat_solution(const double Tc, const double r, const double X, double &Eh) const; /** - * return radio heating rate at certain time + * return radiogenic heating rate at the current time */ double get_radioheating_rate() const; /** - * Update the data for core dynamic simulation, the data will be used - * in the next timestep and for postprocess. + * Update the data of the core dynamic simulation, the data will be used + * in the next timestep and for postprocessing. */ void update_core_data(); diff --git a/source/boundary_temperature/dynamic_core.cc b/source/boundary_temperature/dynamic_core.cc index a5099c1427c..8c50696fca3 100644 --- a/source/boundary_temperature/dynamic_core.cc +++ b/source/boundary_temperature/dynamic_core.cc @@ -23,14 +23,14 @@ #include #include #include -#include -#include +#include +#include +#include #include -#include -#include +#include +#include -#include #include @@ -42,8 +42,8 @@ namespace aspect template double DynamicCore:: - boundary_temperature (const types::boundary_id boundary_indicator, - const Point &/*location*/) const + boundary_temperature (const types::boundary_id boundary_indicator, + const Point &/*location*/) const { switch (boundary_indicator) { @@ -58,6 +58,7 @@ namespace aspect } + template double DynamicCore:: @@ -218,6 +219,7 @@ namespace aspect } + template void DynamicCore::parse_parameters (ParameterHandler &prm) @@ -298,8 +300,8 @@ namespace aspect } prm.leave_subsection (); - L=std::sqrt(3*K0*(std::log(Rho_cen/Rho_0)+1)/(2*M_PI*constants::big_g*Rho_0*Rho_cen)); - D=std::sqrt(3*Cp/(2*M_PI*Alpha*Rho_cen*constants::big_g)); + L=std::sqrt(3*K0*(std::log(Rho_cen/Rho_0)+1)/(2*numbers::PI*constants::big_g*Rho_0*Rho_cen)); + D=std::sqrt(3*Cp/(2*numbers::PI*Alpha*Rho_cen*constants::big_g)); } prm.leave_subsection (); @@ -307,6 +309,8 @@ namespace aspect prm.leave_subsection (); } + + template void DynamicCore::read_data_OES() @@ -334,12 +338,14 @@ namespace aspect << std::endl; } + + template double - DynamicCore::get_OES(double t) const + DynamicCore::get_OES(const double time) const { // The core evolution is quite slow, so the time units used here is billion years. - t/=1.e9*year_in_seconds; + const double t = time / (1.e9*year_in_seconds); double w=0.; for (unsigned i=1; i DynamicCore::DynamicCore() { @@ -361,33 +369,34 @@ namespace aspect core_data.is_initialized = false; } + + template double - DynamicCore::get_initial_Ri(double T) + DynamicCore::get_initial_Ri(const double T) const { - double r0=0., - r1=Rc; - double dT0=get_T(T,r0)-get_solidus(get_X(r0),get_Pressure(r0)), - dT1=get_T(T,r1)-get_solidus(get_X(r1),get_Pressure(r1)); + double r0 = 0.; + double r1 = Rc; + const double dT0 = get_T(T,r0) - get_solidus(get_X(r0),get_pressure(r0)); + const double dT1 = get_T(T,r1) - get_solidus(get_X(r1),get_pressure(r1)); + if (dT0<=0. && dT1<=0.) return Rc; if (dT0>=0. && dT1>=0.) return 0.; for (int i=0; i0 && dT1<0) @@ -401,7 +410,7 @@ namespace aspect template bool - DynamicCore::solve_time_step(double &X, double &T, double &R) + DynamicCore::solve_time_step(double &X, double &T, double &R) const { // When solving the change in core-mantle boundary temperature T, inner core radius R, and // light component (e.g. S, O, Si) composition X, the following relations has to be respected: @@ -419,16 +428,15 @@ namespace aspect // converge well. Instead, we solve the inner core radius by bisection method. int steps=1; - double R_0,R_1,R_2; - // dT is the temperature difference between adiabatic and solidus at + + // dT is the temperature difference between adiabatic temperature and solidus at // inner-outer core boundary. If dT=0 then we found our solution. - double dT0,dT1,dT2; - R_0 = 0.; - R_1 = core_data.Ri; - R_2 = Rc; - dT0 = get_dT(R_0); - dT1 = get_dT(R_1); - dT2 = get_dT(R_2); + double R_0 = 0.; + double R_1 = core_data.Ri; + double R_2 = Rc; + double dT0 = get_dT(R_0); + double dT1 = get_dT(R_1); + double dT2 = get_dT(R_2); if (dT0 >= 0. && dT2 >= 0.) { @@ -498,9 +506,11 @@ namespace aspect return false; } + + template double - DynamicCore::get_Tc(double r) const + DynamicCore::get_Tc(const double r) const { // Using all Q values from last step. // Qs & Qr is constant, while Qg & Ql depends on inner core radius Ri @@ -510,21 +520,26 @@ namespace aspect ) / core_data.Qs; } + + template double - DynamicCore::get_Ts(double r) const + DynamicCore::get_Ts(const double r) const { - return get_solidus(get_X(r),get_Pressure(r)); + return get_solidus(get_X(r),get_pressure(r)); } + template double - DynamicCore::get_dT(double r) const + DynamicCore::get_dT(const double r) const { return get_T(get_Tc(r),r) - get_Ts(r); } + + template void DynamicCore::update_core_data() @@ -537,6 +552,8 @@ namespace aspect get_heat_solution(core_data.Ti,core_data.Ri,core_data.Xi,core_data.Eh); } + + template const internal::CoreData & DynamicCore::get_core_data() const @@ -544,42 +561,49 @@ namespace aspect return core_data; } + + template double - DynamicCore::get_solidus(double X,double p) const + DynamicCore::get_solidus(const double X, const double pressure) const { if (use_bw11) { - // Change from weight percent to mole percent. - double x,x0=32./88.; - if (X double - DynamicCore::get_X(double r) const + DynamicCore::get_X(const double r) const { - double xi_3=std::pow(r/Rc,3); + const double xi_3 = Utilities::fixed_power<3>(r/Rc); return X_init/(1-xi_3+Delta*xi_3); } @@ -587,8 +611,8 @@ namespace aspect void DynamicCore::update() { - core_data.dt = this->get_timestep(); - core_data.H = get_radioheating_rate(); + core_data.dt = this->get_timestep(); + core_data.H = get_radioheating_rate(); // It's a bit tricky here. // Didn't use the initialize() function instead because the postprocess is initialized after boundary temperature. @@ -611,8 +635,8 @@ namespace aspect Plugins::get_plugin_as_type> (this->get_geometry_model()); Rc=spherical_shell_geometry.inner_radius(); - Mc=get_Mass(Rc); - P_Core=get_Pressure(0); + Mc=get_mass(Rc); + P_Core=get_pressure(0); // If the material model is incompressible, we have to get correction for the real core temperature if (this->get_adiabatic_conditions().is_initialized() && !this->get_material_model().is_compressible()) @@ -734,15 +758,13 @@ namespace aspect local_CMB_area += local_face_area; } // now communicate to get the global values - double global_CMB_flux; - double global_CMB_area; - global_CMB_flux = Utilities::MPI::sum (local_CMB_flux, this->get_mpi_communicator()); - global_CMB_area = Utilities::MPI::sum (local_CMB_area, this->get_mpi_communicator()); + const double global_CMB_flux = Utilities::MPI::sum (local_CMB_flux, this->get_mpi_communicator()); + const double global_CMB_area = Utilities::MPI::sum (local_CMB_area, this->get_mpi_communicator()); // Using area averaged heat-flux density times core mantle boundary area to calculate total heat-flux on the 3d sphere. // By doing this, using dynamic core evolution with geometry other than 3d spherical shell becomes possible. - double average_CMB_heatflux_density = global_CMB_flux / global_CMB_area; - core_data.Q = average_CMB_heatflux_density * 4. * M_PI * Rc * Rc; + const double average_CMB_heatflux_density = global_CMB_flux / global_CMB_area; + core_data.Q = average_CMB_heatflux_density * 4. * numbers::PI * Rc * Rc; } core_data.Q_OES = get_OES(this->get_time()); @@ -784,172 +806,207 @@ namespace aspect } } + + template double DynamicCore:: - get_Mass(double r) const + get_mass(const double r) const { - return 4.*M_PI*Rho_cen*(-std::pow(L,2)/2.*r*std::exp(-std::pow(r/L,2))+std::pow(L,3)/4.*std::sqrt(M_PI)*std::erf(r/L)); + return 4.*numbers::PI*Rho_cen*(-std::pow(L,2)/2.*r*std::exp(-std::pow(r/L,2))+std::pow(L,3)/4.*std::sqrt(numbers::PI)*std::erf(r/L)); } + + template double DynamicCore:: - fun_Sn(double B,double R,double n) const + fun_Sn(const double B, const double R, const double n) const { - double S=R/(2.*std::sqrt(M_PI)); - for (unsigned i=1; i<=n; ++i) - S+=B/std::sqrt(M_PI)*std::exp(-std::pow(i,2)/4.)/i*std::sinh(i*R/B); + double S=R/(2.*std::sqrt(numbers::PI)); + for (unsigned int i=1; i<=n; ++i) + S+=B/std::sqrt(numbers::PI)*std::exp(-std::pow(i,2)/4.)/i*std::sinh(i*R/B); return S; } + + template double DynamicCore:: - get_Pressure(double r) const + get_pressure(const double r) const { - return P_CMB-(4*M_PI*constants::big_g*std::pow(Rho_cen,2))/3 + return P_CMB-(4*numbers::PI*constants::big_g*std::pow(Rho_cen,2))/3 *((3*std::pow(r,2)/10.-std::pow(L,2)/5)*std::exp(-std::pow(r/L,2)) -(3*std::pow(Rc,2)/10-std::pow(L,2)/5)*std::exp(-std::pow(Rc/L,2))); } + + template double DynamicCore:: - get_Rho(double r) const + get_rho(const double r) const { return Rho_cen*std::exp(-std::pow(r/L,2)); } + + template double DynamicCore:: - get_g(double r) const + get_g(const double r) const { - return (4*M_PI/3)*constants::big_g*Rho_cen*r*(1-3*std::pow(r,2)/(5*std::pow(L,2))); + return (4*numbers::PI/3)*constants::big_g*Rho_cen*r*(1-3*std::pow(r,2)/(5*std::pow(L,2))); } + + template double DynamicCore:: - get_T(double Tc,double r) const + get_T(const double Tc, const double r) const { return Tc*std::exp((std::pow(Rc,2)-std::pow(r,2))/std::pow(D,2)); } + + template double DynamicCore:: - get_gravity_potential(double r) const + get_gravity_potential(const double r) const { - return 2./3.*M_PI*constants::big_g*Rho_cen*(std::pow(r,2)*(1.-3.*std::pow(r,2) - /(10.*std::pow(L,2)))-std::pow(Rc,2)*(1.-3.*std::pow(Rc,2)/(10.*std::pow(L,2)))); + return 2./3.*numbers::PI*constants::big_g*Rho_cen*(std::pow(r,2)*(1.-3.*std::pow(r,2) + /(10.*std::pow(L,2)))-std::pow(Rc,2)*(1.-3.*std::pow(Rc,2)/(10.*std::pow(L,2)))); } + + template void DynamicCore:: - get_specific_heating(double Tc, double &Qs,double &Es) + get_specific_heating(const double Tc, double &Qs,double &Es) const { - double A=std::sqrt(1./(std::pow(L,-2)+std::pow(D,-2))); - double Is=4.*M_PI*get_T(Tc,0.)*Rho_cen*(-std::pow(A,2)*Rc/2.*std::exp(-std::pow(Rc/A,2))+std::pow(A,3)*std::sqrt(M_PI)/4.*std::erf(Rc/A)); + const double A = std::sqrt(1./(std::pow(L,-2)+std::pow(D,-2))); + const double Is = 4.*numbers::PI*get_T(Tc,0.)*Rho_cen*(-std::pow(A,2)*Rc/2.*std::exp(-std::pow(Rc/A,2))+std::pow(A,3)*std::sqrt(numbers::PI)/4.*std::erf(Rc/A)); Qs=-Cp/Tc*Is; Es=Cp/Tc*(Mc-Is/Tc); } + + template void DynamicCore:: - get_radio_heating(double Tc, double &Qr, double &Er) + get_radio_heating(const double Tc, double &Qr, double &Er) const { double B,It; if (D>L) { B=std::sqrt(1/(1/std::pow(L,2)-1/std::pow(D,2))); - It=4*M_PI*Rho_cen/get_T(Tc,0)*(-std::pow(B,2)*Rc/2*std::exp(-std::pow(Rc/B,2))+std::pow(B,3)/std::sqrt(M_PI)/4*std::erf(Rc/B)); + It=4*numbers::PI*Rho_cen/get_T(Tc,0)*(-std::pow(B,2)*Rc/2*std::exp(-std::pow(Rc/B,2))+std::pow(B,3)/std::sqrt(numbers::PI)/4*std::erf(Rc/B)); } else { B=std::sqrt(1/(std::pow(D,-2)-std::pow(L,-2))); - It=4*M_PI*Rho_cen/get_T(Tc,0)*(std::pow(B,2)*Rc/2*std::exp(std::pow(Rc/B,2))-std::pow(B,2)*fun_Sn(B,Rc,100)/2); + It=4*numbers::PI*Rho_cen/get_T(Tc,0)*(std::pow(B,2)*Rc/2*std::exp(std::pow(Rc/B,2))-std::pow(B,2)*fun_Sn(B,Rc,100)/2); } + Qr=Mc*core_data.H; Er=(Mc/Tc-It)*core_data.H; - } + + template void DynamicCore:: - get_heat_solution(double Tc, double r, double X,double &Eh) + get_heat_solution(const double Tc, const double r, const double X, double &Eh) const { double B,It; if (D>L) { B=std::sqrt(1/(1/std::pow(L,2)-1/std::pow(D,2))); - It=4*M_PI*Rho_cen/get_T(Tc,0)*(-std::pow(B,2)*Rc/2*std::exp(-std::pow(Rc/B,2))+std::pow(B,3)/std::sqrt(M_PI)/4*std::erf(Rc/B)); - It-=4*M_PI*Rho_cen/get_T(Tc,0)*(-std::pow(B,2)*r/2*std::exp(-std::pow(r/B,2))+std::pow(B,3)/std::sqrt(M_PI)/4*std::erf(r/B)); + It = 4*numbers::PI*Rho_cen/get_T(Tc,0)*(-std::pow(B,2)*Rc/2*std::exp(-std::pow(Rc/B,2))+std::pow(B,3)/std::sqrt(numbers::PI)/4*std::erf(Rc/B)); + It -= 4*numbers::PI*Rho_cen/get_T(Tc,0)*(-std::pow(B,2)*r/2*std::exp(-std::pow(r/B,2))+std::pow(B,3)/std::sqrt(numbers::PI)/4*std::erf(r/B)); } else { - B=std::sqrt(1/(std::pow(D,-2)-std::pow(L,-2))); - It=4*M_PI*Rho_cen/get_T(Tc,0)*(std::pow(B,2)*Rc/2*std::exp(std::pow(Rc/B,2))-std::pow(B,2)*fun_Sn(B,Rc,100)/2); - It-=4*M_PI*Rho_cen/get_T(Tc,0)*(std::pow(B,2)*r/2*std::exp(std::pow(r/B,2))-std::pow(B,2)*fun_Sn(B,r,100)/2); + B = std::sqrt(1/(std::pow(D,-2)-std::pow(L,-2))); + It = 4*numbers::PI*Rho_cen/get_T(Tc,0)*(std::pow(B,2)*Rc/2*std::exp(std::pow(Rc/B,2))-std::pow(B,2)*fun_Sn(B,Rc,100)/2); + It -= 4*numbers::PI*Rho_cen/get_T(Tc,0)*(std::pow(B,2)*r/2*std::exp(std::pow(r/B,2))-std::pow(B,2)*fun_Sn(B,r,100)/2); } - double Cc=4*M_PI*std::pow(r,2)*get_Rho(r)*X/(Mc-get_Mass(r)); - Eh=Rh*(It-(Mc-get_Mass(r))/get_T(Tc,r))*Cc; + const double Cc = 4*numbers::PI*std::pow(r,2)*get_rho(r)*X/(Mc-get_mass(r)); + Eh = Rh*(It-(Mc-get_mass(r))/get_T(Tc,r))*Cc; } + + template void DynamicCore:: - get_gravity_heating(double Tc, double r,double X,double &Qg,double &Eg) + get_gravity_heating(const double Tc, const double r, const double X, double &Qg, double &Eg) const { - double Cc=4*M_PI*std::pow(r,2)*get_Rho(r)*X/(Mc-get_Mass(r)); - double C_2=3./16.*std::pow(L,2)-0.5*std::pow(Rc,2)*(1.-3./10.*std::pow(Rc/L,2)); + const double Cc = 4*numbers::PI*std::pow(r,2)*get_rho(r)*X/(Mc-get_mass(r)); + const double C_2 = 3./16.*std::pow(L,2) - 0.5*std::pow(Rc,2)*(1.-3./10.*std::pow(Rc/L,2)); if (r==Rc) Qg=0.; else - Qg=(8./3.*std::pow(M_PI*Rho_cen,2)*constants::big_g*( - ((3./20.*std::pow(Rc,5)-std::pow(L,2)*std::pow(Rc,3)/8.-C_2*std::pow(L,2)*Rc)*std::exp(-std::pow(Rc/L,2)) - +C_2/2.*std::pow(L,3)*std::sqrt(M_PI)*std::erf(Rc/L)) - -((3./20.*std::pow(r,5)-std::pow(L,2)*std::pow(r,3)/8.-C_2*std::pow(L,2)*r)*std::exp(-std::pow(r/L,2)) - +C_2/2.*std::pow(L,3)*std::sqrt(M_PI)*std::erf(r/L))) - -(Mc-get_Mass(r))*get_gravity_potential(r))*Beta_c*Cc; - Eg=Qg/Tc; + { + Qg=(8./3.*std::pow(numbers::PI*Rho_cen,2)*constants::big_g*( + ((3./20.*std::pow(Rc,5)-std::pow(L,2)*std::pow(Rc,3)/8.-C_2*std::pow(L,2)*Rc)*std::exp(-std::pow(Rc/L,2)) + +C_2/2.*std::pow(L,3)*std::sqrt(numbers::PI)*std::erf(Rc/L)) + -((3./20.*std::pow(r,5)-std::pow(L,2)*std::pow(r,3)/8.-C_2*std::pow(L,2)*r)*std::exp(-std::pow(r/L,2)) + +C_2/2.*std::pow(L,3)*std::sqrt(numbers::PI)*std::erf(r/L))) + -(Mc-get_mass(r))*get_gravity_potential(r))*Beta_c*Cc; + } + Eg = Qg/Tc; } + + template void DynamicCore:: - get_adiabatic_heating(double Tc, double &Ek, double &Qk) + get_adiabatic_heating(const double Tc, double &Ek, double &Qk) const { - Ek=16*M_PI*k_c*std::pow(Rc,5)/5/std::pow(D,4); - Qk=8*M_PI*std::pow(Rc,3)*k_c*Tc/std::pow(D,2); + Ek = 16*numbers::PI*k_c*std::pow(Rc,5)/5/std::pow(D,4); + Qk = 8*numbers::PI*std::pow(Rc,3)*k_c*Tc/std::pow(D,2); } + + + template void DynamicCore:: - get_latent_heating(double Tc, double r, double &El, double &Ql) + get_latent_heating(const double Tc, const double r, double &El, double &Ql) const { - Ql=4.*M_PI*std::pow(r,2)*Lh*get_Rho(r); - El=Ql*(get_T(Tc,r)-Tc)/(Tc*get_T(Tc,r)); + Ql = 4.*numbers::PI*std::pow(r,2)*Lh*get_rho(r); + El = Ql*(get_T(Tc,r)-Tc)/(Tc*get_T(Tc,r)); } + + template double DynamicCore:: get_radioheating_rate() const { - double time=this->get_time()+0.5*this->get_timestep(); + const double time=this->get_time()+0.5*this->get_timestep(); + double Ht=0; for (unsigned i=0; i bool DynamicCore:: @@ -961,7 +1018,6 @@ namespace aspect return false; } } - }