Add asserts and tests for multi interpolation

source/simulator/initial_conditions.cc

@@ -289,6 +289,11 @@ namespace aspect
  // A property component mask indicating for each particle world which particle properties need to be interpolated
  std::vector<ComponentMask> property_mask;
 
+ // Mark for each advection field if it has been found in any particle world. We need to keep track of this to:
+ // - make sure all fields tracked by particles are found in at least one particle world
+ // - make sure that we do not interpolate the same field twice from different particle worlds
+ std::vector<bool> advection_field_has_been_found(advection_fields.size(),false);
+
  for (unsigned int world_index = 0; world_index < particle_worlds.size(); ++world_index)
  {
  const Particle::Property::Manager<dim> &particle_property_manager = particle_worlds[world_index].get_property_manager();
@@ -306,6 +311,10 @@ namespace aspect
  // If not: assume we find it in another world.
  if (particle_property_manager.get_data_info().fieldname_exists(particle_property_and_component.first))
  {
+ Assert (advection_field_has_been_found[advection_field] == false,
+ ExcMessage("The same advection field is mapped to particle properties in different particle worlds. This is not supported."));
+ advection_field_has_been_found[advection_field] = true;
+
  const unsigned int particle_property_index = particle_property_manager.get_data_info().get_position_by_field_name(particle_property_and_component.first)
  + particle_property_and_component.second;
 
@@ -332,6 +341,11 @@ namespace aspect
  }
  }
 
+ for (unsigned int advection_field=0; advection_field<advection_fields.size(); ++advection_field)
+ Assert (advection_field_has_been_found[advection_field] == true,
+ ExcMessage("An compositional field is marked as advected by particles, but no particle property exists that is mapped to this compositional field. "
+ "Make sure that the particle property exists and is mapped to the compositional field in the parameter file."));
+
  LinearAlgebra::BlockVector particle_solution;
 
  particle_solution.reinit(system_rhs, false);

tests/particle_multiple_worlds_interpolation.prm

@@ -0,0 +1,102 @@
+# Test that we can activate multiple particle worlds and
+# interpolate properties from different particle worlds into
+# the compositional fields.
+
+set Dimension = 2
+set End time = 0
+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 = 2
+ set Names of fields = anomaly, function
+ set Compositional field methods = particles, particles
+ set Mapped particle properties = anomaly: function, function: initial 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)); 0.0
+ end
+end
+
+subsection Mesh refinement
+ set Initial adaptive refinement = 0
+ set Strategy = composition
+ set Initial global refinement = 1
+ set Time steps between mesh refinement = 0
+end
+
+subsection Postprocess
+ set List of postprocessors = velocity statistics, composition statistics, particles, visualization
+
+ subsection Visualization
+ set Output format = gnuplot
+ end
+
+ subsection Particles
+ set Time between data output = 70
+ set Data output format =
+ end
+end
+
+subsection Particles
+ set List of particle properties = initial composition
+ set Particle generator name = reference cell
+end
+
+subsection Particles 2
+ set List of particle properties = function
+ set Particle generator name = reference cell
+
+ subsection Function
+ set Variable names = x,z
+ set Function constants = pi=3.1415926
+ set Function expression = x
+ end
+end

tests/particle_multiple_worlds_interpolation/screen-output

@@ -0,0 +1,20 @@
+
+Number of active cells: 4 (on 2 levels)
+Number of degrees of freedom: 134 (50+9+25+25+25)
+
+*** Timestep 0: t=0 seconds, dt=0 seconds
+ Skipping temperature solve because RHS is zero.
+ Advecting particles... done.
+ Advecting particles... done.
+ Solving Stokes system... 9+0 iterations.
+
+ Postprocessing:
+ RMS, max velocity: 0.00308 m/s, 0.00465 m/s
+ Compositions min/max/mass: 0.2285/0.6856/0.4527 // 0/0.4999/0.1902
+ Number of advected particles 16, 16
+ Writing graphical output: output-particle_multiple_worlds_interpolation/solution/solution-00000
+
+Termination requested by criterion: end time
+
+
+

tests/particle_multiple_worlds_interpolation/statistics

