Merge pull request #6697 from tjhei/maintenance_oct2025

update parameters

Author: gassmoeller

doc/parameter_view/parameters.xml

@@ -137,7 +137,7 @@ For an in-depth discussion of these issues and a quantitative evaluation of the
 
 **Pattern:** [Selection entropy viscosity|SUPG ]
 
-**Documentation:** Select the method for stabilizing the advection equation. The original method implemented is ’entropy viscosity’ as described in \cite {kronbichler:etal:2012}. SUPG is currently experimental.
+**Documentation:** Select the method for stabilizing the advection equation. The original method implemented is ’entropy viscosity’ as described in {cite}`kronbichler:etal:2012`. SUPG is currently experimental.
 
 (parameters:Discretization/Stabilization_20parameters/Use_20artificial_20viscosity_20smoothing)=
 ### __Parameter name:__ Use artificial viscosity smoothing

doc/sphinx/parameters/Discretization.md

@@ -13,6 +13,14 @@
 
 **Documentation:** By default, every cell needs to contain particles to use this interpolator plugin. If this parameter is set to true, cells are allowed to have no particles. In case both the current cell and its neighbors are empty, the interpolator will return 0 for the current cell’s properties.
 
+(parameters:Particles/Bandwidth)=
+### __Parameter name:__ Bandwidth
+**Default value:** 0.3
+
+**Pattern:** [Double 0.3...MAX_DOUBLE (inclusive)]
+
+**Documentation:** The bandwidth value is used to scale the kernel function when generating the point density function of particles. The bandwidth is measured as a fraction of the cells extent in one spatial dimension. For example, the default bandwidth of 0.3 represents a size equal to 30 percent of the cells size in one spatial dimension.
+
 (parameters:Particles/Integration_20scheme)=
 ### __Parameter name:__ Integration scheme
 **Default value:** rk2
@@ -161,6 +169,14 @@ The following properties are available:
 
 ‘uniform radial’: Generate a uniform distribution of particles over a spherical domain in 2d or 3d. Uniform here means the particles will be generated with an equal spacing in each spherical spatial dimension, i.e., the particles are created at positions that increase linearly with equal spacing in radius, colatitude and longitude around a certain center point. Note that in order to produce a regular distribution the number of generated particles might not exactly match the one specified in the input file.
 
+(parameters:Particles/Particle_20removal_20algorithm)=
+### __Parameter name:__ Particle removal algorithm
+**Default value:** random
+
+**Pattern:** [Selection random|point density function ]
+
+**Documentation:** Algorithm used to delete excess particles from cells. If point density function is chosen, the particle manager will generate a point density function from the locations of each particle and remove the particle whose position is at the maximum of the point density function.
+
 (parameters:Particles/Particle_20weight)=
 ### __Parameter name:__ Particle weight
 **Default value:** 10
@@ -169,6 +185,14 @@ The following properties are available:
 
 **Documentation:** Weight that is associated with the computational load of a single particle. The sum of particle weights will be added to the sum of cell weights to determine the partitioning of the mesh if the ‘repartition’ particle load balancing strategy is selected. The optimal weight depends on the used integrator and particle properties. In general for a more expensive integrator and more expensive properties a larger particle weight is recommended. Before adding the weights of particles, each cell already carries a weight of 1000 to account for the cost of field-based computations.
 
