Fix default probability density function expression

doc/modules/changes/20241021_gassmoeller

@@ -0,0 +1,9 @@
+Changed: The 'random uniform' particle generator no longer uses
+the input parameters in the subsection 'Probability density function'.
+Instead it uses its own subsection 'Random uniform'. Additionally,
+the default parameters of the 'probability density function' generator
+have been fixed to allow not specifying a probability density function.
+If no function is specified the 'probability density function' generator
+now behaves identically to the 'random uniform' generator.
+<br>
+(Rene Gassmoeller, Qianyi Lu, 2024/10/21)

source/particle/generator/probability_density_function.cc

@@ -55,6 +55,26 @@ namespace aspect
           {
             Functions::ParsedFunction<dim>::declare_parameters (prm, 1);
 
+            // This parameter overwrites one of the parameters in ParsedFunction
+            // with a new default value of 1.0. The original default value was 0
+            // everywhere, which is not allowed.
+            prm.declare_entry("Function expression",
+                              "1.0",
+                              Patterns::Anything(),
+                              "The formula that denotes the function you want to evaluate for "
+                              "particular values of the independent variables. This expression "
+                              "may contain any of the usual operations such as addition or "
+                              "multiplication, as well as all of the common functions such as "
+                              "`sin' or `cos'. In addition, it may contain expressions like "
+                              "`if(x>0, 1, -1)' where the expression evaluates to the second "
+                              "argument if the first argument is true, and to the third argument "
+                              "otherwise. For a full overview of possible expressions accepted "
+                              "see the documentation of the muparser library at http://muparser.beltoforion.de/."
+                              "\n\n"
+                              "If the function you are describing represents a vector-valued "
+                              "function with multiple components, then separate the expressions "
+                              "for individual components by a semicolon.");
+
             prm.declare_entry ("Number of particles", "1000",
                                Patterns::Double (0.),
                                "Total number of particles to create (not per processor or per element). "

source/particle/generator/random_uniform.cc

@@ -117,7 +117,10 @@ namespace aspect
       ASPECT_REGISTER_PARTICLE_GENERATOR(RandomUniform,
                                          "random uniform",
                                          "Generates a random uniform distribution of "
-                                         "particles over the entire simulation domain.")
+                                         "particles over the entire simulation domain. "
+                                         "This generator can be understood as the special case "
+                                         "of the 'probability density function' generator where "
+                                         "the probability density is constant over the domain.")
     }
   }
 }

tests/particle_generator_probability_default_value.prm

@@ -0,0 +1,93 @@
+# A description of the van Keken et al. benchmark using smooth initial
+# conditions instead of the discontinuous ones in
+# van-keken-discontinuous.prm. See the manual for more information.
+# This test adds particles to the cookbook parameter file, and specifically tests
+# the probability density function generator when run with default parameters
+# for the probability density function
+
+set Dimension                              = 2
+set End time                               = 70
+set Use years in output instead of seconds = false
+
+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
+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          = 2
+  set Time steps between mesh refinement = 0
+  set Coarsening fraction                = 0.05
+  set Refinement fraction                = 0.3
+end
+
+subsection Postprocess
+  set List of postprocessors = particles
+
+  subsection Particles
+    set Time between data output = 70
+    set Data output format = none
+  end
+end
+
+subsection Particles
+  set Particle generator name = probability density function
+
+  subsection Generator
+    subsection Probability density function
+      set Number of particles = 10
+    end
+  end
+end

tests/particle_generator_probability_default_value/screen-output

@@ -0,0 +1,26 @@
+
+Number of active cells: 16 (on 3 levels)
+Number of degrees of freedom: 349 (162+25+81+81)
+
+*** Timestep 0:  t=0 seconds, dt=0 seconds
+   Skipping temperature solve because RHS is zero.
+   Advecting particles...  done.
+   Solving C_1 system ... 0 iterations.
+   Solving Stokes system (GMG)... 11+0 iterations.
+
+   Postprocessing:
+     Number of advected particles: 10
+
+*** Timestep 1:  t=70 seconds, dt=70 seconds
+   Skipping temperature solve because RHS is zero.
+   Advecting particles...  done.
+   Solving C_1 system ... 9 iterations.
+   Solving Stokes system (GMG)... 11+0 iterations.
+
+   Postprocessing:
+     Number of advected particles: 10
+
+Termination requested by criterion: end time
+
+
+

tests/particle_generator_probability_default_value/statistics

@@ -0,0 +1,15 @@
+# 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 composition solver 1
+# 10: Iterations for Stokes solver
+# 11: Velocity iterations in Stokes preconditioner
+# 12: Schur complement iterations in Stokes preconditioner
+# 13: Number of advected particles
+0 0.000000000000e+00 0.000000000000e+00 16 187 81 81 0 0 10 12 12 10 
+1 7.000000000000e+01 7.000000000000e+01 16 187 81 81 0 9 10 12 12 10