@@ -0,0 +1,23 @@
+# 1: Time step number
+# 2: Time (seconds)
+# 3: Time step size (seconds)
+# 4: Number of mesh cells
+# 5: Number of Stokes degrees of freedom
+# 6: Number of temperature degrees of freedom
+# 7: Number of degrees of freedom for all compositions
+# 8: Iterations for temperature solver
+# 9: Iterations for Stokes solver
+# 10: Velocity iterations in Stokes preconditioner
+# 11: Schur complement iterations in Stokes preconditioner
+# 12: RMS velocity (m/s)
+# 13: Max. velocity (m/s)
+# 14: Minimal value for composition anomaly
+# 15: Maximal value for composition anomaly
+# 16: Global mass for composition anomaly
+# 17: Minimal value for composition function
+# 18: Maximal value for composition function
+# 19: Global mass for composition function
+# 20: Number of advected particles
+# 21: Number of advected particles (World 2)
+# 22: Visualization file name
+0 0.000000000000e+00 0.000000000000e+00 4 59 25 50 0 8 10 10 3.07552704e-03 4.65482045e-03 2.28550000e-01 6.85650000e-01 4.52704222e-01 0.00000000e+00 4.99913963e-01 1.90184427e-01 16 16 output-particle_multiple_worlds_interpolation/solution/solution-00000 

tests/particle_multiple_worlds_interpolation_fail_1.prm

@@ -0,0 +1,116 @@
+# Test that we can activate multiple particle worlds and
+# interpolate properties from different particle worlds into
+# the compositional fields. This test is supposed to fail
+# because one of the necessary particle properties does not
+# exist in any particle world.
+#
+# EXPECT FAILURE
+
+set Dimension = 2
+set End time = 1
+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
+
+############### Parameters describing the temperature field
+# Note: The temperature plays no role in this model
+
+subsection Initial temperature model
+ set Model name = function
+
+ subsection Function
+ set Function expression = 0
+ end
+end
+
+############### Parameters describing the compositional field
+# Note: The compositional field is what drives the flow
+# in this example
+
+subsection Compositional fields
+ set Number of fields = 2
+ set Names of fields = anomaly, function
+ set Compositional field methods = particles, particles
+ set Mapped particle properties = anomaly: function, function: initial 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)); 0.0
+ end
+end
+
+############### Parameters describing the discretization
+
+subsection Mesh refinement
+ set Initial adaptive refinement = 0
+ set Strategy = composition
+ set Initial global refinement = 4
+ set Time steps between mesh refinement = 0
+ set Coarsening fraction = 0.05
+ set Refinement fraction = 0.3
+end
+
+############### Parameters describing what to do with the solution
+
+subsection Postprocess
+ set List of postprocessors = velocity statistics, composition statistics, particles
+
+ subsection Particles
+ set Time between data output = 70
+ set Data output format = vtu
+ end
+end
+
+subsection Particles
+ set List of particle properties = #initial composition
+ set Particle generator name = reference cell
+end
+
+subsection Particles 2
+ set List of particle properties = function
+ set Particle generator name = reference cell
+
+ subsection Function
+ set Variable names = x,z
+ set Function constants = pi=3.1415926
+ set Function expression = x
+ end
+end

tests/particle_multiple_worlds_interpolation_fail_1/screen-output

@@ -0,0 +1,21 @@
+
+Number of active cells: 256 (on 5 levels)
+Number of degrees of freedom: 5,734 (2,178+289+1,089+1,089+1,089)
+
+*** Timestep 0: t=0 seconds, dt=0 seconds
+ Skipping temperature solve because RHS is zero.
+ Advecting particles... done.
+ Advecting particles... done.
+
+An error occurred in file <initial_conditions.cc> in function
+(line in output replaced by default.sh script)
+The violated condition was: 
+ advection_field_has_been_found[advection_field] == true
+Additional information: 
+ An compositional field is marked as advected by particles, but no
+ particle property exists that is mapped to this compositional field.
+ Make sure that the particle property exists and is mapped to the
+ compositional field in the parameter file.
+
+Stacktrace:
+(rest of the output replaced by default.sh script)