+(parameters:Particles/Point_20density_20kernel_20function)=
+### __Parameter name:__ Point density kernel function
+**Default value:** cutoff c1 dealii
+
+**Pattern:** [Selection cutoff c1 dealii|cutoff w1 dealii|uniform|triangular|gaussian ]
+
+**Documentation:** The kernel function is summed at each particle location to generate a point density function of the particle locations according to a process known as kernel density estimation. Because kernel density estimation sums the value of a kernel function centered on each point of interest to every other point in the dataset, the only parameter of each kernel function is the distance between the particles, and each kernel function only returns a single value depending on this distance. The return value of each function is also scaled by the selected bandwidth value.The gaussian function uses the gaussian distribution to generate an output from the input distance. The output of the triangular function decreases at a constant rate with increasing distance between the particles. The uniform function returns a constant value as long as the distance between particles is less than the selected bandwidth.The cutoff w1 and cutoff c1 dealii options call the deal.II functions called cutoffW1 and cutoffC1 respectively. These are functions whose return values decrease with distance. A more detailed explanation on these two function are available in the deal.II documentation.
+
 (parameters:Particles/Update_20ghost_20particles)=
 ### __Parameter name:__ Update ghost particles
 **Default value:** true

doc/sphinx/parameters/Particles.md

@@ -13,6 +13,14 @@
 
 **Documentation:** By default, every cell needs to contain particles to use this interpolator plugin. If this parameter is set to true, cells are allowed to have no particles. In case both the current cell and its neighbors are empty, the interpolator will return 0 for the current cell’s properties.
 
+(parameters:Particles_202/Bandwidth)=
+### __Parameter name:__ Bandwidth
+**Default value:** 0.3
+
+**Pattern:** [Double 0.3...MAX_DOUBLE (inclusive)]
+
+**Documentation:** The bandwidth value is used to scale the kernel function when generating the point density function of particles. The bandwidth is measured as a fraction of the cells extent in one spatial dimension. For example, the default bandwidth of 0.3 represents a size equal to 30 percent of the cells size in one spatial dimension.
+
 (parameters:Particles_202/Integration_20scheme)=
 ### __Parameter name:__ Integration scheme
 **Default value:** rk2
@@ -153,6 +161,14 @@ The following properties are available:
 
 ‘uniform radial’: Generate a uniform distribution of particles over a spherical domain in 2d or 3d. Uniform here means the particles will be generated with an equal spacing in each spherical spatial dimension, i.e., the particles are created at positions that increase linearly with equal spacing in radius, colatitude and longitude around a certain center point. Note that in order to produce a regular distribution the number of generated particles might not exactly match the one specified in the input file.
 
+(parameters:Particles_202/Particle_20removal_20algorithm)=
+### __Parameter name:__ Particle removal algorithm
+**Default value:** random
+
+**Pattern:** [Selection random|point density function ]
+
+**Documentation:** Algorithm used to delete excess particles from cells. If point density function is chosen, the particle manager will generate a point density function from the locations of each particle and remove the particle whose position is at the maximum of the point density function.
+
 (parameters:Particles_202/Particle_20weight)=
 ### __Parameter name:__ Particle weight
 **Default value:** 10
@@ -161,6 +177,14 @@ The following properties are available:
 
 **Documentation:** Weight that is associated with the computational load of a single particle. The sum of particle weights will be added to the sum of cell weights to determine the partitioning of the mesh if the ‘repartition’ particle load balancing strategy is selected. The optimal weight depends on the used integrator and particle properties. In general for a more expensive integrator and more expensive properties a larger particle weight is recommended. Before adding the weights of particles, each cell already carries a weight of 1000 to account for the cost of field-based computations.
 
+(parameters:Particles_202/Point_20density_20kernel_20function)=
+### __Parameter name:__ Point density kernel function
+**Default value:** cutoff c1 dealii
+
+**Pattern:** [Selection cutoff c1 dealii|cutoff w1 dealii|uniform|triangular|gaussian ]
+
+**Documentation:** The kernel function is summed at each particle location to generate a point density function of the particle locations according to a process known as kernel density estimation. Because kernel density estimation sums the value of a kernel function centered on each point of interest to every other point in the dataset, the only parameter of each kernel function is the distance between the particles, and each kernel function only returns a single value depending on this distance. The return value of each function is also scaled by the selected bandwidth value.The gaussian function uses the gaussian distribution to generate an output from the input distance. The output of the triangular function decreases at a constant rate with increasing distance between the particles. The uniform function returns a constant value as long as the distance between particles is less than the selected bandwidth.The cutoff w1 and cutoff c1 dealii options call the deal.II functions called cutoffW1 and cutoffC1 respectively. These are functions whose return values decrease with distance. A more detailed explanation on these two function are available in the deal.II documentation.
+
 (parameters:Particles_202/Update_20ghost_20particles)=
 ### __Parameter name:__ Update ghost particles
 **Default value:** true

