make sure harmonic average particle interpolator does not average neg…

source/particle/interpolator/harmonic_average.cc

@@ -81,7 +81,11 @@ namespace aspect
 
                 for (unsigned int i = 0; i < particle_properties.size(); ++i)
                   if (selected_properties[i])
-                    cell_properties[i] += 1/particle_properties[i];
+                    {
+                      AssertThrow(particle_properties[i] > 0,
+                                  ExcMessage ("All particle property values values must be greater than 0 for harmonic averaging!"));
+                      cell_properties[i] += 1/particle_properties[i];
+                    }
               }
 
             for (unsigned int i = 0; i < n_particle_properties; ++i)

tests/particle_interpolator_harmonic_average_solkz.prm

@@ -1,7 +1,13 @@
-# A test based on the SolKz benchmark that performs bilinear
-# least squares interpolation from particle properties to
+# A test based on the SolKz benchmark that performs a cellwise
+# average using harmonic averaging from particle properties to
 # compositional fields that are flagged as 'particle' advected fields.
 
+# Note that since the test performs harmonic averaging, the density
+# values being averaged all have to be positive, which required a
+# modification of the original density from the benchmark. Therefore,
+# the error in pressure is large (since the pressure no longer
+# corresponds to the analytical solution).
+
 set Additional shared libraries            = ./libparticle_interpolator_harmonic_average_solkz.so
 
 ############### Global parameters
@@ -72,12 +78,14 @@ subsection Compositional fields
   set Mapped particle properties = density_comp:function[0], viscosity_comp:function[1]
 end
 
+# The initial condition is adapted from the original solkz benchmark by adding a constant
+# value of 1 to keep all properties positive (which is necessary for harmonic averaging)
 subsection Initial composition model
   set Model name = function
   subsection Function
     set Variable names      = x,z
     set Function constants  = pi=3.1415926, eta_b=1e6
-    set Function expression = -1 * sin(2*z) * cos(3*pi*x); exp(log(eta_b)*z)
+    set Function expression = 1 - sin(2*z) * cos(3*pi*x); exp(log(eta_b)*z)
   end
 end
 
@@ -96,11 +104,13 @@ subsection Postprocess
     set Maximum particles per cell = 4096
     set Update ghost particles = true
 
+    # The initial condition is adapted from the original solkz benchmark by adding a constant
+    # value of 1 to keep all properties positive (which is necessary for harmonic averaging)
     subsection Function
       set Number of components = 2
       set Variable names      = x, z
       set Function constants  = pi=3.1415926, eta_b=1e6
-      set Function expression =  -1 * sin(2*z) * cos(3*pi*x); exp(log(eta_b)*z)
+      set Function expression = 1 - sin(2*z) * cos(3*pi*x); exp(log(eta_b)*z)
     end
 
     set Particle generator name = reference cell

tests/particle_interpolator_harmonic_average_solkz/screen-output

@@ -11,7 +11,7 @@ Number of degrees of freedom: 8,164 (2,178+289+1,089+2,304+2,304)
 
    Postprocessing:
      Number of advected particles:  4096
-     Errors u_L1, p_L1, u_L2, p_L2: 1.041640e-06, 1.716573e-03, 1.849065e-06, 3.085735e-03
+     Errors u_L1, p_L1, u_L2, p_L2: 4.231884e-07, 2.467717e-01, 6.352882e-07, 2.849793e-01
 
 Termination requested by criterion: end time