Adding a particle property for water content

Hi. I am trying to write a plugin in particle property that will contain water content in the particle. In the viscoplastic material model I will initialize the water content of each composition. For initialization of water content in each particle it will check the volume fractions of the compositions in it and their initial water content. In subsequent time steps, it will:
a) Ensure that it is not first time step,
b) Check the water content from previous time step for each composition in that particle,
c) Check the pressure and temperature information of the particle in current time step,
d) Input text files for each composition where the maximum water content allowed for that composition is stated for given pressure temperature and,
e) Check whether previous time step value is greater or less than this and update the current water content for that composition in the particle accordingly, then finally update the water content of the particle based on volume fractions of the compositions present.
f) After this I intend to use this water content of the particle for calculation of viscosity and density by including a water content dependence.
I understand accessing the time step as this->get_timestep().
This is how I intend to make the update work.

std::vector initial_water_capacity; // to include water content of compositions from previous time step
std::vector water_capacity; // to update water content of compositions
double pressure = solution[this->introspection().component_indices.pressure];
double temperature = solution[this->introspection().component_indices.temperature]
for (unsigned int j = 0; j < this->n_compositional_fields(); j++)
{ std::vector vecT, vecP, vecw;
double Tm, Pr, wa;

          std::ifstream inputFile("composition_j_water.txt");

          while(inputFile >> Tm >> Pr >> wa)
          {
              vecT.push_back(Tm);
              vecP.push_back(Pr);
              vecw.push_back(wa);
          }
          double Tmax = *max_element(vecT.begin(), vecT.end()); 
          double Pmax = *max_element(vecP.begin(), vecP.end());
          double distance[vecT.size()];

          for(int i(0); i < vecT.size(); i++)
          {
              distance[i] = sqrt(pow((vecT[i]-temperature)/Tmax,2.0)+pow((vecP[i]-pressure)/Pmax,2.0));
          } 

          double dist_min(distance[0]);
          for(int i(0); i < vecT.size(); i++)
          {    
              if(distance[i] < dist_min)
                  dist_min = distance[i];
          }

          double wnew;
          for(int i(0); i < vecT.size(); i++)
          {
              if(distance[i] == dist_min)
              {
                  if(initial_water_capacity[j] >= vecw[i])
                      wnew = vecw[i];
                  else
                      wnew = initial_water_capacity[j];
              }
          }

          water_capacity[j]= wnew;

}

I don’t understand how I can access previous time step information for the particles, if taking text file inputs is going to work, and how I can include a water content dependence in viscoplastic plugin for viscosity and density for the particles from this plugin.
I am new to ASPECT and would appreciate the help very much.

Regards,
Mainak

Mainak:
I’ve gone through your question again, but I must admit that I get lost in it (and suspect that others did as well) because there was so much detail to it. You explained the big picture and showed a bunch of code, but in the end, it is not clear to me what the actual question is. Can you cook it down a one concrete question at a time? For example, I think that some of your question are of the following kind:

  • I am writing a function that updates the value of a particle property. In it, I need access to the time. How do I get to that information? (Then maybe show the outline of the function, without the actual details of it as the details are not important to understand the question.)
  • I am trying to update the value of a particle property and that involves reading in a file. How do I do that?

In the end, I believe, but don’t actually know, that the details of what specifically you want to do (namely, have particles carry a water content property) are not relevant to what it is you want to ask, but these details make it quite difficult to understand what the actual question is.

Best
W.

Sorry for my vague broad description. I would say two of my main questions are:

  1. Can we take a table as input from a text file for particles in a new property that I would create in particle/property/
  2. As pT_path works for a particle so that we can access the pressure temperature information at each time-step of a particle given its id. We can plot this in output. Similarly I want to create this new particle-property where we can access information about this property for a particle at all time-steps. I need to access the previous time-step’s information of this particle-property at current time-step inside the code for this property in order to update the current time-step’s value of this property for a particle given its id. Rather than pushing this information to output or postprocess, can we access it from inside the code?

