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