doc/sphinx/parameters/Particles_202.md

@@ -1085,7 +1085,7 @@ In two space dimensions, $\mathbf n$ is simply a vector that is horizontal (we c
 
 In three space dimensions, given two horizontal, perpendicular, unit length, but otherwise arbitrarily chosen vectors $\mathbf u,\mathbf v$, we can express $\mathbf n = (\cos \alpha)\mathbf u + (\sin\alpha)\mathbf v$ where $\alpha$ maximizes the expression $\begin{align*}  f(\alpha) = \mathbf n^T \sigma_c \mathbf n  = (\mathbf u^T \sigma_c \mathbf u)(\cos\alpha)^2    +2(\mathbf u^T \sigma_c \mathbf v)(\cos\alpha)(\sin\alpha)    +(\mathbf v^T \sigma_c \mathbf v)(\sin\alpha)^2.\end{align*}$
 
-The maximum of $f(\alpha)$ is attained where $f^\prime(\alpha)=0$. Evaluating the derivative and using trigonometric identities, one finds that $\alpha$ has to satisfy the equation $\begin{align*}  \tan(2\alpha) = \frac{2.0\mathbf u^T \sigma_c \mathbf v}                          {\mathbf u^T \sigma_c \mathbf u                            - \mathbf v^T \sigma_c \mathbf v}.\end{align*}$Since the transform $\alpha\mapsto\alpha+\pi$ flips the direction of $\mathbf n$, we only need to seek a solution to this equation in the interval $\alpha\in[0,\pi)$. These are given by $\alpha_1=\frac 12 \arctan \frac{\mathbf u^T \sigma_c \mathbf v}{\mathbf u^T \sigma_c \mathbf u - \mathbf v^T \sigma_c \mathbf v}$ and $\alpha_2=\alpha_1+\frac{\pi}{2}$, one of which will correspond to a minimum and the other to a maximum of $f(\alpha)$. One checks the sign of $f^{\prime\prime}(\alpha)=-2(\mathbf u^T \sigma_c \mathbf u - \mathbf v^T \sigma_c \mathbf v)\cos(2\alpha) - 2 (\mathbf u^T \sigma_c \mathbf v) \sin(2\alpha)$ for each of these to determine the $\alpha$ that maximizes $f(\alpha)$, and from this immediately arrives at the correct form for the maximum horizontal stress $\mathbf n$.
+The maximum of $f(\alpha)$ is attained where $f^\prime(\alpha)=0$. Evaluating the derivative and using trigonometric identities, one finds that $\alpha$ has to satisfy the equation $\begin{align*}  \tan(2\alpha) = \frac{2.0\mathbf u^T \sigma_c \mathbf v}                          {\mathbf u^T \sigma_c \mathbf u                            - \mathbf v^T \sigma_c \mathbf v}.\end{align*}$Since the transform $\alpha\mapsto\alpha+\pi$ flips the direction of $\mathbf n$, we only need to seek a solution to this equation in the interval $\alpha\in[0,\pi)$. These are given by $\alpha_1=\frac 12 \arctan \frac{2\mathbf u^T \sigma_c \mathbf v}{\mathbf u^T \sigma_c \mathbf u - \mathbf v^T \sigma_c \mathbf v}$ and $\alpha_2=\alpha_1+\frac{\pi}{2}$, one of which will correspond to a minimum and the other to a maximum of $f(\alpha)$. One checks the sign of $f^{\prime\prime}(\alpha)=-2(\mathbf u^T \sigma_c \mathbf u - \mathbf v^T \sigma_c \mathbf v)\cos(2\alpha) - 4 (\mathbf u^T \sigma_c \mathbf v) \sin(2\alpha)$ for each of these to determine the $\alpha$ that maximizes $f(\alpha)$, and from this immediately arrives at the correct form for the maximum horizontal stress $\mathbf n$.
 
 The description above computes a 3d *direction* vector $\mathbf n$. If one were to scale this vector the same way as done in 2d, i.e., with the magnitude of the stress in this direction, one will typically get vectors whose length is principally determined by the hydrostatic pressure at a given location simply because the hydrostatic pressure is the largest component of the overall stress. On the other hand, the hydrostatic pressure does not determine any principal direction because it is an isotropic, anti-compressive force. As a consequence, there are often points in simulations (e.g., at the center of convection rolls) where the stress has no dominant horizontal direction, and the algorithm above will then in essence choose a random direction because the stress is approximately equal in all horizontal directions. If one scaled the output by the magnitude of the stress in this direction (i.e., approximately equal to the hydrostatic pressure at this point), one would get randomly oriented vectors at these locations with significant lengths.
 
@@ -1095,8 +1095,6 @@ Fig.~\ref{fig:max-horizontal-compressive-stress} shows a simple example for this
 
 \begin{figure}  \includegraphics[width=0.3\textwidth]    {viz/plugins/maximum_horizontal_compressive_stress/temperature.png}  \hfill  \includegraphics[width=0.3\textwidth]    {viz/plugins/maximum_horizontal_compressive_stress/velocity.png}  \hfill  \includegraphics[width=0.3\textwidth]    {viz/plugins/maximum_horizontal_compressive_stress/horizontal-stress.png}  \caption{\it Illustration of the ‘maximum horizontal     compressive stress’ visualization plugin. The left     figure shows a ridge-like temperature anomaly. Together     with no-slip boundary along all six boundaries, this     results in two convection rolls (center). The maximal     horizontal compressive strength at the bottom center     of the domain is perpendicular to the ridge because     the flow comes together there from the left and right,     yielding a compressive force in left-right direction.     At the top of the model, the flow separates outward,     leading to a *negative* compressive stress     in left-right direction; because there is no flow     in front-back direction, the compressive strength     in front-back direction is zero, making the along-ridge     direction the dominant one. At the center of the     convection rolls, both horizontal directions yield     the same stress; the plugin therefore chooses an     essentially arbitrary horizontal vector, but then     uses a zero magnitude given that the difference     between the maximal and minimal horizontal stress     is zero at these points.}  \label{fig:max-horizontal-compressive-stress}\end{figure}
 
-Note that this plugin does not take into account elastic stresses and therefore cannot be used when elasticity is switched on.
-
 Physical units: $\text{Pa}$.
 
 ‘melt fraction’: A visualization output object that generates output for the melt fraction at the temperature and pressure of the current point. If the material model computes a melt fraction, this is the quantity that will be visualized. Otherwise, a specific parametrization for batch melting (as described in the following) will be used. It does not take into account latent heat. If there are no compositional fields, or no fields called ’pyroxenite’,  this postprocessor will visualize the melt fraction of peridotite (calculated using the anhydrous model of Katz, 2003). If there is a compositional field called ’pyroxenite’, the postprocessor assumes that this compositional field is the content of pyroxenite, and will visualize the melt fraction for a mixture of peridotite and pyroxenite (using the melting model of Sobolev, 2011 for pyroxenite). All the parameters that were used in these calculations can be changed in the input file, the most relevant maybe being the mass fraction of Cpx in peridotite in the Katz melting model (Mass fraction cpx), which right now has a default of 15\%. The corresponding $p$-$T$-diagrams can be generated by running the tests melt\_postprocessor\_peridotite and melt\_postprocessor\_pyroxenite.