Hello folks,
Hope you’re doing good! I found peculiar behavior of a variable used in two customized plugins. Let me break down.
- The purpose is to track solid mantle convection (including a less dense crust layer, more dense layer below crust and peridotite mantle), melt production and segregation, radioactive decay heating from U238, U235, Th232 and K40. So I use 5 fields, field porosity for melt, field peridotite for mantle, field crust for crust, field ibc for dense layer, a passive representative field U238 for all radioactive isotopes. For now I just loaded field porosity, peridotite and U238. The setup can be seen in ibc.prm which is based on melt_global.prm.
- In ibc.prm, heating models include adiabatic heating out, radioactive decay heating. Model adiabatic heating out is derived from adiabatic_heating.cc. The reason to create model adiabatic heating out is to test if the peculiar variable heating_model_outputs.heating_source_terms is properly allocated and initialized. In plugin adiabatic_heating_out.cc, a couple of test lines can be seen. The address of heating_model_outputs.heating_source_terms [0] will be printed. I also noticed that heating_model_outputs.heating_source_terms[q] is equal to some contributions, not +=, since it’s the first heat source term:
if (!simplified_adiabatic_heating)
heating_model_outputs.heating_source_terms[q] = (material_model_inputs.velocity[q] * material_model_inputs.pressure_gradient[q])
* material_model_outputs.thermal_expansion_coefficients[q]
* material_model_inputs.temperature[q];
else
heating_model_outputs.heating_source_terms[q] = (material_model_inputs.velocity[q]
* this->get_gravity_model().gravity_vector(material_model_inputs.position[q]))
* material_model_outputs.thermal_expansion_coefficients[q]
* material_model_inputs.temperature[q]
* material_model_outputs.densities[q];
- Now comes to 2nd customized plugin any_decay.cc which is derived from radioactive_decay.cc.In evaluate function of any_decay.cc, the heat produced by radioactive decay is added to heat source term:
heating_model_outputs.heating_source_terms[q] += decay_heating.
The peculiar behavior happens here: adding a number will trigger floating point exception (FPE). FPE happens at first element of heating_model_outputs.heating_source_terms. I can confirm that the address of heating_model_outputs.heating_source_terms[0] is exactly the same as heating_model_outputs.heating_source_terms[0] in adiabatic_heating_out.cc. And the size of heating_model_outputs.heating_source_terms in both plugins are the same. Th FPE happens even you just read number from heating_model_outputs.heating_source_terms[0]:
double val = heating_model_outputs.heating_source_terms[q];
std::cout << " Value = " << val << std::endl;
The code is 0x7ff4000000000000, a signaling NaN. Why is heating_model_outputs.heating_source_terms[0] changed and where? Is there any other manipulation between adiabatic_heating_out.cc and any_decay.cc? This is wired. Any suggestions?
Thanks,
Mingming
P.S. (1) adiabatic_heating_out.cc
#include </home/dealii/aspect/model_input/src/adiabatic_heating_out.h>
#include <aspect/heating_model/adiabatic_heating.h>
#include <aspect/gravity_model/interface.h>
namespace aspect
{
namespace HeatingModel
{
template <int dim>
bool
AdiabaticHeatingOut<dim>::use_simplified_adiabatic_heating() const
{
return simplified_adiabatic_heating;
}
template <int dim>
void
AdiabaticHeatingOut<dim>::
evaluate (const MaterialModel::MaterialModelInputs<dim> &material_model_inputs,
const MaterialModel::MaterialModelOutputs<dim> &material_model_outputs,
HeatingModel::HeatingModelOutputs &heating_model_outputs) const
{
for (unsigned int q = 0; q < heating_model_outputs.heating_source_terms.size(); ++q)
{heating_model_outputs.heating_source_terms[q] = 10.0;
std::cout<<"begining:heating_model_outputs.heating_source_terms[q] "<<heating_model_outputs.heating_source_terms[q]<<std::endl;}
const double forced_value = 10.0;
for (unsigned int q = 0; q < heating_model_outputs.heating_source_terms.size(); ++q)
{
std::memcpy(&heating_model_outputs.heating_source_terms[q], &forced_value, sizeof(double));
std::cout<<"forced:heating_model_outputs.heating_source_terms[q] "<<heating_model_outputs.heating_source_terms[q]<<std::endl;std::cout<<"heating_model_outputs.heating_source_terms[q] "<<heating_model_outputs.heating_source_terms[q]<<std::endl;
}
Assert(heating_model_outputs.heating_source_terms.size() == material_model_inputs.n_evaluation_points(),
ExcMessage ("Heating outputs need to have the same number of entries as the material model inputs."));
std::cout<<"size of heating_model_outputs.heating_source_terms "<<heating_model_outputs.heating_source_terms.size()<<std::endl;
for (unsigned int q=0; q<heating_model_outputs.heating_source_terms.size(); ++q)
{
if (!simplified_adiabatic_heating)
heating_model_outputs.heating_source_terms[q] = (material_model_inputs.velocity[q] * material_model_inputs.pressure_gradient[q])
* material_model_outputs.thermal_expansion_coefficients[q]
* material_model_inputs.temperature[q];
else
heating_model_outputs.heating_source_terms[q] = (material_model_inputs.velocity[q]
* this->get_gravity_model().gravity_vector(material_model_inputs.position[q]))
* material_model_outputs.thermal_expansion_coefficients[q]
* material_model_inputs.temperature[q]
* material_model_outputs.densities[q];
heating_model_outputs.lhs_latent_heat_terms[q] = 0.0;
std::cout<<"This is adiabatic heating out"<<std::endl;
std::cout<<"heating_model_outputs.heating_source_terms[q] "<<heating_model_outputs.heating_source_terms[q]<<std::endl;
std::cout<<"heating_model_outputs.lhs_latent_heat_terms[q] "<<heating_model_outputs.lhs_latent_heat_terms[q]<<std::endl;
Assert(q < heating_model_outputs.heating_source_terms.size(), ExcMessage("Index out of bounds!"));
}
std::cout << "Address of heating_source_terms in adiabatic plugin: "
<< static_cast<const void *>(heating_model_outputs.heating_source_terms.data())
<< std::endl;
}
template <int dim>
void
AdiabaticHeatingOut<dim>::declare_parameters (ParameterHandler &prm)
{
prm.enter_subsection("Heating model");
{
prm.enter_subsection("Adiabatic heating out");
{
prm.declare_entry ("Use simplified adiabatic heating", "false",
Patterns::Bool (),
"A flag indicating whether the adiabatic heating should be simplified "
"from $\\alpha T (\\mathbf u \\cdot \\nabla p)$ to "
"$ \\alpha \\rho T (\\mathbf u \\cdot \\mathbf g) $.");
}
prm.leave_subsection();
}
prm.leave_subsection();
}
template <int dim>
void
AdiabaticHeatingOut<dim>::parse_parameters (ParameterHandler &prm)
{
prm.enter_subsection("Heating model");
{
prm.enter_subsection("Adiabatic heating out");
{
simplified_adiabatic_heating = prm.get_bool ("Use simplified adiabatic heating");
}
prm.leave_subsection();
}
prm.leave_subsection();
}
}
}
// explicit instantiations
namespace aspect
{
namespace HeatingModel
{
ASPECT_REGISTER_HEATING_MODEL(AdiabaticHeatingOut,
"adiabatic heating out",
"Implementation of a standard and a simplified model of "
"adiabatic heating.")
}
}
(2)adiabatic_heating_out.h
#ifndef _aspect_heating_model_adiabatic_heating_h
#define _aspect_heating_model_adiabatic_heating_h
#include <aspect/heating_model/interface.h>
#include <aspect/simulator_access.h>
namespace aspect
{
namespace HeatingModel
{
using namespace dealii;
/*
* @ingroup HeatingModels
*/
template <int dim>
class AdiabaticHeatingOut : public Interface<dim>, public ::aspect::SimulatorAccess<dim>
{
public:
bool
use_simplified_adiabatic_heating() const;
void
evaluate (const MaterialModel::MaterialModelInputs<dim> &material_model_inputs,
const MaterialModel::MaterialModelOutputs<dim> &material_model_outputs,
HeatingModel::HeatingModelOutputs &heating_model_outputs) const override;
/**
* @name Functions used in dealing with run-time parameters
* @{
*/
static
void
declare_parameters (ParameterHandler &prm);
void
parse_parameters (ParameterHandler &prm) override;
/**
* @}
*/
private:
bool simplified_adiabatic_heating;
};
}
}
#endif
(3)any_decay.cc
#include </home/dealii/aspect/model_input/src/any_decay.h>
#include <aspect/heating_model/interface.h>
#include <aspect/simulator_access.h>
#include <deal.II/base/parameter_handler.h>
#include <deal.II/base/function_lib.h>
#include <deal.II/base/point.h>
#include <aspect/gravity_model/interface.h>
//To include U, Th, K radioactive decay
//Assuming radioactive elements will always stay with ilmenite and will not mix with mantle or differentiate or relative differentiate wrt each other during partial melting
namespace aspect
{
namespace HeatingModel
{
template <int dim>
void AnyDecayHeating<dim>::evaluate(const MaterialModel::MaterialModelInputs<dim> &in,
const MaterialModel::MaterialModelOutputs<dim> &material_model_outputs,
HeatingModel::HeatingModelOutputs &heating_model_outputs) const
{
std::cout<<"size of heating_source_terms "<<heating_model_outputs.heating_source_terms.size()<<std::endl;
for (unsigned int q = 0; q < heating_model_outputs.heating_source_terms.size(); ++q)
{
uint64_t bits;
std::memcpy(&bits, &heating_model_outputs.heating_source_terms[q], sizeof(double));
// Check if exponent bits are all 1s (NaN or Inf)
if ((bits & 0x7ff0000000000000) == 0x7ff0000000000000)
{
std::cerr << "Invalid value at q=" << q << ", bits = 0x"
<< std::hex << bits << std::dec << std::endl;
std::abort();
}
std::cout << " Entry " << q << " address = "
<< static_cast<const void *>(&heating_model_outputs.heating_source_terms[q]) << std::endl;
double val = heating_model_outputs.heating_source_terms[q]; // FPE here if still corrupted
std::cout << " Value = " << val << std::endl;
}
// Get current time
const double current_time = this->get_time();
const double current_dt = this->get_timestep();
const unsigned int U238_idx = this->introspection().compositional_index_for_name("U238");
std::cout<<"U238_idx "<<U238_idx<<std::endl;
std::cout<<"initial_U238 "<<initial_U238<<std::endl;
std::cout<<"decay_const_U238 "<<decay_const_U238<<std::endl;
std::cout<<"current_time "<<current_time<<std::endl;
std::cout<<"current_dt "<<current_dt<<std::endl;
std::cout<<"points "<<in.position.size()<<std::endl;
std::cout<<"size of heating_source_terms "<<heating_model_outputs.heating_source_terms.size()<<std::endl;
std::cout<<"velocity "<<in.velocity[0]<<std::endl;
std::cout<<"thermal_expansion_coefficients "<<material_model_outputs.thermal_expansion_coefficients[0]<<std::endl;
std::cout<<"temperature "<<in.temperature[0]<<std::endl;
std::cout<<"pressure_gradient "<<in.pressure_gradient[0]<<std::endl;
std::cout<<"densities "<<material_model_outputs.densities[0]<<std::endl;
std::cout<<"this->get_gravity_model().gravity_vector(material_model_inputs.position[q]) "<<this->get_gravity_model().gravity_vector(in.position[0])<<std::endl;
for (unsigned int q=0; q<heating_model_outputs.heating_source_terms.size(); ++q)
{double aa = (in.velocity[q] * in.pressure_gradient[q])* material_model_outputs.thermal_expansion_coefficients[q] * in.temperature[q];
std::cout<<"heating_model_outputs.heating_source_terms[q] "<<aa<<std::endl;
}
std::cout << "Address of heating_source_terms in radioactive decay heating: "
<< static_cast<const void *>(heating_model_outputs.heating_source_terms.data())
<< std::endl;
for (unsigned int p = 0; p < heating_model_outputs.heating_source_terms.size(); ++p)
{
std::cout << "Reading entry " << p << std::endl;
double val = heating_model_outputs.heating_source_terms[p]; // Possibly triggers FPE
std::cout << " Address = " << static_cast<const void *>(&heating_model_outputs.heating_source_terms[p]) << std::endl;
std::cout << " Value = " << val << std::endl;
}
// Loop over all quadrature points
for (unsigned int q = 0; q < in.position.size(); ++q)
{
// Retrieve the concentrations of U238, U235, Th232 and K40 from compositional fields
//composition[i][0] is porosity, composition[i][1] is peridotite, composition[i][2] is crust, composition[i][3] ibc and composition[i][4] isrepresentative radioactive element
//composition[q][U238_idx] is a scalar field initiallized in prm.
const double concentration_U238 = in.composition[q][U238_idx] * initial_U238 * std::exp(-decay_const_U238 * current_time);
const double concentration_U235 = in.composition[q][U238_idx] * initial_U235 * std::exp(-decay_const_U235 * current_time);
const double concentration_Th232 = in.composition[q][U238_idx] * initial_Th232 * std::exp(-decay_const_Th232 * current_time);
const double concentration_K40 = in.composition[q][U238_idx] * initial_K40 * std::exp(-decay_const_K40 * current_time);
std::cout<<"place/point "<<q<<std::endl;
std::cout<<"bulk density "<<material_model_outputs.densities[q]<<std::endl;
// Compute the heat production from radioactive decay
// I need bulk density here
double decay_heating = (decay_const_U238 * concentration_U238 * q_U238 * current_dt * material_model_outputs.densities[q]) +
(decay_const_U235 * concentration_U235 * q_U235 * current_dt * material_model_outputs.densities[q]) +
(decay_const_Th232 * concentration_Th232 * q_Th232 * current_dt * material_model_outputs.densities[q]) +
(decay_const_K40 * concentration_K40 * q_K40 * current_dt * material_model_outputs.densities[q]);
std::cout<<"decay heating "<<decay_heating<<std::endl;
// Add the radioactive heating to the heating outputs
//decay_heating = 0.1;
//std::cout<<"decay heating "<<decay_heating<<std::endl;
heating_model_outputs.heating_source_terms[q] += decay_heating;
}
if (heating_model_outputs.lhs_latent_heat_terms.size() != 0)
{
if (Latent_heat_contribution)
{
std::cout << "[AnyDecayHeating] Note: phase change is expected to be handled by another plugin." << std::endl;
// Don't touch lhs_latent_heat_terms here
}
else
{
// If no other plugin provides phase change, we must zero it to avoid crash
for (unsigned int q = 0; q < heating_model_outputs.lhs_latent_heat_terms.size(); ++q)
heating_model_outputs.lhs_latent_heat_terms[q] = 0.0;
}
}
std::cout<<"Sweep once"<<std::endl;
}
template <int dim>
void AnyDecayHeating<dim>::declare_parameters(ParameterHandler &prm)
{
prm.enter_subsection("Heating model");
{
prm.enter_subsection("Radioactive decay heating");
{
//Default values are calculated based on present day concentrations in urKREEP,
//decay constants are listed in private members,
//present day concentrations are from Korotev2000 and Wieczorek2000
prm.declare_entry("Initial concentration of U238", "6.85e-6",
Patterns::Double(0.),
"Initial concentration of U238 (kg/kg)");
prm.declare_entry("Initial concentration of U235", "2.15e-6",
Patterns::Double(0.),
"Initial concentration of U235 (kg/kg)");
prm.declare_entry("Initial concentration of Th232", "1.54e-5",
Patterns::Double(0.),
"Initial concentration of Th232 (kg/kg)");
prm.declare_entry("Initial concentration of K40", "9.15e-6",
Patterns::Double(0.),
"Initial concentration of K40 (kg/kg)");
prm.declare_entry("Latent heat contribution", "false",
Patterns::Bool(),
"Latent heat contribution");
}
prm.leave_subsection();
}
prm.leave_subsection();
}
template <int dim>
void AnyDecayHeating<dim>::parse_parameters(ParameterHandler &prm)
{
prm.enter_subsection("Heating model");
{
prm.enter_subsection("Radioactive decay heating");
{
initial_U238 = prm.get_double("Initial concentration of U238");
initial_U235 = prm.get_double("Initial concentration of U235");
initial_Th232 = prm.get_double("Initial concentration of Th232");
initial_K40 = prm.get_double("Initial concentration of K40");
Latent_heat_contribution = prm.get_bool("Latent heat contribution");
}
prm.leave_subsection();
}
prm.leave_subsection();
}
}
}
namespace aspect
{
namespace HeatingModel
{
ASPECT_REGISTER_HEATING_MODEL(AnyDecayHeating,
"radioactive decay heating",
"A heating model that computes radioactive decay heating "
"from any number of radioactive elements based on one"
"representative element presented by one compositional field.")
}
}
(4)any_decay.h
#ifndef ASPECT_ANY_DECAY_H
#define ASPECT_ANY_DECAY_H
#include <aspect/heating_model/interface.h>
#include <aspect/simulator_access.h>
namespace aspect
{
namespace HeatingModel
{
using namespace dealii;
/*
* A heating model plugin that implements heating from radioactive decay.
* This model can read compositional fields such as U238 and compute the
* associated decay heating using user-defined initial concentrations.
* @ingroup HeatingModels
*/
template <int dim>
class AnyDecayHeating : public Interface<dim>, public ::aspect::SimulatorAccess<dim>
{
public:
void evaluate(const MaterialModel::MaterialModelInputs<dim> &material_model_inputs,
const MaterialModel::MaterialModelOutputs<dim> &material_model_outputs,
HeatingModel::HeatingModelOutputs &heating_model_outputs
) const override;
static void declare_parameters(ParameterHandler &prm);
virtual void parse_parameters(ParameterHandler &prm) override;
private:
double yr = 3600*24*365;//1 year in seconds
double decay_const_U238 = std::log(2.0)/(4.47e9*yr);//tau_U238 = 4.47e9 [year], time must be in seconds --> [1/s]
double decay_const_U235 = std::log(2.0)/(7.04e8*yr);//tau_U235 = 7.04e8 [year], time must be in seconds --> [1/s]
double decay_const_Th232 = std::log(2.0)/(1.4e10*yr);//tau_Th232 = 1.4e10 [year], time must be in seconds --> [1/s]
double decay_const_K40 = std::log(2.0)/(1.25e9*yr);//tau_K40 = 1.25e9 [year], time must be in seconds --> [1/s]
double initial_U238;
double initial_U235;
double initial_Th232;
double initial_K40;
double q_U238 = 9.46e-5;// [W/kg]
double q_U235 = 5.69e-4;//[W/kg]
double q_Th232 = 2.64e-5;//[W/kg]
double q_K40 = 2.92e-5;//[W/kg]
bool Latent_heat_contribution;
};
}
}
#endif
ibc.prm (18.3 KB)
