Enable multiple particle worlds.

include/aspect/postprocess/particles.h

@@ -213,36 +213,38 @@ namespace aspect
  std::vector<std::string> output_formats;
 
  /**
- * A list of pairs (time, pvtu_filename) that have so far been written
- * and that we will pass to DataOutInterface::write_pvd_record
+ * A map between particle world name and list of pairs of
+ * (time, pvtu_filename) that have so far been written.
+ * We will pass these lists to DataOutInterface::write_pvd_record
  * to create a description file that can make the association
  * between simulation time and corresponding file name (this
  * is done because there is no way to store the simulation
- * time inside the .pvtu or .vtu files).
+ * time inside the .pvtu or .vtu files). We store one list
+ * per particle world, because each particle world will
+ * have its own output directory and description file.
  */
- std::vector<std::pair<double,std::string>> times_and_pvtu_file_names;
+ std::map<std::string,std::vector<std::pair<double,std::string>>> times_and_pvtu_file_names;
 
  /**
- * A corresponding variable that we use for the .visit files created
- * by DataOutInterface::write_visit_record. The second part of a
- * pair contains all files that together form a time step.
- */
- std::vector<std::pair<double,std::vector<std::string>>> times_and_vtu_file_names;
-
- /**
- * A list of list of filenames, sorted by timestep, that correspond to
- * what has been created as output. This is used to create a descriptive
+ * A map between particle world name and list of list of
+ * filenames that corresponds to
+ * what has been created as output. This list is sorted by
+ * filename and is used to create a descriptive
  * .visit file for the entire simulation.
+ * We store one list per particle world, because each particle
+ * world will have its own output directory and description file.
  */
- std::vector<std::vector<std::string>> output_file_names_by_timestep;
+ std::map<std::string,std::vector<std::vector<std::string>>> output_file_names_by_timestep;
 
  /**
- * A set of data related to XDMF file sections describing the HDF5
- * heavy data files created. These contain things such as the
- * dimensions and names of data written at all steps during the
- * simulation.
+ * A map between particle world name and a list of data for
+ * the XDMF file sections describing the HDF5 files created.
+ * These XDMF data contain things such as the dimensions
+ * and names of data written at all steps during the simulation.
+ * We store one list per particle world, because each particle
+ * world will have its own output directory and description file.
  */
- std::vector<XDMFEntry> xdmf_entries;
+ std::map<std::string,std::vector<XDMFEntry>> xdmf_entries;
 
  /**
  * VTU file output supports grouping files from several CPUs into one
@@ -305,11 +307,16 @@ namespace aspect
  *
  * @param data_out The DataOut object that was used to write the
  * solutions.
+ * @param description_file_prefix The stem of the filename to be written.
+ * @param solution_file_directory The directory where the solution files
+ * are written.
  * @param solution_file_prefix The stem of the filename to be written.
  * @param filenames List of filenames for the current output from all
  * processors.
  */
  void write_description_files (const internal::ParticleOutput<dim> &data_out,
+ const std::string &description_file_prefix,
+ const std::string &solution_file_directory,
  const std::string &solution_file_prefix,
  const std::vector<std::string> &filenames);
  };

source/particle/world.cc

