diff --git a/include/aspect/utilities.h b/include/aspect/utilities.h index 22f7b03d998..52cdceaadc2 100644 --- a/include/aspect/utilities.h +++ b/include/aspect/utilities.h @@ -825,30 +825,64 @@ namespace aspect std_cxx11::unique_ptr > lookup; }; + /** - * This function computes a factor which can be used to make sure that the - * Jacobian remains positive definite. + * This function computes the weighted average $\bar y$ of $y_i$ for a weighted p norm. This + * leads for a general p to: + * $\bar y = \left(\frac{\sum_{i=1}^k w_i y_i^p}{\sum_{i=1}^k w_i}\right)^{\frac{1}{p}}$. + * When p = 0 we take the geometric average: + * $\bar y = \exp\left(\frac{\sum_{i=1}^k w_i \log\left(y_i\right)}{\sum_{i=1}^k w_i}\right)$, + * and when $p \le -1000$ or $p \ge 1000$ we take the minimum and maximum norm respectively. + * This means that the smallest and largest value is respectively taken taken. * - * The goal of this function is to find a factor $\alpha$ so that - * $2\eta(\varepsilon(\bm u)) I \otimes I + \alpha\left[a \otimes b + b \otimes a\right]$ remains a - * positive definite matrix. Here, $a=\varepsilon(\bm u)$ is the @p strain_rate - * and $b=\frac{\partial\eta(\varepsilon(\bm u),p)}{\partial \varepsilon}$ is the derivative of the viscosity - * with respect to the strain rate and is given by @p dviscosities_dstrain_rate. Since the viscosity $\eta$ - * must be positive, there is always a value of $\alpha$ (possibly small) so that the result is a positive - * definite matrix. In the best case, we want to choose $\alpha=1$ because that corresponds to the full Newton step, - * and so the function never returns anything larger than one. + * This function has been set up to be very tolerant to strange values, such as negative weights. + * The only things we require in for the general p is that the sum of the weights and the sum of + * the weights times the values to the power p may not be smaller or equal to zero. Futhermore, + * when a value is zero, the exponent is smaller then one and the correspondent weight is non-zero, + * this corresponds to no resistance in a parallel system. This means that this 'path' will be followed, + * and we return zero. * - * The factor is defined by: - * $\frac{2\eta(\varepsilon(\bm u))}{\left[1-\frac{b:a}{\|a\| \|b\|} \right]^2\|a\|\|b\|}$. Alpha is - * reset to a maximum of one, and if it is smaller then one, a safety_factor scales the alpha to make - * sure that the 1-alpha won't get to close to zero. + * The implemented special cases (which are minimum (p <= -1000), harmonic average (p = -1), geometric + * average (p = 0), arithmetic average (p = 1), quadratic average (RMS) (p = 2), cubic average (p = 3) + * and maximum (p >= 1000) ) is, except for the harmonic and quadratic averages even more tolerant of + * negative values, because they only require the sum of weights to be non-zero. */ - template - double compute_spd_factor(const double eta, - const SymmetricTensor<2,dim> &strain_rate, - const SymmetricTensor<2,dim> &dviscosities_dstrain_rate, - const double safety_factor); + double weighted_p_norm_average (const std::vector &weights, + const std::vector &values, + const double p); + + /** + * This function computes the derivative ($\frac{\partial\bar y}{\partial x}$) of an average + * of the values $y_i(x)$ with regard to $x$, using $\frac{\partial y_i(x)}{\partial x}$. + * This leads for a general p to: + * $\frac{\partial\bar y}{\partial x} = + * \frac{1}{p}\left(\frac{\sum_{i=1}^k w_i y_i^p}{\sum_{i=1}^k w_i}\right)^{\frac{1}{p}-1} + * \frac{\sum_{i=1}^k w_i p y_i^{p-1} y'_i}{\sum_{i=1}^k w_i}$. + * When p = 0 we take the geometric average as a reference, which results in: + * $\frac{\partial\bar y}{\partial x} = + * \exp\left(\frac{\sum_{i=1}^k w_i \log\left(y_i\right)}{\sum_{i=1}^k w_i}\right) + * \frac{\sum_{i=1}^k\frac{w_i y'_i}{y_i}}{\sum_{i=1}^k w_i}$ + * and when $p \le -1000$ or $p \ge 1000$ we take the min and max norm respectively. + * This means that the derivative is taken which has the min/max value. + * + * This function has, like the function weighted_p_norm_average been set up to be very tolerant to + * strange values, such as negative weights. The only things we require in for the general p is that + * the sum of the weights and the sum of the weights times the values to the power p may not be smaller + * or equal to zero. Futhermore, when a value is zero, the exponent is smaller then one and the + * correspondent weight is non-zero, this corresponds to no resistance in a parallel system. This means + * that this 'path' will be followed, and we return the corresponding derivative. + * + * The implemented special cases (which are minimum (p <= -1000), harmonic average (p = -1), geometric + * average (p = 0), arithmetic average (p = 1), and maximum (p >= 1000) ) is, except for the harmonic + * average even more tolerant of negative values, because they only require the sum of weights to be non-zero. + */ + template + T derivative_of_weighted_p_norm_average (const double averaged_parameter, + const std::vector &weights, + const std::vector &values, + const std::vector &derivatives, + const double p); } } diff --git a/source/utilities.cc b/source/utilities.cc index 725ee178322..2f68c6b8a12 100644 --- a/source/utilities.cc +++ b/source/utilities.cc @@ -1887,28 +1887,286 @@ namespace aspect } - template - double compute_spd_factor(const double eta, - const SymmetricTensor<2,dim> &strain_rate, - const SymmetricTensor<2,dim> &dviscosities_dstrain_rate, - const double safety_factor) + + + double + weighted_p_norm_average ( const std::vector &weights, + const std::vector &values, + const double p) { - double alpha = 0; - const double norm_a_b = std::sqrt((strain_rate*strain_rate)*(dviscosities_dstrain_rate*dviscosities_dstrain_rate)); - const double contract_a_b = (strain_rate*dviscosities_dstrain_rate); - const double one_minus_part = 1 - (contract_a_b / norm_a_b); - const double denom = one_minus_part * one_minus_part * norm_a_b; - if (denom == 0) - alpha = 1.0; + //TODO: prevent devision by zero for all + double averaged_parameter = 0.0; + + // first look at the special cases which can be done faster + if (p <= -1000) + { + // Minimum + double min_value = 0; + unsigned int first_element_with_nonzero_weight = 0; + for (; first_element_with_nonzero_weight < weights.size(); ++first_element_with_nonzero_weight) + if (weights[first_element_with_nonzero_weight] > 0) + { + min_value = values[first_element_with_nonzero_weight]; + break; + } + Assert (first_element_with_nonzero_weight < weights.size(), + ExcMessage ("There are only zero (or smaller) weights in the weights vector.")); + + for (unsigned int i=first_element_with_nonzero_weight+1; i < weights.size(); ++i) + if (weights[i] != 0) + if (values[i] < min_value) + min_value = values[i]; + + return min_value; + } + else if (p == -1) + { + // Harmonic average + for (unsigned int i=0; i< weights.size(); ++i) + { + /** + * if the value is zero, we get a division by zero. To prevent this + * we look at what should happen in this case. When a value is zero, + * and the correspondent weight is non-zero, this corresponds to no + * resistance in a parallel system. This means that this will dominate, + * and we should return zero. If the value is zero and the weight is + * zero, we just ignore it. + */ + if (values[i] == 0 && weights[i] > 0) + return 0; + else if (values[i] != 0) + averaged_parameter += weights[i]/values[i]; + } + + Assert (averaged_parameter > 0, ExcMessage ("The sum of the weights/values may not be smaller or equal to zero.")); + const double sum_of_weights = std::accumulate(weights.begin(), weights.end(), 0); + return sum_of_weights/averaged_parameter; + } + else if (p == 0) + { + // Geometric average + for (unsigned int i=0; i < weights.size(); ++i) + averaged_parameter += weights[i]*std::log(values[i]); + const double sum_of_weights = std::accumulate(weights.begin(), weights.end(), 0); + Assert (sum_of_weights != 0, + ExcMessage ("The sum of the weights may not be equal to zero, because we need to divide through it.")); + return std::exp(averaged_parameter/sum_of_weights); + } + else if (p == 1) + { + // Arithmetic average + for (unsigned int i=0; i< weights.size(); ++i) + averaged_parameter += weights[i]*values[i]; + const double sum_of_weights = std::accumulate(weights.begin(), weights.end(), 0); + Assert (sum_of_weights != 0, + ExcMessage ("The sum of the weights may not be equal to zero, because we need to divide through it.")); + return averaged_parameter/sum_of_weights; + } + else if (p == 2) + { + // Quadratic average (RMS) + for (unsigned int i=0; i< weights.size(); ++i) + averaged_parameter += weights[i]*values[i]*values[i]; + const double sum_of_weights = std::accumulate(weights.begin(), weights.end(), 0); + Assert (sum_of_weights != 0, + ExcMessage ("The sum of the weights may not be equal to zero, because we need to divide through it.")); + Assert (averaged_parameter/sum_of_weights > 0, ExcMessage ("The sum of the weights is smaller or equal to zero.")); + return std::sqrt(averaged_parameter/sum_of_weights); + } + else if (p == 3) + { + // Cubic average + for (unsigned int i=0; i< weights.size(); ++i) + averaged_parameter += weights[i]*values[i]*values[i]*values[i]; + const double sum_of_weights = std::accumulate(weights.begin(), weights.end(), 0); + Assert (sum_of_weights != 0, + ExcMessage ("The sum of the weights may not be equal to zero, because we need to divide through it.")); + return std::cbrt(averaged_parameter/sum_of_weights); + } + else if (p >= 1000) + { + // Maximum + double max_value = 0; + unsigned int first_element_with_nonzero_weight = 0; + for (; first_element_with_nonzero_weight < weights.size(); ++first_element_with_nonzero_weight) + if (weights[first_element_with_nonzero_weight] > 0) + { + max_value = values[first_element_with_nonzero_weight]; + break; + } + Assert (first_element_with_nonzero_weight < weights.size(), + ExcMessage ("There are only zero (or smaller) weights in the weights vector.")); + + for (unsigned int i=first_element_with_nonzero_weight+1; i < weights.size(); ++i) + if (weights[i] != 0) + if (values[i] > max_value) + max_value = values[i]; + + return max_value; + } else { - alpha = (2.0*eta)/denom; - if (alpha >= 1.0) - alpha = 1.0; - else - alpha = std::max(0.0,safety_factor*alpha); + for (unsigned int i=0; i< weights.size(); ++i) + { + /** + * When a value is zero or smaller, the exponent is smaller then one and the + * correspondent weight is non-zero, this corresponds to no resistance in a + * parallel system. This means that this 'path' will be followed, and we + * return zero. + */ + if (values[i] <= 0 && p < 0) + return 0; + averaged_parameter += weights[i] * std::pow(values[i],p); + } + const double sum_of_weights = std::accumulate(weights.begin(), weights.end(), 0); + Assert (sum_of_weights > 0, ExcMessage ("The sum of the weights may not be smaller or equal to zero.")); + Assert (averaged_parameter > 0, + ExcMessage ("The sum of the weights times the values to the power p may not be smaller or equal to zero.")); + return std::pow(averaged_parameter/sum_of_weights, 1/p); + } + } + + + + template + T + derivative_of_weighted_p_norm_average (const double /*averaged_parameter*/, + const std::vector &weights, + const std::vector &values, + const std::vector &derivatives, + const double p) + { + // TODO: use averaged_parameter to speed up computation? + // TODO: add special cases p = 2 and p = 3 + double averaged_parameter_derivative_part_1 = 0.0; + T averaged_parameter_derivative_part_2 = T(); + + // first look at the special cases which can be done faster + if (p <= -1000) + { + // Minimum + double min_value = 0; + unsigned int element_with_minimum_value = 0; + unsigned int first_element_with_nonzero_weight = 0; + for (; first_element_with_nonzero_weight < weights.size(); ++first_element_with_nonzero_weight) + if (weights[first_element_with_nonzero_weight] > 0) + { + min_value = values[first_element_with_nonzero_weight]; + element_with_minimum_value = first_element_with_nonzero_weight; + break; + } + Assert (first_element_with_nonzero_weight < weights.size(), + ExcMessage ("There are only zero (or smaller) weights in the weights vector.")); + + for (unsigned int i=first_element_with_nonzero_weight+1; i < weights.size(); ++i) + if (weights[i] != 0) + if (values[i] < min_value) + { + min_value = values[i]; + element_with_minimum_value = i; + } + return derivatives[element_with_minimum_value]; + } + else if (p == -1) + { + // Harmonic average + for (unsigned int i=0; i< weights.size(); ++i) + { + /** + * if the value is zero, we get a division by zero. To prevent this + * we look at what should happen in this case. When a value is zero, + * and the correspondent weight is non-zero, this corresponds to no + * resistance in a parallel system. This means that this will dominate, + * and we should return this derivative. If the value is zero and the + * weight is zero, we just ignore it. + */ + if (values[i] == 0 && weights[i] > 0) + return derivatives[i]; + else if (values[i] != 0) + { + averaged_parameter_derivative_part_1 += weights[i] / values[i]; + averaged_parameter_derivative_part_2 += weights[i] * (1/(values[i] * values[i])) * derivatives[i]; + } + } + const double sum_of_weights = std::accumulate(weights.begin(), weights.end(), 0); + Assert (sum_of_weights > 0, ExcMessage ("The sum of the weights may not be smaller or equal to zero.")); + return std::pow(averaged_parameter_derivative_part_1/sum_of_weights,-2) * averaged_parameter_derivative_part_2/sum_of_weights; + } + else if (p == 0) + { + // Geometric average + for (unsigned int i=0; i < weights.size(); ++i) + { + averaged_parameter_derivative_part_1 += weights[i]*std::log(values[i]); + averaged_parameter_derivative_part_2 += weights[i]*(1/values[i])*derivatives[i]; + } + const double sum_of_weights = std::accumulate(weights.begin(), weights.end(), 0); + Assert (sum_of_weights != 0, + ExcMessage ("The sum of the weights may not be equal to zero, because we need to divide through it.")); + return std::exp(averaged_parameter_derivative_part_1/sum_of_weights) * averaged_parameter_derivative_part_2/sum_of_weights; + } + else if (p == 1) + { + // Arithmetic average + for (unsigned int i=0; i< weights.size(); ++i) + averaged_parameter_derivative_part_2 += weights[i]*derivatives[i]; + const double sum_of_weights = std::accumulate(weights.begin(), weights.end(), 0); + Assert (sum_of_weights != 0, + ExcMessage ("The sum of the weights may not be equal to zero, because we need to divide through it.")); + return averaged_parameter_derivative_part_2/sum_of_weights; + } + else if (p >= 1000) + { + // Maximum + double max_value = 0; + unsigned int element_with_maximum_value = 0; + unsigned int first_element_with_nonzero_weight = 0; + for (; first_element_with_nonzero_weight < weights.size(); ++first_element_with_nonzero_weight) + if (weights[first_element_with_nonzero_weight] > 0) + { + max_value = values[first_element_with_nonzero_weight]; + element_with_maximum_value = first_element_with_nonzero_weight; + break; + } + Assert (first_element_with_nonzero_weight < weights.size(), + ExcMessage ("There are only zero (or smaller) weights in the weights vector.")); + + for (unsigned int i=first_element_with_nonzero_weight+1; i < weights.size(); ++i) + if (weights[i] != 0) + if (values[i] > max_value) + { + max_value = values[i]; + element_with_maximum_value = i; + } + + return derivatives[element_with_maximum_value]; + } + else + { + // The general case: We can simplify the equation by stating that (1/p) * p = 1 + //TODO: This can probably be optimized by using: + //averaged_parameter_derivative_part_2 += weights[i]*values_p[i]*(1/values[i])*derivatives[i]; and + //averaged_parameter_derivative = averaged_parameter * (1/averaged_parameter_derivative_part_1) * averaged_parameter_derivative_part_2; + for (unsigned int i=0; i< weights.size(); ++i) + { + /** + * When a value is zero or smaller, the exponent is smaller then one and the + * correspondent weight is non-zero, this corresponds to no resistance in a + * parallel system. This means that this 'path' will be followed, and we + * return that derivative. + */ + if (values[i] <= 0 && p < 0) + return derivatives[i]; + averaged_parameter_derivative_part_1 += weights[i] * std::pow(values[i],p); + averaged_parameter_derivative_part_2 += weights[i] * std::pow(values[i],p-1) * derivatives[i]; + } + const double sum_of_weights = std::accumulate(weights.begin(), weights.end(), 0); + Assert (sum_of_weights > 0, ExcMessage ("The sum of the weights may not be smaller or equal to zero.")); + Assert (averaged_parameter_derivative_part_1/sum_of_weights > 0, + ExcMessage ("The sum of the weights times the values to the power p may not be smaller or equal to zero.")); + return std::pow(averaged_parameter_derivative_part_1/sum_of_weights,(1/p)-1) * averaged_parameter_derivative_part_2/sum_of_weights; + // TODO: find a way to check if value is finite for any type? Or just leave this kind of checking up to the user? } - return alpha; } // Explicit instantiations @@ -1945,14 +2203,25 @@ namespace aspect template std_cxx11::array,1> orthogonal_vectors (const Tensor<1,2> &v); template std_cxx11::array,2> orthogonal_vectors (const Tensor<1,3> &v); - template double compute_spd_factor(const double eta, - const SymmetricTensor<2,2> &strain_rate, - const SymmetricTensor<2,2> &dviscosities_dstrain_rate, - const double safety_factor); - template double compute_spd_factor(const double eta, - const SymmetricTensor<2,3> &strain_rate, - const SymmetricTensor<2,3> &dviscosities_dstrain_rate, - const double safety_factor); - + template double + derivative_of_weighted_p_norm_average (const double averaged_parameter, + const std::vector &weights, + const std::vector &values, + const std::vector &derivatives, + const double p); + + template dealii::SymmetricTensor<2, 2, double> + derivative_of_weighted_p_norm_average (const double averaged_parameter, + const std::vector &weights, + const std::vector &values, + const std::vector > &derivatives, + const double p); + + template dealii::SymmetricTensor<2, 3, double> + derivative_of_weighted_p_norm_average (const double averaged_parameter, + const std::vector &weights, + const std::vector &values, + const std::vector > &derivatives, + const double p); } } diff --git a/tests/derivative_of_weighted_p_norm_average.cc b/tests/derivative_of_weighted_p_norm_average.cc new file mode 100644 index 00000000000..8f5bd584c3d --- /dev/null +++ b/tests/derivative_of_weighted_p_norm_average.cc @@ -0,0 +1,22 @@ +#include + +int f() +{ + using namespace aspect; + + std::vector weights = {1,2}; + std::vector values = {3,4}; + std::vector derivatives = {5,6}; + std::vector p_norms = {-1000,-2.5,-2,-1,0,1,2,2.5,3,4,1000}; + + for (unsigned int i = 0; i < p_norms.size(); i++) + { + std::cout << "p = " << p_norms[i] << ", average = " << aspect::Utilities::derivative_of_weighted_p_norm_average(0,weights,values,derivatives,p_norms[i]) << std::endl; + } + + exit(0); + return 42; +} +// run this function by initializing a global variable by it +int i = f(); + diff --git a/tests/derivative_of_weighted_p_norm_average.prm b/tests/derivative_of_weighted_p_norm_average.prm new file mode 100644 index 00000000000..bcddea0d59b --- /dev/null +++ b/tests/derivative_of_weighted_p_norm_average.prm @@ -0,0 +1 @@ +set Additional shared libraries = tests/libderivative_of_weighted_p_norm_average.so diff --git a/tests/derivative_of_weighted_p_norm_average/screen-output b/tests/derivative_of_weighted_p_norm_average/screen-output new file mode 100644 index 00000000000..bffa67d7b72 --- /dev/null +++ b/tests/derivative_of_weighted_p_norm_average/screen-output @@ -0,0 +1,15 @@ +----------------------------------------------------------------------------- +----------------------------------------------------------------------------- + +Loading shared library <./libderivative_of_weighted_p_norm_average.so> +p = -1000, average = 5 +p = -2.5, average = 5.61923 +p = -2, average = 5.62637 +p = -1, average = 5.64 +p = 0, average = 5.65326 +p = 1, average = 5.66667 +p = 2, average = 5.68052 +p = 2.5, average = 5.68765 +p = 3, average = 5.69491 +p = 4, average = 5.70974 +p = 1000, average = 6 diff --git a/tests/weighted_p_norm_average.cc b/tests/weighted_p_norm_average.cc new file mode 100644 index 00000000000..ad9d4cbbc21 --- /dev/null +++ b/tests/weighted_p_norm_average.cc @@ -0,0 +1,28 @@ +//#include +//#include +//#include +//#include +#include + +//#include +//#include + +int f() +{ + using namespace aspect; + + std::vector weights = {1,1,2,2,3,3}; + std::vector values = {6,5,4,3,2,1}; + std::vector p_norms = {-1000,-2.5,-2,-1,0,1,2,2.5,3,4,1000}; + + for (unsigned int i = 0; i < p_norms.size(); i++) + { + std::cout << "p = " << p_norms[i] << ", average = " << aspect::Utilities::weighted_p_norm_average(weights,values,p_norms[i]) << std::endl; + } + + exit(0); + return 42; +} +// run this function by initializing a global variable by it +int i = f(); + diff --git a/tests/weighted_p_norm_average.prm b/tests/weighted_p_norm_average.prm new file mode 100644 index 00000000000..b721ab11236 --- /dev/null +++ b/tests/weighted_p_norm_average.prm @@ -0,0 +1 @@ +set Additional shared libraries = tests/libweighted_p_norm_average.so diff --git a/tests/weighted_p_norm_average/screen-output b/tests/weighted_p_norm_average/screen-output new file mode 100644 index 00000000000..725d386157c --- /dev/null +++ b/tests/weighted_p_norm_average/screen-output @@ -0,0 +1,15 @@ +----------------------------------------------------------------------------- +----------------------------------------------------------------------------- + +Loading shared library <./libweighted_p_norm_average.so> +p = -1000, average = 1 +p = -2.5, average = 1.59237 +p = -2, average = 1.6974 +p = -1, average = 1.98895 +p = 0, average = 2.38899 +p = 1, average = 2.83333 +p = 2, average = 3.24037 +p = 2.5, average = 3.41824 +p = 3, average = 3.57872 +p = 4, average = 3.85347 +p = 1000, average = 6