Thanks,
Mainak

Mainak:
To your questions:

  • Yes, you can take a table and initialize particle properties that way. I’m just unclear why you want to do that. In most cases, you want to initialize particle properties in a way that depends on where a particle is, or what the conditions at that location are. That’s, for example, what we do in the Particle::Property::Position and Particle::Property::PTPath plugins. My question for you would then be what it is you are storing in the file, and how that will be used for initializing particle properties? (Of course, you can always read in a file in the constructor of the initialize() function of a particle property class, store the information from the file in some member variable, and then use that to initialize an individual particle’s properties.) As a general rule, it’s often useful to read through the existing files in the source/particle/property/ directory and see what’s already there.

  • As for accessing the previous time step’s property values: That’s what the update_particle_property() function is there for. Take a look at the intergrated_strain.cc file, for example.

Best
W.

Thank you for the directions. Now, a particle has its composition tracked when “composition” is selected in list of particle properties. Let’s say there are 3 compositions in my model, and each composition has its corresponding “water content” which would be tracked by this new plugin. In my understanding, via particle interpolator, all particle property values of the particles in a cell are averaged and given to material model as input for calculation of, say, viscosity at the evaluation points of cells. The information of these 3 compositions is taken from the cell averaged particles and given as input to the material model in this way. Now I wish to supply the information about this “water content” property to material model. It would contain information about the water content within 3 compositions and would have the same dimensions as composition information. Can this new plugin information be supplied to material model in this way or in some other way, please let me know.

Thanks,
Mainak

And for the initialization, yes the values would be selected based on the compositions present in that position, as in each composition would have an initial water content prescribed.

Mainak,
yes, I think you have the idea correct. My recommendation would be to start with a simple model problem with just a single composition and an otherwise trivial setup where you try out what you describe and get a feel for how this is supposed to work.
Best
W.

I did these processes in the particle domain and then accessed their information in the cell domain for material property calculations. It is working as I had hoped for, thanks again. Finally, in the cell domain, I would like to access the information of adjacent cells vertically above and below the current cell and pass information to them from the current cell. Is there a way these cells’ information can be accessed from the material model plugin, say within the evaluate function call, please let me know.

Thanks,

Mainak

Mainak: If you have a cell, you can ask cell->neighbor(f) where f is the number of a face and the call returns the neighbor behind the face. Of course, you should think of the mesh as, in general, unstructured, and so notions such as “vertically above/below” is something you should define. In particular, “vertical” implies a direction that is worth thinking about.
Best
W.

I created a lookup table of temperature and pressure inside this plugin. The pressure and temperature at the particle position are then matched from this table. It is running ok for some time but suddenly gives this error

An error occurred in line <1684> of file </home/mainakchakraborty/bin/tmp/unpack/deal.II-v9.3.3/source/lac/trilinos_sparse_matrix.cc> in function
void dealii::TrilinosWrappers::SparseMatrix::add(dealii::TrilinosWrappers::SparseMatrix::size_type, dealii::TrilinosWrappers::SparseMatrix::size_type, const size_type*, const TrilinosScalar*, bool, bool)
The violated condition was:
ierr == 0
Additional information:
An error with error number 2 occurred while calling a Trilinos
function

Mainak – there is really not any way we can know what the problem is. In essence, you’re question is “My car doesn’t start, and there is a red light” kind of question :slight_smile:

You’ll have to provide more information (what you have done, whether or not you use any plugins and if so what their source code is, what the input file is), etc. In an idea world, you would run your program in a debugger and find out where the error originates, showing the whole backtrace, so we can see where the function that generates the error is being called.

Best
W.

My bad…

In this new plugin a lookup table of this format is taken, the columns representing temperature, pressure, and maximum-water-content respectively:

773 30e9 12
1023 30e9 12
1273 30e9 8
1373 30e9 0.1
1573 30e9 0.1
1773 30e9 0.1
773 40e9 12
1023 40e9 12
1273 40e9 8
1373 40e9 0.1
1573 40e9 0.1
1773 40e9 0.1

This is how the code for this looks like:

#include "particle_composition_water.h"
#include <aspect/initial_composition/interface.h>
#include <aspect/adiabatic_conditions/interface.h>

namespace aspect
{
  namespace Particle
  {
    namespace Property
    {
      template <int dim>
      void
      CompositionWater<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(6);
     }

      template <int dim>
      void
      CompositionWater<dim>::update_one_particle_property(const unsigned int data_position,
                                                     const Point<dim> &position,
                                                     const Vector<double> &solution,
                                                     const std::vector<Tensor<1,dim> > &,
						     const ArrayView<double> &data) const
      {
        const double tp = solution[this->introspection().component_indices.temperature];
        const double ps = solution[this->introspection().component_indices.pressure];
        const double ct = (this->get_time())/year_in_seconds;
	    const double depth = this->get_geometry_model().depth(position);

        double Tmax = (1773>tp ? 1773:tp); 
	double Pmax = (40e9>ps ? 40e9:ps);

	const std::string filename = data_directory+data_filename;
	std::istringstream lookup_ta(Utilities::read_and_distribute_file_content(filename, this->get_mpi_communicator()));

             for (unsigned int i = 0; i < this->n_compositional_fields(); i++)
             { 
                if (ct == 0e6)
                   data[data_position+i]=6;
               else
               {
		  double wnew = 0.0;
		  double dist1, dist2, distm;
		  double Tm1, Pr1, wa1;
		  int count = 0;
		while (lookup_ta >> Tm1 >> Pr1 >> wa1)
		{
		  if (count==0)
		  {
	    	    dist1 = sqrt(pow((Tm1-tp)/Tmax,2.0)+pow((Pr1-ps)/Pmax,2.0));
	    	    distm = dist1;
		    if(data[data_position+i] >= wa1)
		       wnew = wa1;
		    else
		       wnew = data[data_position+i];
		  }
		  else
		  {
		    dist2 = sqrt(pow((Tm1-tp)/Tmax,2.0)+pow((Pr1-ps)/Pmax,2.0));
		    if (dist2<distm)
		    {
			distm = dist2;
			if(data[data_position+i] >= wa1)
			   wnew = wa1;
			else
			   wnew = data[data_position+i];
		    }
		  }
		  count++;
		}
		  data[data_position+i]= wnew;
	this->get_pcout() << " \n " << depth << " " << position[0] << " " << data[data_position+i];
             }
            }
	lookup_ta.clear();
      }

      template <int dim>
      UpdateTimeFlags
      CompositionWater<dim>::need_update() const
      {
        return update_time_step;
      }

      template <int dim>
      UpdateFlags
      CompositionWater<dim>::get_needed_update_flags () const
      {
        return update_values;
      }