@@ -755,7 +755,7 @@ namespace aspect
  void
  World<dim>::declare_parameters (ParameterHandler &prm)
  {
- constexpr unsigned int number_of_particle_worlds = 1;
+ constexpr unsigned int number_of_particle_worlds = ASPECT_MAX_NUM_PARTICLE_WORLDS;
  for (unsigned int world_index = 0; world_index < number_of_particle_worlds; ++world_index)
  {
  if (world_index == 0)

source/postprocess/particles.cc

@@ -259,6 +259,8 @@ namespace aspect
  template <int dim>
  void
  Particles<dim>::write_description_files (const internal::ParticleOutput<dim> &data_out,
+ const std::string &description_file_prefix,
+ const std::string &solution_file_directory,
  const std::string &solution_file_prefix,
  const std::vector<std::string> &filenames)
  {
@@ -267,23 +269,23 @@ namespace aspect
  this->get_time());
  const std::string pvtu_filename = (solution_file_prefix +
  ".pvtu");
- std::ofstream pvtu_file (this->get_output_directory() + "particles/" +
+ std::ofstream pvtu_file (this->get_output_directory() + solution_file_directory +
  pvtu_filename);
  data_out.write_pvtu_record (pvtu_file, filenames);
 
  // now also generate a .pvd file that matches simulation
  // time and corresponding .pvtu record
- times_and_pvtu_file_names.emplace_back(time_in_years_or_seconds, "particles/"+pvtu_filename);
+ times_and_pvtu_file_names[description_file_prefix].emplace_back(time_in_years_or_seconds, solution_file_directory + pvtu_filename);
 
- const std::string pvd_filename = (this->get_output_directory() + "particles.pvd");
+ const std::string pvd_filename = (this->get_output_directory() + description_file_prefix + ".pvd");
  std::ofstream pvd_file (pvd_filename);
 
- DataOutBase::write_pvd_record (pvd_file, times_and_pvtu_file_names);
+ DataOutBase::write_pvd_record (pvd_file, times_and_pvtu_file_names[description_file_prefix]);
 
  // finally, do the same for VisIt via the .visit file for this
  // time step, as well as for all time steps together
  const std::string visit_filename = (this->get_output_directory()
- + "particles/"
+ + solution_file_directory
  + solution_file_prefix
  + ".visit");
  std::ofstream visit_file (visit_filename);
@@ -297,19 +299,19 @@ namespace aspect
  filenames_with_path.reserve(filenames.size());
  for (const auto &filename : filenames)
  {
- filenames_with_path.push_back("particles/" + filename);
+ filenames_with_path.push_back(solution_file_directory + filename);
  }
 
- output_file_names_by_timestep.push_back (filenames_with_path);
+ output_file_names_by_timestep[description_file_prefix].push_back (filenames_with_path);
  }
 
  std::ofstream global_visit_file (this->get_output_directory() +
- "particles.visit");
+ description_file_prefix + ".visit");
 
  std::vector<std::pair<double, std::vector<std::string>>> times_and_output_file_names;
- for (unsigned int timestep=0; timestep<times_and_pvtu_file_names.size(); ++timestep)
- times_and_output_file_names.emplace_back(times_and_pvtu_file_names[timestep].first,
- output_file_names_by_timestep[timestep]);
+ for (unsigned int timestep=0; timestep<times_and_pvtu_file_names[description_file_prefix].size(); ++timestep)
+ times_and_output_file_names.emplace_back(times_and_pvtu_file_names[description_file_prefix][timestep].first,
+ output_file_names_by_timestep[description_file_prefix][timestep]);
  DataOutBase::write_visit_record (global_visit_file, times_and_output_file_names);
  }
 
@@ -330,54 +332,55 @@ namespace aspect
  for (unsigned int particle_world = 0; particle_world < this->n_particle_worlds(); ++particle_world)
  {
  const Particle::World<dim> &world = this->get_particle_world(particle_world);
- statistics.add_value("Number of advected particles",world.n_global_particles());
 
- // If it's not time to generate an output file
- // return early with the number of particles that were advected
- if (this->get_time() < last_output_time + output_interval)
- {
- write_output = false;
- if (particle_world > 0)
- {
- number_of_advected_particles += ", ";
- }
- number_of_advected_particles += Utilities::int_to_string(world.n_global_particles());
+ const std::string statistics_column_name = (particle_world == 0 ?
+ "Number of advected particles" :
+ "Number of advected particles (World " + Utilities::int_to_string(particle_world+1) + ")");
 
- continue;
- }
+ statistics.add_value(statistics_column_name,world.n_global_particles());
 
- // If we do not write output
- // return early with the number of particles that were advected
- if (output_formats.size() == 0 || output_formats[0] == "none")
+ if (particle_world > 0)
  {
- // Up the next time we need output. This is relevant to correctly
- // write output after a restart if the format is changed.
- set_last_output_time (this->get_time());
+ number_of_advected_particles += ", ";
+ }
+ number_of_advected_particles += Utilities::int_to_string(world.n_global_particles());
+ }
 
-  write_output = false;
-  if (particle_world > 0)
-  {
-  number_of_advected_particles += ", ";
-  }
-  number_of_advected_particles += Utilities::int_to_string(world.n_global_particles());
+ // If it's not time to generate an output file
+ // return early with the number of particles that were advected
+ if (this->get_time() < last_output_time + output_interval)
+ {
+ write_output = false;
+ }
 
- continue;
- }
+ // If we do not write output
+ // return early with the number of particles that were advected
+ if (output_formats.size() == 0 || output_formats[0] == "none")
+ {
+ // Up the next time we need output. This is relevant to correctly
+ // write output after a restart if the format is changed.
+ set_last_output_time (this->get_time());
+
+ write_output = false;
+ }
 
- Assert(write_output, ExcMessage("Trying to write output while it was set to not do so."));
+ if (write_output == false)
+ return std::make_pair("Number of advected particles", number_of_advected_particles);
 
-  if (output_file_number == numbers::invalid_unsigned_int)
-  output_file_number = 0;
-  else
-  ++output_file_number;
+ if (output_file_number == numbers::invalid_unsigned_int)
+ output_file_number = 0;
+ else
+ ++output_file_number;
 
+ for (unsigned int particle_world = 0; particle_world < this->n_particle_worlds(); ++particle_world)
+ {
+ const Particle::World<dim> &world = this->get_particle_world(particle_world);
  std::string particles_output_base_name = "particles";
  if (particle_world > 0)
  {
  particles_output_base_name += "-" + Utilities::int_to_string(particle_world+1);
  }
 
-
  // Create the particle output
  const bool output_hdf5 = std::find(output_formats.begin(), output_formats.end(),"hdf5") != output_formats.end();
  internal::ParticleOutput<dim> data_out;
@@ -415,8 +418,8 @@ namespace aspect
  particle_file_name,
  time_in_years_or_seconds,
  this->get_mpi_communicator());
