[cig-commits] [commit] master: Add radioactive decay plugin (8afeaa7)
cig_noreply at geodynamics.org
cig_noreply at geodynamics.org
Thu May 22 15:34:55 PDT 2014
Repository : https://github.com/geodynamics/aspect
On branch : master
Link : https://github.com/geodynamics/aspect/compare/4a5e27e5690a0e0c1e13e19100af7398e37f36b8...d31f4e9435a1c12781f5b673b672cbd29c41c2c2
>---------------------------------------------------------------
commit 8afeaa7dae74630bd764dc6d133d1987c3381ee1
Author: Siqi Zhang <siqi.zhang at mq.edu.au>
Date: Mon May 19 16:14:07 2014 -0500
Add radioactive decay plugin
Add radioactive decay plugin.
Add postprocess plugin to output change in average heating rate.
Add function in simulator_access to get the pointer of heating model.
>---------------------------------------------------------------
8afeaa7dae74630bd764dc6d133d1987c3381ee1
.../{function.h => radioactive_decay.h} | 16 ++-
...pography.h => radioactive_heating_statistics.h} | 11 +-
include/aspect/simulator_access.h | 5 +
source/heating_model/radioactive_decay.cc | 152 +++++++++++++++++++++
.../postprocess/radioactive_heating_statistics.cc | 140 +++++++++++++++++++
source/simulator/simulator_access.cc | 7 +
6 files changed, 319 insertions(+), 12 deletions(-)
diff --git a/include/aspect/heating_model/function.h b/include/aspect/heating_model/radioactive_decay.h
similarity index 79%
copy from include/aspect/heating_model/function.h
copy to include/aspect/heating_model/radioactive_decay.h
index f6679e8..443fe6c 100644
--- a/include/aspect/heating_model/function.h
+++ b/include/aspect/heating_model/radioactive_decay.h
@@ -20,8 +20,8 @@
/* $Id$ */
-#ifndef __aspect__heating_model_function_h
-#define __aspect__heating_model_function_h
+#ifndef __aspect__heating_model_radioactive_decay_h
+#define __aspect__heating_model_radioactive_decay_h
#include <aspect/heating_model/interface.h>
#include <aspect/simulator_access.h>
@@ -41,13 +41,13 @@ namespace aspect
* @ingroup HeatingModels
*/
template <int dim>
- class Function : public Interface<dim>, public SimulatorAccess<dim>
+ class Radioactive_decay : public Interface<dim>, public SimulatorAccess<dim>
{
public:
/**
* Constructor.
*/
- Function ();
+ Radioactive_decay ();
/**
* Return the specific heating rate as calculated by the
@@ -85,9 +85,13 @@ namespace aspect
private:
/**
- * A function object representing the components of the velocity.
+ * Parameters for decaying radioactive heating.
*/
- Functions::ParsedFunction<dim> heating_model_function;
+ double time;
+ unsigned int num_radio_heating_elements;
+ std::vector<double> half_decay_time;
+ std::vector<double> radioactive_heating_rate;
+ std::vector<double> radioactive_initial_consentration;
};
}
}
diff --git a/include/aspect/postprocess/dynamic_topography.h b/include/aspect/postprocess/radioactive_heating_statistics.h
similarity index 74%
copy from include/aspect/postprocess/dynamic_topography.h
copy to include/aspect/postprocess/radioactive_heating_statistics.h
index b14efd3..2c1496f 100644
--- a/include/aspect/postprocess/dynamic_topography.h
+++ b/include/aspect/postprocess/radioactive_heating_statistics.h
@@ -20,29 +20,28 @@
/* $Id$ */
-#ifndef __aspect__postprocess_surface_topography_h
-#define __aspect__postprocess_surface_topography_h
+#ifndef __aspect__postprocess_radioactive_heating_statistics_h
+#define __aspect__postprocess_radioactive_heating_statistics_h
#include <aspect/postprocess/interface.h>
#include <aspect/simulator_access.h>
-
namespace aspect
{
namespace Postprocess
{
/**
- * A postprocessor that computes dynamic topography at the surface.
+ * A postprocessor that computes some statistics about the radioactive heating.
*
* @ingroup Postprocessing
*/
template <int dim>
- class DynamicTopography : public Interface<dim>, public ::aspect::SimulatorAccess<dim>
+ class Radioactive_Heating_Statistics : public Interface<dim>, public ::aspect::SimulatorAccess<dim>
{
public:
/**
- * Evaluate the solution for the dynamic topography.
+ * Evaluate the solution for some radioactive heating statistics.
*/
virtual
std::pair<std::string,std::string>
diff --git a/include/aspect/simulator_access.h b/include/aspect/simulator_access.h
index 2635af0..a1f9a02 100644
--- a/include/aspect/simulator_access.h
+++ b/include/aspect/simulator_access.h
@@ -404,6 +404,11 @@ namespace aspect
const std::set<types::boundary_id> &
get_fixed_temperature_boundary_indicators () const;
+ /**
+ * Return a pointer to the heating model.
+ */
+ const HeatingModel::Interface<dim> &
+ get_heating_model () const;
/**
* A convenience function that copies the values of the compositional
diff --git a/source/heating_model/radioactive_decay.cc b/source/heating_model/radioactive_decay.cc
new file mode 100644
index 0000000..01a8d3d
--- /dev/null
+++ b/source/heating_model/radioactive_decay.cc
@@ -0,0 +1,152 @@
+/*
+ Copyright (C) 2014 by the authors of the ASPECT code.
+
+ This file is part of ASPECT.
+
+ ASPECT is free software; you can redistribute it and/or modify
+ it under the terms of the GNU General Public License as published by
+ the Free Software Foundation; either version 2, or (at your option)
+ any later version.
+
+ ASPECT is distributed in the hope that it will be useful,
+ but WITHOUT ANY WARRANTY; without even the implied warranty of
+ MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ GNU General Public License for more details.
+
+ You should have received a copy of the GNU General Public License
+ along with ASPECT; see the file doc/COPYING. If not see
+ <http://www.gnu.org/licenses/>.
+*/
+/* $Id$ */
+
+
+#include <aspect/heating_model/radioactive_decay.h>
+#include <aspect/global.h>
+
+namespace aspect
+{
+ namespace HeatingModel
+ {
+ template <int dim>
+ Radioactive_decay<dim>::Radioactive_decay ()
+ {}
+
+
+
+ template <int dim>
+ double
+ Radioactive_decay<dim>::
+ specific_heating_rate (const double,
+ const double,
+ const std::vector<double> &,
+ const Point<dim> &p) const
+ {
+ double timedependent_radioactive_heating_rate=0;
+ if(num_radio_heating_elements!=0)
+ for(unsigned i_radio=0;i_radio<num_radio_heating_elements;i_radio++)
+ timedependent_radioactive_heating_rate+=
+ radioactive_heating_rate[i_radio]
+ *radioactive_initial_consentration[i_radio]*1e-6
+ *std::pow(0.5,time/half_decay_time[i_radio]);
+ return (timedependent_radioactive_heating_rate);
+ }
+
+
+ template <int dim>
+ void
+ Radioactive_decay<dim>::update ()
+ {
+ time = this->get_time();
+ // we get time passed as seconds (always) but may want
+ // to reinterpret it in years
+ if (this->convert_output_to_years())
+ time/=year_in_seconds;
+ }
+
+
+ template <int dim>
+ void
+ Radioactive_decay<dim>::declare_parameters (ParameterHandler &prm)
+ {
+ prm.enter_subsection("Heating model");
+ {
+ prm.enter_subsection("Radioactive decay");
+ {
+ prm.declare_entry("Number of elements","0",
+ Patterns::Integer(0),
+ "Number of radioactive elements");
+ prm.declare_entry("Heating rate","",
+ Patterns::List (Patterns::Double ()),
+ "Heating rate of different element (W/kg)");
+ prm.declare_entry("Half decay time","",
+ Patterns::List (Patterns::Double (0)),
+ "Half decay time. Units: (Seconds), or (Years) if set 'use years instead of seconds'.");
+ prm.declare_entry("Initial consentration","",
+ Patterns::List (Patterns::Double (0)),
+ "Initial consentration of different elments (ppm)");
+ }
+ prm.leave_subsection();
+ }
+ prm.leave_subsection();
+ }
+
+
+ template <int dim>
+ void
+ Radioactive_decay<dim>::parse_parameters (ParameterHandler &prm)
+ {
+ prm.enter_subsection("Heating model");
+ {
+ prm.enter_subsection("Radioactive decay");
+ {
+
+ num_radio_heating_elements= prm.get_integer ("Number of elements");
+ radioactive_heating_rate=Utilities::string_to_double
+ (Utilities::split_string_list
+ (prm.get("Heating rate")));
+ AssertThrow(radioactive_heating_rate.size()==num_radio_heating_elements,
+ ExcMessage("Number of heating rate entities does not match "
+ "the number of radioactive elements."));
+ radioactive_initial_consentration=Utilities::string_to_double
+ (Utilities::split_string_list
+ (prm.get("Initial consentration")));
+ AssertThrow(radioactive_initial_consentration.size()==num_radio_heating_elements,
+ ExcMessage("Number of initial consentration entities does not match "
+ "the number of radioactive elements."));
+ half_decay_time=Utilities::string_to_double
+ (Utilities::split_string_list
+ (prm.get("Half decay time")));
+ AssertThrow(half_decay_time.size()==num_radio_heating_elements,
+ ExcMessage("Number of half decay time entities does not match "
+ "the number of radioactive elements."));
+ }
+ prm.leave_subsection();
+ }
+ prm.leave_subsection();
+ }
+
+ }
+}
+
+// explicit instantiations
+namespace aspect
+{
+ namespace HeatingModel
+ {
+ ASPECT_REGISTER_HEATING_MODEL(Radioactive_decay,
+ "radioactive decay",
+ "Implementation of a model in which the heating "
+ "rate is given in terms of an explicit formula "
+ "that is elaborated in the parameters in section "
+ "``Heating model|Function''. "
+ "\n\n"
+ "The formula is interpreted as having units "
+ "W/kg."
+ "\n\n"
+ "Since the symbol $t$ indicating time "
+ "may appear in the formulas for the heating rate"
+ ", it is interpreted as having units "
+ "seconds unless the global parameter ``Use "
+ "years in output instead of seconds'' is set.")
+ }
+}
diff --git a/source/postprocess/radioactive_heating_statistics.cc b/source/postprocess/radioactive_heating_statistics.cc
new file mode 100644
index 0000000..606c182
--- /dev/null
+++ b/source/postprocess/radioactive_heating_statistics.cc
@@ -0,0 +1,140 @@
+/*
+ Copyright (C) 2011 - 2014 by the authors of the ASPECT code.
+
+ This file is part of ASPECT.
+
+ ASPECT is free software; you can redistribute it and/or modify
+ it under the terms of the GNU General Public License as published by
+ the Free Software Foundation; either version 2, or (at your option)
+ any later version.
+
+ ASPECT is distributed in the hope that it will be useful,
+ but WITHOUT ANY WARRANTY; without even the implied warranty of
+ MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ GNU General Public License for more details.
+
+ You should have received a copy of the GNU General Public License
+ along with ASPECT; see the file doc/COPYING. If not see
+ <http://www.gnu.org/licenses/>.
+*/
+/* $Id$ */
+
+
+#include <aspect/postprocess/radioactive_heating_statistics.h>
+#include <aspect/simulator_access.h>
+
+#include <deal.II/base/quadrature_lib.h>
+#include <deal.II/fe/fe_values.h>
+
+#include <boost/archive/text_oarchive.hpp>
+#include <boost/archive/text_iarchive.hpp>
+
+
+namespace aspect
+{
+ namespace Postprocess
+ {
+ template <int dim>
+ std::pair<std::string,std::string>
+ Radioactive_Heating_Statistics<dim>::execute (TableHandler &statistics)
+ {
+ const HeatingModel &heating_model=this->get_heating_model();
+
+ // create a quadrature formula based on the compositional element alone.
+ // be defensive about determining that a compositional field actually exists
+ AssertThrow (this->introspection().base_elements.compositional_fields
+ != numbers::invalid_unsigned_int,
+ ExcMessage("This postprocessor cannot be used without compositional fields."));
+ const QGauss<dim> quadrature_formula (this->get_fe().base_element(this->introspection().base_elements.temperature).degree+1);
+ const unsigned int n_q_points = quadrature_formula.size();
+
+ FEValues<dim> fe_values (this->get_mapping(),
+ this->get_fe(),
+ quadrature_formula,
+ update_values |
+ update_quadrature_points |
+ update_JxW_values);
+
+ std::vector<double> temperature_values (n_q_points);
+ std::vector<double> pressure_values (n_q_points);
+ std::vector<std::vector<double> > compositional_values (this->n_compositional_fields(),std::vector<double> (n_q_points));
+ std::vector<double> composition_values_at_q_point (this->n_compositional_fields());
+ std::vector<Point<dim>> quadrature_points(n_q_points);
+
+ typename DoFHandler<dim>::active_cell_iterator
+ cell = this->get_dof_handler().begin_active(),
+ endc = this->get_dof_handler().end();
+
+ double local_radioactive_heating_integrals;
+ double local_volume;
+
+ // compute the integral quantities by quadrature
+ for (; cell!=endc; ++cell)
+ if (cell->is_locally_owned())
+ {
+ fe_values.reinit (cell);
+ for(unsigned int c=0;c<
+ fe_values[this->introspection().extractors.temperature].get_function_values (this->get_solution(),
+ temperature_values);
+ fe_values[this->introspection().extractors.pressure].get_function_values (this->get_solution(),
+ pressure_values);
+ for(unsigned c=0;c<this->n_compositional_fields();c++)
+ fe_values[this->introspection().extractors.compositional_fields[c]].get_function_values (this->get_solution(),
+ compositional_values);
+ quadrature_points=fe_values.get_quadrature_points();
+
+ for (unsigned int q=0; q<n_q_points; ++q)
+ {
+ for(unsigned c=0;this->n_compositional_fields();c++)
+ composition_values_at_q_point[c]=compositional_values[c][q];
+ local_radioactive_heating_integrals += (heating_model.specific_heating_rate(temperature_values[q],
+ pressure_values[q],
+ composition_values_at_q_point,
+ quadrature_points[q])*fe_values.JxW(q);
+ local_volume+=fe_values.JxW(q);
+ }
+ }
+ // compute the sum over all processors
+ double global_radioactive_heating_integrals;
+ double global_volume;
+ Utilities::MPI::sum (local_radioactive_heating_integrals,
+ this->get_mpi_communicator(),
+ global_radioactive_heating_integrals);
+ Utilities::MPI::sum (local_volume,
+ this->get_mpi_communicator(),
+ global_volume);
+
+ // finally produce something for the statistics file
+ const std:string name="Average radioactive heating rate (W/kg) ";
+ statistics.add_value (name, global_radioactive_heating_integrals/global_volume);
+
+ // also make sure that the other columns filled by the this object
+ // all show up with sufficient accuracy and in scientific notation
+
+ statistics.set_precision (name, 8);
+ statistics.set_scientific (name, true);
+
+ std::ostringstream output;
+ output.precision(4);
+ output << global_radioactive_heating_integrals/global_volume
+ output << " (W//kg) ";
+ }
+
+ return std::pair<std::string, std::string> ("Average radioactive heating rate: ",
+ output.str());
+ }
+ }
+}
+
+
+// explicit instantiations
+namespace aspect
+{
+ namespace Postprocess
+ {
+ ASPECT_REGISTER_POSTPROCESSOR(Radioactive_Heating_Statistics,
+ "radioactive heating statistics",
+ "A postprocessor that computes some statistics about "
+ "radioactive heating, averaged by volume. ")
+ }
+}
diff --git a/source/simulator/simulator_access.cc b/source/simulator/simulator_access.cc
index e0fec9b..0d15f1f 100644
--- a/source/simulator/simulator_access.cc
+++ b/source/simulator/simulator_access.cc
@@ -328,6 +328,13 @@ namespace aspect
}
template <int dim>
+ const HeatingModel::Interface<dim> &
+ SimulatorAccess<dim>::get_heating_model () const
+ {
+ return *simulator->heating_model.get();
+ }
+
+ template <int dim>
void
SimulatorAccess<dim>::get_composition_values_at_q_point (const std::vector<std::vector<double> > &composition_values,
const unsigned int q,
More information about the CIG-COMMITS
mailing list