      template <int dim>
      std::vector<std::pair<std::string, unsigned int> >
      CompositionWater<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++)
          {
            const std::string field_name = this->introspection().name_for_compositional_index(i);
            property_information.emplace_back(field_name,1);
          }
        return property_information;
      }

      template <int dim>
      void
      CompositionWater<dim>::declare_parameters (ParameterHandler &prm)
      {
        prm.enter_subsection("Postprocess");
        {
          prm.enter_subsection("Particles");
            {
              prm.enter_subsection("Ascii file");
              {
                prm.declare_entry ("Data directory",
                                   "$ASPECT_SOURCE_DIR/data/particle/generator/ascii/",
                                   Patterns::DirectoryName (),
                                   "The name of a directory that contains the particle data. This path "
                                   "may either be absolute (if starting with a '/') or relative to "
                                   "the current directory. The path may also include the special "
                                   "text '$ASPECT_SOURCE_DIR' which will be interpreted as the path "
                                   "in which the ASPECT source files were located when ASPECT was "
                                   "compiled. This interpretation allows, for example, to reference "
                                   "files located in the `data/' subdirectory of ASPECT. ");
                prm.declare_entry ("Data file name", "particle.dat",
                                   Patterns::Anything (),
                                   "The name of the particle file.");
                prm.leave_subsection();
              }
              prm.leave_subsection();
          }
          prm.leave_subsection();
        }
      }


      template <int dim>
      void
      CompositionWater<dim>::parse_parameters (ParameterHandler &prm)
      {
        prm.enter_subsection("Postprocess");
        {
          prm.enter_subsection("Particles");
            {
              prm.enter_subsection("Ascii file");
              {
                data_directory = Utilities::expand_ASPECT_SOURCE_DIR(prm.get ("Data directory"));

                data_filename    = prm.get ("Data file name");
              }
              prm.leave_subsection();
            }
            prm.leave_subsection();
        }
        prm.leave_subsection();
      }
    }
  }
}

namespace aspect
{
  namespace Particle
  {
    namespace Property
    {
      ASPECT_REGISTER_PARTICLE_PROPERTY(CompositionWater,
                                        "composition water",
                                        "Implementation of a plugin in which the particle "
                                        "water content is defined by lookup table contents ")
    }
  }
}

For the compositions an initial water component is declared. In subsequent time-steps it is updated based on the condition- if the value given in the table is less, then the value at data[data_position+i] is updated to the value given in table. Inside the update particle property section with this->get_pcout() the outputs are printed. They are coming like this- depth, x-extent and the updated water content at that location:

1.46991e+06 2.81657e+06 0.1
1.47803e+06 2.84295e+06 0.1
1.45423e+06 2.82213e+06 0.1
1.4527e+06 2.84307e+06 0.1
1.47279e+06 2.84958e+06 0.1
1.46632e+06 2.84872e+06 0.1
1.45124e+06 2.82424e+06 0.1
1.47476e+06 2.84011e+06 0.1
1.45857e+06 2.8246e+06 0.1
1.47408e+06 2.84109e+06 0.1

There are many more output lines.

This code is running up to this point. As the depth in my model is given 2.9e6, it is running for almost half of my model then it stops suddenly. After quite some time it was giving the error mentioned here some times.

I should add… it is working for much smaller and simpler models that I ran. But for larger scales this stops running. Currently I am running aspect in my local installation. When I get some cluster access I would check if it can handle larger geophysical models, or there is a certain scope limitation for trilinos there also.

Mainak

Another peculiar observation… all these models run only up to 5 time-steps, regardless of time-steppings. Meaning, if I set the time-steps as 1000 years it runs up to 5000 years and if the same was 10,000 years, it runs up to 50,000 years. That is because that is set for mesh refinement time, with

subsection Mesh refinement
set Initial global refinement = 3
set Initial adaptive refinement = 2
set Time steps between mesh refinement = 5
set Refinement fraction = 0.9
set Coarsening fraction = 0.10
set Strategy = composition threshold, viscosity
set Refinement criteria merge operation = max
set Run postprocessors on initial refinement = false

subsection Composition threshold
set Compositional field thresholds = 0.25, 2, 0.25, 2
end
end

subsection Postprocess
set List of postprocessors = visualization, temperature statistics, composition statistics, particles

subsection Visualization
set List of output variables = material properties, named additional outputs, strain rate
set Time between graphical output = 3
end

subsection Particles
set Number of particles = 500000
set Minimum particles per cell = 5
set Time between data output = 0
set Data output format = vtu
set List of particle properties = velocity, composition water, pT path
set Interpolation scheme = cell average
set Update ghost particles = true
set Particle generator name = random uniform

subsection Ascii file 
  set Data directory  = /home/mainakchakraborty/bin/aspect/data/particle/generator/ascii/
  set Data file name  = lookup.txt
end

end

end

So it stops after mesh refinement.