- xdmf_entries.push_back(new_xdmf_entry);
- data_out.write_xdmf_file(xdmf_entries, this->get_output_directory() + xdmf_filename,
+ xdmf_entries[particles_output_base_name].push_back(new_xdmf_entry);
+ data_out.write_xdmf_file(xdmf_entries[particles_output_base_name], this->get_output_directory() + xdmf_filename,
  this->get_mpi_communicator());
  }
  else if (output_format == "vtu")
@@ -433,7 +436,11 @@ namespace aspect
  filenames.push_back (particle_file_prefix
  + "." + Utilities::int_to_string(i, 4)
  + ".vtu");
- write_description_files (data_out, particle_file_prefix, filenames);
+ write_description_files (data_out,
+ particles_output_base_name,
+ particles_output_base_name + "/",
+ particle_file_prefix,
+ filenames);
  }
 
  const unsigned int n_processes = Utilities::MPI::n_mpi_processes(this->get_mpi_communicator());
@@ -548,14 +555,14 @@ namespace aspect
  screen_output = particle_output;
 
  // record the file base file name in the output file
- statistics.add_value ("Particle file name",
+ const std::string statistics_column_name = (particle_world == 0 ?
+ "Particle file name" :
+ "Particle file name (" + Utilities::int_to_string(particle_world+1) + ")");
+ statistics.add_value (statistics_column_name,
  particle_output);
  }
 
- if (write_output)
- return std::make_pair("Writing particle output:", screen_output);
- else
- return std::make_pair("Number of advected particles:", number_of_advected_particles);
+ return std::make_pair("Writing particle output:", screen_output);
  }
 
 

source/simulator/parameters.cc

@@ -359,7 +359,7 @@ namespace aspect
  "Name of the world builder file. If empty, the world builder is not initialized.");
 
  prm.declare_entry ("Number of particle worlds", "1",
- Patterns::Integer(0,1),
+ Patterns::Integer(0, ASPECT_MAX_NUM_PARTICLE_WORLDS),
  "The number of particle worlds to be created. The maximum number of particle worlds "
  "is set by the CMake variable `ASPECT_MAX_NUM_PARTICLE_WORLDS` and is by default 2.");
 

tests/particle_multiple_worlds.prm

@@ -0,0 +1,107 @@
+# Test that we can activate multiple particle worlds with
+# different particle properties and output directories.
+
+# MPI: 2
+
+set Dimension = 2
+set End time = 70
+set Use years in output instead of seconds = false
+set Number of particle worlds = 2
+
+subsection Geometry model
+ set Model name = box
+
+ subsection Box
+ set X extent = 0.9142
+ set Y extent = 1.0000
+ end
+end
+
+subsection Boundary velocity model
+ set Tangential velocity boundary indicators = left, right
+ set Zero velocity boundary indicators = bottom, top
+end
+
+subsection Material model
+ set Model name = simple
+
+ subsection Simple model
+ set Reference density = 1010
+ set Viscosity = 1e2
+ set Thermal expansion coefficient = 0
+ set Density differential for compositional field 1 = -10
+ end
+end
+
+subsection Gravity model
+ set Model name = vertical
+
+ subsection Vertical
+ set Magnitude = 10
+ end
+end
+
+subsection Initial temperature model
+ set Model name = function
+
+ subsection Function
+ set Function expression = 0
+ end
+end
+
+subsection Compositional fields
+ set Number of fields = 1
+ set Names of fields = anomaly
+end
+
+subsection Initial composition model
+ set Model name = function
+
+ subsection Function
+ set Variable names = x,z
+ set Function constants = pi=3.1415926
+ set Function expression = 0.5*(1+tanh((0.2+0.02*cos(pi*x/0.9142)-z)/0.02))
+ end
+end
+
+subsection Mesh refinement
+ set Initial adaptive refinement = 0
+ set Strategy = composition
+ set Initial global refinement = 4
+ set Time steps between mesh refinement = 0
+end
+
+subsection Postprocess
+ set List of postprocessors = velocity statistics, composition statistics, particles
+
+ subsection Particles
+ set Time between data output = 70
+ set Data output format = gnuplot
+ end
+end
+
+subsection Particles
+ set List of particle properties = initial composition
+
+ subsection Generator
+ subsection Probability density function
+ set Number of particles = 10
+ end
+ end
+end
+
+subsection Particles 2
+ set List of particle properties = function
+
+ subsection Function
+ set Variable names = x,z
+ set Function constants = pi=3.1415926
+ set Function expression = 0.5*(1+tanh((0.2+0.02*cos(pi*x/0.9142)-z)/0.02))
+ end
+
+ subsection Generator
+ subsection Probability density function
+ set Number of particles = 10
+ end
+ end
+end