diff --git a/include/aspect/boundary_temperature/dynamic_core.h b/include/aspect/boundary_temperature/dynamic_core.h index e9d4dc25a54..11c4d9cd457 100644 --- a/include/aspect/boundary_temperature/dynamic_core.h +++ b/include/aspect/boundary_temperature/dynamic_core.h @@ -344,7 +344,7 @@ namespace aspect /** * Max iterations for the core energy balance solver. */ - int max_steps; + unsigned int max_steps; /** * Temperature correction value for adiabatic @@ -409,8 +409,8 @@ namespace aspect double get_Ts(const double r) const; /** - * Compute the core solidus at certain pressure @p pressure - * and light element concentration @p X (in wt.%). + * Compute the core solidus at a given light element concentration @p X (in wt.%) + * and pressure @p pressure. */ double get_solidus(const double X, const double pressure) const; @@ -434,7 +434,7 @@ namespace aspect /** * Calculate Sn(B,R), referring to \cite NPB+04 . */ - double fun_Sn(const double B, const double R, const double n) const; + double fun_Sn(const double B, const double R, const unsigned int n) const; /** * Calculate density at given radius @p r. diff --git a/source/boundary_temperature/dynamic_core.cc b/source/boundary_temperature/dynamic_core.cc index 8c50696fca3..593a9f41803 100644 --- a/source/boundary_temperature/dynamic_core.cc +++ b/source/boundary_temperature/dynamic_core.cc @@ -38,275 +38,232 @@ namespace aspect { namespace BoundaryTemperature { - template - double - DynamicCore:: - boundary_temperature (const types::boundary_id boundary_indicator, - const Point &/*location*/) const + DynamicCore::DynamicCore() { - switch (boundary_indicator) - { - case 0: - return inner_temperature; - case 1: - return outer_temperature; - default: - Assert (false, ExcMessage ("Unknown boundary indicator.")); - return std::numeric_limits::quiet_NaN(); - } + is_first_call = true; + core_data.is_initialized = false; } template - double - DynamicCore:: - minimal_temperature (const std::set &/*fixed_boundary_ids*/) const + void + DynamicCore::update() { - return std::min (inner_temperature, outer_temperature); - } + 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. + // It is not available at the time initialize() function of boundary temperature is called. + if (is_first_call==true) + { + AssertThrow(this->get_postprocess_manager().template has_matching_active_plugin>(), + ExcMessage ("Dynamic core boundary condition has to work with dynamic core statistics postprocessor.")); + const Postprocess::CoreStatistics &core_statistics + = this->get_postprocess_manager().template get_matching_active_plugin>(); + // The restart data is stored in 'core statistics' postprocessor. + // If restart from checkpoint, extract data from there. + core_data = core_statistics.get_core_data(); - template - double - DynamicCore:: - maximal_temperature (const std::set &/*fixed_boundary_ids*/) const - { - return std::max (inner_temperature, outer_temperature); - } + // Read data of other energy source + read_data_OES(); + const GeometryModel::SphericalShell &spherical_shell_geometry = + Plugins::get_plugin_as_type> (this->get_geometry_model()); + Rc = spherical_shell_geometry.inner_radius(); + Mc = get_mass(Rc); + P_Core = get_pressure(0); - template - void - DynamicCore::declare_parameters (ParameterHandler &prm) - { - prm.enter_subsection("Boundary temperature model"); + // 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()) + { + Point p1; + p1(0) = spherical_shell_geometry.inner_radius(); + dTa = this->get_adiabatic_conditions().temperature(p1) + - this->get_adiabatic_surface_temperature(); + } + else + dTa = 0.; + + // Setup initial core data from prm input. + // If resumed from checkpoint, core_data is read from postprocess instead of set from prm file. + // (The boundary_temperature doesn't seem to support restart/resume, the data has to passed and + // stored in the postprocessor 'core statistics') + if (!core_data.is_initialized) + { + core_data.Ti = inner_temperature + dTa; + core_data.Ri = get_initial_Ri(core_data.Ti); + core_data.Xi = get_X(core_data.Ri); + + core_data.Q = 0.; + core_data.dt = 0.; + core_data.dT_dt = init_dT_dt; + core_data.dR_dt = init_dR_dt; + core_data.dX_dt = init_dX_dt; + update_core_data(); + core_data.is_initialized = true; + std::stringstream output; + output<get_pcout() << output.str(); + } + is_first_call = false; + } + + // Calculate core mantle boundary heat flow { - prm.enter_subsection("Dynamic core"); - { - prm.declare_entry ("Outer temperature", "0.", - Patterns::Double (), - "Temperature at the outer boundary (lithosphere water/air). Units: \\si{\\kelvin}."); - prm.declare_entry ("Inner temperature", "6000.", - Patterns::Double (), - "Temperature at the inner boundary (core mantle boundary) at the " - "beginning. Units: \\si{\\kelvin}."); - prm.declare_entry ("dT over dt", "0.", - Patterns::Double (), - "Initial CMB temperature changing rate. " - "Units: \\si{\\kelvin}/year."); - prm.declare_entry ("dR over dt", "0.", - Patterns::Double (), - "Initial inner core radius changing rate. " - "Units: \\si{\\kilo\\meter}/year."); - prm.declare_entry ("dX over dt", "0.", - Patterns::Double (), - "Initial light composition changing rate. " - "Units: 1/year."); - prm.declare_entry ("Core density", "12.5e3", - Patterns::Double (), - "Density of the core. " - "Units: \\si{\\kilogram\\per\\meter\\cubed}."); - prm.declare_entry ("Gravity acceleration", "9.8", - Patterns::Double (), - "Gravitation acceleration at CMB. " - "Units: \\si{\\meter\\per\\second\\squared}."); - prm.declare_entry ("CMB pressure", "0.14e12", - Patterns::Double (), - "Pressure at CMB. Units: \\si{\\pascal}."); - prm.declare_entry ("Initial light composition", "0.01", - Patterns::Double (0.), - "Initial light composition (eg. S,O) concentration " - "in weight fraction."); - prm.declare_entry ("Max iteration", "30000", - Patterns::Integer (0), - "The max iterations for nonlinear core energy solver."); - prm.declare_entry ("Core heat capacity", "840.", - Patterns::Double (0.), - "Heat capacity of the core. " - "Units: \\si{\\joule\\per\\kelvin\\per\\kilogram}."); - prm.declare_entry ("K0", "4.111e11", - Patterns::Double (0.), - "Core compressibility at zero pressure. " - "See \\cite{NPB+04} for more details."); - prm.declare_entry ("Rho0", "7.019e3", - Patterns::Double (0.), - "Core density at zero pressure. " - "Units: \\si{\\kilogram\\per\\meter\\cubed}. " - "See \\cite{NPB+04} for more details."); - prm.declare_entry ("Alpha", "1.35e-5", - Patterns::Double (0.), - "Core thermal expansivity. Units: \\si{\\per\\kelvin}."); - prm.declare_entry ("Lh", "750e3", - Patterns::Double (0.), - "The latent heat of core freeze. " - "Units: \\si{\\joule\\per\\kilogram}."); - prm.declare_entry ("Rh","-27.7e6", - Patterns::Double (), - "The heat of reaction. " - "Units: \\si{\\joule\\per\\kilogram}."); - prm.declare_entry ("Beta composition", "1.1", - Patterns::Double (0.), - "Compositional expansion coefficient $Beta_c$. " - "See \\cite{NPB+04} for more details."); - prm.declare_entry ("Delta","0.5", - Patterns::Double (0., 1.), - "Partition coefficient of the light element."); - prm.declare_entry ("Core conductivity", "60.", - Patterns::Double (0.), - "Core heat conductivity $k_c$. Units: \\si{\\watt\\per\\meter\\per\\kelvin}."); - prm.enter_subsection("Geotherm parameters"); - { - prm.declare_entry ("Tm0","1695.", - Patterns::Double (0.), - "Melting curve (\\cite{NPB+04} eq. (40)) parameter Tm0. Units: \\si{\\kelvin}."); - prm.declare_entry ("Tm1","10.9", - Patterns::Double (), - "Melting curve (\\cite{NPB+04} eq. (40)) parameter Tm1. " - "Units: \\si{\\per\\tera\\pascal}."); - prm.declare_entry ("Tm2","-8.0", - Patterns::Double (), - "Melting curve (\\cite{NPB+04} eq. (40)) parameter Tm2. " - "Units: \\si{\\per\\tera\\pascal\\squared}."); - prm.declare_entry ("Theta","0.11", - Patterns::Double (), - "Melting curve (\\cite{NPB+04} eq. (40)) parameter Theta."); - prm.declare_entry ("Composition dependency","true", - Patterns::Bool (), - "If melting curve dependent on composition."); - prm.declare_entry ("Use BW11","false", - Patterns::Bool (), - "If using the Fe-FeS system solidus from Buono \\& Walker (2011) instead."); - } - prm.leave_subsection (); + const Quadrature &quadrature_formula = this->introspection().face_quadratures.temperature; + FEFaceValues fe_face_values (this->get_mapping(), + this->get_fe(), + quadrature_formula, + update_gradients | update_values | + update_normal_vectors | + update_quadrature_points | update_JxW_values); - prm.enter_subsection("Radioactive heat source"); - { - prm.declare_entry ("Number of radioactive heating elements","0", - Patterns::Integer (0), - "Number of different radioactive heating elements in core"); - prm.declare_entry ("Heating rates","", - Patterns::List (Patterns::Double ()), - "Heating rates of different elements (W/kg)"); - prm.declare_entry ("Half life times","", - Patterns::List (Patterns::Double ()), - "Half decay times of different elements (Ga)"); - prm.declare_entry ("Initial concentrations","", - Patterns::List (Patterns::Double ()), - "Initial concentrations of different elements (ppm)"); - } - prm.leave_subsection (); + std::vector> temperature_gradients (quadrature_formula.size()); + std::vector> composition_values (this->n_compositional_fields(),std::vector (quadrature_formula.size())); - prm.enter_subsection("Other energy source"); - { - prm.declare_entry ("File name","", - Patterns::Anything(), - "Data file name for other energy source into the core. " - "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.)" - "Format [Time(Gyr) Energy rate(W)]"); - } - prm.leave_subsection (); + //std::map local_boundary_fluxes; + double local_CMB_flux = 0.; + double local_CMB_area = 0.; - } - prm.leave_subsection (); - } - prm.leave_subsection (); - } + types::boundary_id CMB_id = 0; + typename MaterialModel::Interface::MaterialModelInputs in(fe_face_values.n_quadrature_points, this->n_compositional_fields()); + typename MaterialModel::Interface::MaterialModelOutputs out(fe_face_values.n_quadrature_points, this->n_compositional_fields()); + // Do not request viscosity or reaction rates + in.requested_properties = MaterialModel::MaterialProperties::equation_of_state_properties | + MaterialModel::MaterialProperties::thermal_conductivity; + // for every surface face on which it makes sense to compute a + // heat flux and that is owned by this processor, + // integrate the normal heat flux given by the formula + // j = - k * n . grad T + // + // for the spherical shell geometry, note that for the inner boundary, + // the normal vector points *into* the core, i.e. we compute the flux + // *out* of the mantle, not into it. we fix this when we add the local + // contribution to the global flux - template - void - DynamicCore::parse_parameters (ParameterHandler &prm) - { - prm.enter_subsection("Boundary temperature model"); - { - prm.enter_subsection("Dynamic core"); + for (const auto &cell : this->get_dof_handler().active_cell_iterators()) + if (cell->is_locally_owned()) + for (const unsigned int f : cell->face_indices()) + if (cell->at_boundary(f)) + if (cell->face(f)->boundary_id() == CMB_id) + { + fe_face_values.reinit (cell, f); + + in.reinit(fe_face_values, cell, this->introspection(), this->get_solution()); + + fe_face_values[this->introspection().extractors.temperature].get_function_gradients (this->get_solution(), + temperature_gradients); + + this->get_material_model().evaluate(in, out); + + + double local_normal_flux = 0; + double local_face_area = 0; + for (unsigned int q=0; qget_material_model().is_compressible()==false) + { + const double alpha = out.thermal_expansion_coefficients[q]; + const double cp = out.specific_heat[0]; + const double gravity = this->get_gravity_model().gravity_vector(in.position[q]).norm(); + if (cell->face(f)->boundary_id()==0) + adiabatic_flux = - alpha * gravity / cp; + else if (cell->face(f)->boundary_id()==1) + adiabatic_flux = alpha * gravity / cp; + } + + local_normal_flux += -thermal_conductivity * + (temperature_gradients[q] * fe_face_values.normal_vector(q) + + adiabatic_flux) * fe_face_values.JxW(q); + local_face_area += fe_face_values.JxW(q); + + } + local_CMB_flux += local_normal_flux; + local_CMB_area += local_face_area; + } + // now communicate to get the global values + 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. + 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()); + + if ((core_data.Q + core_data.Q_OES) * core_data.dt!=0.) { - // verify that the geometry is in fact a spherical shell since only - // for this geometry do we know for sure what boundary indicators it - // uses and what they mean - AssertThrow (Plugins::plugin_type_matches>(this->get_geometry_model()), - ExcMessage ("This boundary model is only implemented if the geometry is " - "a spherical shell.")); + double X1,R1 = core_data.Ri,T1; + solve_time_step(X1,T1,R1); + if (core_data.dt != 0) + { + core_data.dR_dt = (R1-core_data.Ri)/core_data.dt; + core_data.dT_dt = (T1-core_data.Ti)/core_data.dt; + core_data.dX_dt = (X1-core_data.Xi)/core_data.dt; + } + else + { + core_data.dR_dt = 0.; + core_data.dT_dt = 0.; + core_data.dX_dt = 0.; + } + core_data.Xi = X1; + core_data.Ri = R1; + core_data.Ti = T1; + } - inner_temperature = prm.get_double ("Inner temperature"); - outer_temperature = prm.get_double ("Outer temperature"); - init_dT_dt = prm.get_double ("dT over dt") / year_in_seconds; - init_dR_dt = prm.get_double ("dR over dt") / year_in_seconds * 1.e3; - init_dX_dt = prm.get_double ("dX over dt") / year_in_seconds; - Rho_cen = prm.get_double ("Core density"); - g = prm.get_double ("Gravity acceleration"); - P_CMB = prm.get_double ("CMB pressure"); - X_init = prm.get_double ("Initial light composition"); - max_steps = prm.get_integer ("Max iteration"); - Cp = prm.get_double ("Core heat capacity"); - CpRho = Cp*Rho_cen; + inner_temperature = core_data.Ti - dTa; + update_core_data(); + if ((core_data.Q + core_data.Q_OES + core_data.Qr) * core_data.dt != 0.) + { + std::stringstream output; + output<get_pcout() << output.str(); + } + } - //\cite{NPB+04} - K0 = prm.get_double ("K0"); - Alpha = prm.get_double ("Alpha"); - Rho_0 = prm.get_double ("Rho0"); - Lh = prm.get_double ("Lh"); - Beta_c = prm.get_double ("Beta composition"); - k_c = prm.get_double ("Core conductivity"); - Delta = prm.get_double ("Delta"); - Rh = prm.get_double ("Rh"); - prm.enter_subsection("Geotherm parameters"); - { - Tm0 = prm.get_double ("Tm0"); - Tm1 = prm.get_double ("Tm1"); - Tm2 = prm.get_double ("Tm2"); - Theta = prm.get_double ("Theta"); - composition_dependency - = prm.get_bool("Composition dependency"); - use_bw11 = prm.get_bool("Use BW11"); - } - prm.leave_subsection (); - prm.enter_subsection("Radioactive heat source"); - { - n_radioheating_elements = prm.get_integer ("Number of radioactive heating elements"); - heating_rate = Utilities::string_to_double - (Utilities::split_string_list - (prm.get("Heating rates"))); - AssertThrow(n_radioheating_elements==heating_rate.size(), - ExcMessage("Number of heating rate entities does not match " - "the number of radioactive elements.")); - half_life = Utilities::string_to_double - (Utilities::split_string_list - (prm.get("Half life times"))); - AssertThrow(n_radioheating_elements==half_life.size(), - ExcMessage("Number of half life time entities does not match " - "the number of radioactive elements.")); - initial_concentration = Utilities::string_to_double - (Utilities::split_string_list - (prm.get("Initial concentrations"))); - AssertThrow(n_radioheating_elements==initial_concentration.size(), - ExcMessage("Number of initial concentration entities does not match " - "the number of radioactive elements.")); - } - prm.leave_subsection (); + template + double + DynamicCore:: + minimal_temperature (const std::set &/*fixed_boundary_ids*/) const + { + return std::min (inner_temperature, outer_temperature); + } - prm.enter_subsection("Other energy source"); - { - name_OES = prm.get("File name"); - } - prm.leave_subsection (); - 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 (); - } - prm.leave_subsection (); + template + double + DynamicCore:: + maximal_temperature (const std::set &/*fixed_boundary_ids*/) const + { + return std::max (inner_temperature, outer_temperature); } @@ -316,7 +273,7 @@ namespace aspect DynamicCore::read_data_OES() { data_OES.clear(); - if (name_OES.size()==0) + if (name_OES.size() == 0) return; std::istringstream in(Utilities::read_and_distribute_file_content(name_OES, this->get_mpi_communicator())); @@ -331,7 +288,7 @@ namespace aspect data_OES.push_back(data_read); } } - if (data_OES.size()!=0) + if (data_OES.size() != 0) this->get_pcout() << "Other energy source is in use ( " << data_OES.size() << " data points is read)." @@ -346,7 +303,7 @@ namespace aspect { // The core evolution is quite slow, so the time units used here is billion years. const double t = time / (1.e9*year_in_seconds); - double w=0.; + double w = 0.; for (unsigned i=1; i=data_OES[i-1].t && t - DynamicCore::DynamicCore() - { - is_first_call = true; - core_data.is_initialized = false; - } - - - template double DynamicCore::get_initial_Ri(const double T) const @@ -384,17 +332,17 @@ namespace aspect return Rc; if (dT0>=0. && dT1>=0.) return 0.; - for (int i=0; i double - DynamicCore::get_X(const double r) const - { - const double xi_3 = Utilities::fixed_power<3>(r/Rc); - return X_init/(1-xi_3+Delta*xi_3); - } - - template - void - DynamicCore::update() - { - 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. - // It is not available at the time initialize() function of boundary temperature is called. - if (is_first_call==true) - { - AssertThrow(this->get_postprocess_manager().template has_matching_active_plugin>(), - ExcMessage ("Dynamic core boundary condition has to work with dynamic core statistics postprocessor.")); - - const Postprocess::CoreStatistics &core_statistics - = this->get_postprocess_manager().template get_matching_active_plugin>(); - // The restart data is stored in 'core statistics' postprocessor. - // If restart from checkpoint, extract data from there. - core_data = core_statistics.get_core_data(); - - // Read data of other energy source - read_data_OES(); - - const GeometryModel::SphericalShell &spherical_shell_geometry = - Plugins::get_plugin_as_type> (this->get_geometry_model()); - - Rc=spherical_shell_geometry.inner_radius(); - 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()) - { - Point p1; - p1(0) = spherical_shell_geometry.inner_radius(); - dTa = this->get_adiabatic_conditions().temperature(p1) - - this->get_adiabatic_surface_temperature(); - } - else - dTa = 0.; - - // Setup initial core data from prm input. - // If resumed from checkpoint, core_data is read from postprocess instead of set from prm file. - // (The boundary_temperature doesn't seem to support restart/resume, the data has to passed and - // stored in the postprocessor 'core statistics') - if (!core_data.is_initialized) - { - core_data.Ti=inner_temperature + dTa; - core_data.Ri = get_initial_Ri(core_data.Ti); - core_data.Xi = get_X(core_data.Ri); - - core_data.Q=0.; - core_data.dt=0.; - core_data.dT_dt=init_dT_dt; - core_data.dR_dt=init_dR_dt; - core_data.dX_dt=init_dX_dt; - update_core_data(); - core_data.is_initialized = true; - std::stringstream output; - output<get_pcout() << output.str(); - } - is_first_call = false; - } - - // Calculate core mantle boundary heat flow - { - const Quadrature &quadrature_formula = this->introspection().face_quadratures.temperature; - FEFaceValues fe_face_values (this->get_mapping(), - this->get_fe(), - quadrature_formula, - update_gradients | update_values | - update_normal_vectors | - update_quadrature_points | update_JxW_values); - - std::vector> temperature_gradients (quadrature_formula.size()); - std::vector> composition_values (this->n_compositional_fields(),std::vector (quadrature_formula.size())); - - //std::map local_boundary_fluxes; - double local_CMB_flux = 0.; - double local_CMB_area = 0.; - - types::boundary_id CMB_id = 0; - - typename MaterialModel::Interface::MaterialModelInputs in(fe_face_values.n_quadrature_points, this->n_compositional_fields()); - typename MaterialModel::Interface::MaterialModelOutputs out(fe_face_values.n_quadrature_points, this->n_compositional_fields()); - // Do not request viscosity or reaction rates - in.requested_properties = MaterialModel::MaterialProperties::equation_of_state_properties | - MaterialModel::MaterialProperties::thermal_conductivity; - - // for every surface face on which it makes sense to compute a - // heat flux and that is owned by this processor, - // integrate the normal heat flux given by the formula - // j = - k * n . grad T - // - // for the spherical shell geometry, note that for the inner boundary, - // the normal vector points *into* the core, i.e. we compute the flux - // *out* of the mantle, not into it. we fix this when we add the local - // contribution to the global flux - - for (const auto &cell : this->get_dof_handler().active_cell_iterators()) - if (cell->is_locally_owned()) - for (const unsigned int f : cell->face_indices()) - if (cell->at_boundary(f)) - if (cell->face(f)->boundary_id() == CMB_id) - { - fe_face_values.reinit (cell, f); - - in.reinit(fe_face_values, cell, this->introspection(), this->get_solution()); - - fe_face_values[this->introspection().extractors.temperature].get_function_gradients (this->get_solution(), - temperature_gradients); - - this->get_material_model().evaluate(in, out); - - - double local_normal_flux = 0; - double local_face_area = 0; - for (unsigned int q=0; qget_material_model().is_compressible()==false) - { - const double alpha = out.thermal_expansion_coefficients[q]; - const double cp = out.specific_heat[0]; - const double gravity = this->get_gravity_model().gravity_vector(in.position[q]).norm(); - if (cell->face(f)->boundary_id()==0) - adiabatic_flux = - alpha * gravity / cp; - else if (cell->face(f)->boundary_id()==1) - adiabatic_flux = alpha * gravity / cp; - } - - local_normal_flux += -thermal_conductivity * - (temperature_gradients[q] * fe_face_values.normal_vector(q) - + adiabatic_flux) * fe_face_values.JxW(q); - local_face_area += fe_face_values.JxW(q); - - } - local_CMB_flux += local_normal_flux; - local_CMB_area += local_face_area; - } - // now communicate to get the global values - 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. - 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()); - - if ((core_data.Q + core_data.Q_OES) * core_data.dt!=0.) - { - double X1,R1=core_data.Ri,T1; - solve_time_step(X1,T1,R1); - if (core_data.dt!=0) - { - core_data.dR_dt=(R1-core_data.Ri)/core_data.dt; - core_data.dT_dt=(T1-core_data.Ti)/core_data.dt; - core_data.dX_dt=(X1-core_data.Xi)/core_data.dt; - } - else - { - core_data.dR_dt=0.; - core_data.dT_dt=0.; - core_data.dX_dt=0.; - } - core_data.Xi=X1; - core_data.Ri=R1; - core_data.Ti=T1; - } - - inner_temperature = core_data.Ti - dTa; - update_core_data(); - if ((core_data.Q + core_data.Q_OES + core_data.Qr) * core_data.dt!=0.) - { - std::stringstream output; - output<get_pcout() << output.str(); - } + DynamicCore::get_X(const double r) const + { + const double xi_3 = Utilities::fixed_power<3>(r/Rc); + return X_init/(1-xi_3+Delta*xi_3); } @@ -821,11 +570,17 @@ namespace aspect template double DynamicCore:: - fun_Sn(const double B, const double R, const double n) const + fun_Sn(const double B, const double R, const unsigned int n) const { - double S=R/(2.*std::sqrt(numbers::PI)); + // TODO: sqrt_pi could be made constexpr, but std::sqrt is not a constexpr function + // for the MacOS tester at the moment + const double sqrt_pi = std::sqrt(numbers::PI); + double S = R/(2.*sqrt_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); + { + const double it = static_cast(i); + S += (B/sqrt_pi) * (std::exp(-it*it/4.)/it) * std::sinh(it*R/B); + } return S; } @@ -892,8 +647,8 @@ namespace aspect 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); + Qs = -Cp/Tc*Is; + Es = Cp/Tc*(Mc-Is/Tc); } @@ -903,20 +658,20 @@ namespace aspect DynamicCore:: get_radio_heating(const double Tc, double &Qr, double &Er) const { - double B,It; + double It = numbers::signaling_nan(); if (D>L) { - B=std::sqrt(1/(1/std::pow(L,2)-1/std::pow(D,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,3)/std::sqrt(numbers::PI)/4*std::erf(Rc/B)); + const double B = std::sqrt(1/(1/std::pow(L,2)-1/std::pow(D,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,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*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); + const double 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); } - Qr=Mc*core_data.H; - Er=(Mc/Tc-It)*core_data.H; + Qr = Mc*core_data.H; + Er = (Mc/Tc-It)*core_data.H; } @@ -926,16 +681,16 @@ namespace aspect DynamicCore:: get_heat_solution(const double Tc, const double r, const double X, double &Eh) const { - double B,It; + double It = numbers::signaling_nan(); if (D>L) { - B=std::sqrt(1/(1/std::pow(L,2)-1/std::pow(D,2))); + const double B = std::sqrt(1./(1./std::pow(L,2)-1./std::pow(D,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,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))); + const double 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); } @@ -953,15 +708,15 @@ namespace aspect 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.; + Qg = 0.; else { - 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; + 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; @@ -996,11 +751,11 @@ namespace aspect DynamicCore:: get_radioheating_rate() const { - 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; + double Ht = 0; for (unsigned i=0; i + double + DynamicCore:: + boundary_temperature (const types::boundary_id boundary_indicator, + const Point &/*location*/) const + { + switch (boundary_indicator) + { + case 0: + return inner_temperature; + case 1: + return outer_temperature; + default: + Assert (false, ExcMessage ("Unknown boundary indicator.")); + } + + return std::numeric_limits::quiet_NaN(); + } + + + + template + void + DynamicCore::declare_parameters (ParameterHandler &prm) + { + prm.enter_subsection("Boundary temperature model"); + { + prm.enter_subsection("Dynamic core"); + { + prm.declare_entry ("Outer temperature", "0.", + Patterns::Double (), + "Temperature at the outer boundary (lithosphere water/air). Units: \\si{\\kelvin}."); + prm.declare_entry ("Inner temperature", "6000.", + Patterns::Double (), + "Temperature at the inner boundary (core mantle boundary) at the " + "beginning. Units: \\si{\\kelvin}."); + prm.declare_entry ("dT over dt", "0.", + Patterns::Double (), + "Initial CMB temperature changing rate. " + "Units: \\si{\\kelvin}/year."); + prm.declare_entry ("dR over dt", "0.", + Patterns::Double (), + "Initial inner core radius changing rate. " + "Units: \\si{\\kilo\\meter}/year."); + prm.declare_entry ("dX over dt", "0.", + Patterns::Double (), + "Initial light composition changing rate. " + "Units: 1/year."); + prm.declare_entry ("Core density", "12.5e3", + Patterns::Double (), + "Density of the core. " + "Units: \\si{\\kilogram\\per\\meter\\cubed}."); + prm.declare_entry ("Gravity acceleration", "9.8", + Patterns::Double (), + "Gravitation acceleration at CMB. " + "Units: \\si{\\meter\\per\\second\\squared}."); + prm.declare_entry ("CMB pressure", "0.14e12", + Patterns::Double (), + "Pressure at CMB. Units: \\si{\\pascal}."); + prm.declare_entry ("Initial light composition", "0.01", + Patterns::Double (0.), + "Initial light composition (eg. S,O) concentration " + "in weight fraction."); + prm.declare_entry ("Max iteration", "30000", + Patterns::Integer (0), + "The max iterations for nonlinear core energy solver."); + prm.declare_entry ("Core heat capacity", "840.", + Patterns::Double (0.), + "Heat capacity of the core. " + "Units: \\si{\\joule\\per\\kelvin\\per\\kilogram}."); + prm.declare_entry ("K0", "4.111e11", + Patterns::Double (0.), + "Core compressibility at zero pressure. " + "See \\cite{NPB+04} for more details."); + prm.declare_entry ("Rho0", "7.019e3", + Patterns::Double (0.), + "Core density at zero pressure. " + "Units: \\si{\\kilogram\\per\\meter\\cubed}. " + "See \\cite{NPB+04} for more details."); + prm.declare_entry ("Alpha", "1.35e-5", + Patterns::Double (0.), + "Core thermal expansivity. Units: \\si{\\per\\kelvin}."); + prm.declare_entry ("Lh", "750e3", + Patterns::Double (0.), + "The latent heat of core freeze. " + "Units: \\si{\\joule\\per\\kilogram}."); + prm.declare_entry ("Rh","-27.7e6", + Patterns::Double (), + "The heat of reaction. " + "Units: \\si{\\joule\\per\\kilogram}."); + prm.declare_entry ("Beta composition", "1.1", + Patterns::Double (0.), + "Compositional expansion coefficient $Beta_c$. " + "See \\cite{NPB+04} for more details."); + prm.declare_entry ("Delta","0.5", + Patterns::Double (0., 1.), + "Partition coefficient of the light element."); + prm.declare_entry ("Core conductivity", "60.", + Patterns::Double (0.), + "Core heat conductivity $k_c$. Units: \\si{\\watt\\per\\meter\\per\\kelvin}."); + prm.enter_subsection("Geotherm parameters"); + { + prm.declare_entry ("Tm0","1695.", + Patterns::Double (0.), + "Melting curve (\\cite{NPB+04} eq. (40)) parameter Tm0. Units: \\si{\\kelvin}."); + prm.declare_entry ("Tm1","10.9", + Patterns::Double (), + "Melting curve (\\cite{NPB+04} eq. (40)) parameter Tm1. " + "Units: \\si{\\per\\tera\\pascal}."); + prm.declare_entry ("Tm2","-8.0", + Patterns::Double (), + "Melting curve (\\cite{NPB+04} eq. (40)) parameter Tm2. " + "Units: \\si{\\per\\tera\\pascal\\squared}."); + prm.declare_entry ("Theta","0.11", + Patterns::Double (), + "Melting curve (\\cite{NPB+04} eq. (40)) parameter Theta."); + prm.declare_entry ("Composition dependency","true", + Patterns::Bool (), + "If melting curve dependent on composition."); + prm.declare_entry ("Use BW11","false", + Patterns::Bool (), + "If using the Fe-FeS system solidus from Buono \\& Walker (2011) instead."); + } + prm.leave_subsection (); + + prm.enter_subsection("Radioactive heat source"); + { + prm.declare_entry ("Number of radioactive heating elements","0", + Patterns::Integer (0), + "Number of different radioactive heating elements in core"); + prm.declare_entry ("Heating rates","", + Patterns::List (Patterns::Double ()), + "Heating rates of different elements (W/kg)"); + prm.declare_entry ("Half life times","", + Patterns::List (Patterns::Double ()), + "Half decay times of different elements (Ga)"); + prm.declare_entry ("Initial concentrations","", + Patterns::List (Patterns::Double ()), + "Initial concentrations of different elements (ppm)"); + } + prm.leave_subsection (); + + prm.enter_subsection("Other energy source"); + { + prm.declare_entry ("File name","", + Patterns::Anything(), + "Data file name for other energy source into the core. " + "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.)" + "Format [Time(Gyr) Energy rate(W)]"); + } + prm.leave_subsection (); + + } + prm.leave_subsection (); + } + prm.leave_subsection (); + } + + + + template + void + DynamicCore::parse_parameters (ParameterHandler &prm) + { + prm.enter_subsection("Boundary temperature model"); + { + prm.enter_subsection("Dynamic core"); + { + // verify that the geometry is in fact a spherical shell since only + // for this geometry do we know for sure what boundary indicators it + // uses and what they mean + AssertThrow (Plugins::plugin_type_matches>(this->get_geometry_model()), + ExcMessage ("This boundary model is only implemented if the geometry is " + "a spherical shell.")); + + inner_temperature = prm.get_double ("Inner temperature"); + outer_temperature = prm.get_double ("Outer temperature"); + init_dT_dt = prm.get_double ("dT over dt") / year_in_seconds; + init_dR_dt = prm.get_double ("dR over dt") / year_in_seconds * 1.e3; + init_dX_dt = prm.get_double ("dX over dt") / year_in_seconds; + Rho_cen = prm.get_double ("Core density"); + g = prm.get_double ("Gravity acceleration"); + P_CMB = prm.get_double ("CMB pressure"); + X_init = prm.get_double ("Initial light composition"); + max_steps = prm.get_integer ("Max iteration"); + Cp = prm.get_double ("Core heat capacity"); + CpRho = Cp*Rho_cen; + + //\cite{NPB+04} + K0 = prm.get_double ("K0"); + Alpha = prm.get_double ("Alpha"); + Rho_0 = prm.get_double ("Rho0"); + Lh = prm.get_double ("Lh"); + Beta_c = prm.get_double ("Beta composition"); + k_c = prm.get_double ("Core conductivity"); + Delta = prm.get_double ("Delta"); + Rh = prm.get_double ("Rh"); + + prm.enter_subsection("Geotherm parameters"); + { + Tm0 = prm.get_double ("Tm0"); + Tm1 = prm.get_double ("Tm1"); + Tm2 = prm.get_double ("Tm2"); + Theta = prm.get_double ("Theta"); + composition_dependency = prm.get_bool("Composition dependency"); + use_bw11 = prm.get_bool("Use BW11"); + } + prm.leave_subsection (); + + prm.enter_subsection("Radioactive heat source"); + { + n_radioheating_elements = prm.get_integer ("Number of radioactive heating elements"); + heating_rate = Utilities::string_to_double + (Utilities::split_string_list + (prm.get("Heating rates"))); + AssertThrow(n_radioheating_elements==heating_rate.size(), + ExcMessage("Number of heating rate entities does not match " + "the number of radioactive elements.")); + half_life = Utilities::string_to_double + (Utilities::split_string_list + (prm.get("Half life times"))); + AssertThrow(n_radioheating_elements==half_life.size(), + ExcMessage("Number of half life time entities does not match " + "the number of radioactive elements.")); + initial_concentration = Utilities::string_to_double + (Utilities::split_string_list + (prm.get("Initial concentrations"))); + AssertThrow(n_radioheating_elements==initial_concentration.size(), + ExcMessage("Number of initial concentration entities does not match " + "the number of radioactive elements.")); + } + prm.leave_subsection (); + + prm.enter_subsection("Other energy source"); + { + name_OES = prm.get("File name"); + } + prm.leave_subsection (); + + 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 (); + } + prm.leave_subsection (); + } } }