Particle property Solid Comp

data/particle/generator/ascii/particle_property_composition.txt

@@ -0,0 +1,7 @@
+# Columns: x_coordinate y_coordinate
+
+5 210
+5 200
+5 190
+5 100
+5 300

doc/modules/changes/20170515_schools

@@ -1,6 +1,8 @@
 New: A new particle property is added which indicates the presence of melt 
 greater than the melt transport threshold at the particles position. If melt 
 is not present a 0 is recorded. If melt is present a 1 is recorded. Only works
-if there is a compositonal field named "porosity".
+if there is a compositonal field named "porosity". A second new particle property
+is also added which indicates the peridotite depletion value. This can be used
+to track the petrological evolution of material.
 <br>
 (Joe Schools, Rene Gassmoeller, 2017/05/15)

include/aspect/particle/property/composition.h

@@ -0,0 +1,121 @@
+/*
+ Copyright (C) 2015 - 2017 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/>.
+ */
+
+#ifndef _aspect_particle_property_composition_h
+#define _aspect_particle_property_composition_h
+
+#include <aspect/particle/property/interface.h>
+#include <aspect/simulator_access.h>
+
+namespace aspect
+{
+  namespace Particle
+  {
+    namespace Property
+    {
+      /**
+       * Implementation of a plugin in which the particle
+       * property is defined by the compositional fields in
+       * the model. This can be used to track solid composition
+       * evolution over time.
+       *
+       * @ingroup ParticleProperties
+       */
+      template <int dim>
+      class Composition : public Interface<dim>, public ::aspect::SimulatorAccess<dim>
+      {
+        public:
+          /**
+           * Initialization function. This function is called once at the
+           * creation of every particle for every property to initialize its
+           * value.
+           *
+           * @param [in] position The current particle position.
+           * @param [in,out] particle_properties The properties of the particle
+           * that is initialized within the call of this function. The purpose
+           * of this function should be to extend this vector by a number of
+           * properties.
+           */
+          virtual
+          void
+          initialize_one_particle_property (const Point<dim> &position,
+                                            std::vector<double> &particle_properties) const;
+
+          /**
+           * Update function. This function is called every time an update is
+           * request by need_update() for every particle for every property.
+           *
+           * @param [in] data_position An unsigned integer that denotes which
+           * component of the particle property vector is associated with the
+           * current property. For properties that own several components it
+           * denotes the first component of this property, all other components
+           * fill consecutive entries in the @p particle_properties vector.
+           *
+           * @param [in] position The current particle position.
+           *
+           * @param [in] solution The values of the solution variables at the
+           * current particle position.
+           *
+           * @param [in] gradients The gradients of the solution variables at
+           * the current particle position.
+           *
+           * @param [in,out] particle_properties The properties of the particle
+           * that is updated within the call of this function.
+           */
+          virtual
+          void
+          update_one_particle_property (const unsigned int data_position,
+                                        const Point<dim> &position,
+                                        const Vector<double> &solution,
+                                        const std::vector<Tensor<1,dim> > &gradients,
+                                        const ArrayView<double> &particle_properties) const;
+
+          /**
+           * This implementation tells the particle manager that
+           * we need to update particle properties over time.
+           */
+          UpdateTimeFlags
+          need_update () const;
+
+          /**
+           * Return which data has to be provided to update the property.
+           * The pressure and temperature need the values of their variables.
+           */
+          virtual
+          UpdateFlags
+          get_needed_update_flags () const;
+
+          /**
+           * Set up the information about the names and number of components
+           * this property requires.
+           *
+           * @return A vector that contains pairs of the property names and the
+           * number of components this property plugin defines.
+           */
+          virtual
+          std::vector<std::pair<std::string, unsigned int> >
+          get_property_information() const;
+      };
+    }
+  }
+}
+
+#endif
+

source/particle/property/composition.cc

@@ -0,0 +1,112 @@
+/*
+  Copyright (C) 2015 - 2017 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/>.
+ */
+
+
+#include <aspect/particle/property/composition.h>
+
+namespace aspect
+{
+  namespace Particle
+  {
+    namespace Property
+    {
+      template <int dim>
+      void
+      Composition<dim>::initialize_one_particle_property(const Point<dim> &position,
+                                                         std::vector<double> &data) const
+      {
+        for (unsigned int i = 0; i < this->n_compositional_fields(); i++)
+          data.push_back(this->get_initial_composition_manager().initial_composition(position,i));
+      }
+
+      template <int dim>
+      void
+      Composition<dim>::update_one_particle_property(const unsigned int data_position,
+                                                     const Point<dim> &,
+                                                     const Vector<double> &solution,
+                                                     const std::vector<Tensor<1,dim> > &,
+                                                     const ArrayView<double> &data) const
+      {
+        for (unsigned int i = 0; i < this->n_compositional_fields(); i++)
+          {
+            const unsigned int solution_component = this->introspection().component_indices.compositional_fields[i];
+            data[data_position+i] = solution[solution_component];
+          }
+      }
+
+      template <int dim>
+      UpdateTimeFlags
+      Composition<dim>::need_update() const
+      {
+        return update_time_step;
+      }
+
+      template <int dim>
+      UpdateFlags
+      Composition<dim>::get_needed_update_flags () const
+      {
+        return update_values;
+      }
+
+      template <int dim>
+      std::vector<std::pair<std::string, unsigned int> >
+      Composition<dim>::get_property_information() const
+      {
+
+        AssertThrow(this->n_compositional_fields() > 0,
+                    ExcMessage("You have requested the particle property <composition>, "
+                               "but the number of compositional fields is 0. "
+                               "Please add compositional fields to your model, or remove "
+                               "this particle property."));
+
+        std::vector<std::pair<std::string,unsigned int> > property_information;
+
+
+
+        for (unsigned int i = 0; i < this->n_compositional_fields(); i++)
+          {
+            std::ostringstream field_name;
+            field_name << this->introspection().name_for_compositional_index(i);
+            property_information.push_back(std::make_pair(field_name.str(),1));
+          }
+        return property_information;
+      }
+
+    }
+  }
+}
+
+// explicit instantiations
+namespace aspect
+{
+  namespace Particle
+  {
+    namespace Property
+    {
+      ASPECT_REGISTER_PARTICLE_PROPERTY(Composition,
+                                        "composition",
+                                        "Implementation of a plugin in which the particle "
+                                        "property is defined by the compositional fields in "
+                                        "the model. This can be used to track solid composition"
+                                        "evolution over time.")
+    }
+  }
+}
+

tests/particle_property_composition.cc

@@ -0,0 +1 @@
+#include <../benchmarks/solitary_wave/solitary_wave.cc>