From 4336599d15d6279293212d5e7bd43912a4ef6d62 Mon Sep 17 00:00:00 2001 From: Jean-Paul Pelteret Date: Wed, 24 Apr 2019 17:57:58 +0200 Subject: [PATCH] Add functions to create SD subsitution maps from scalar types --- .../sd/symengine_scalar_operations.h | 484 ++++++++++++++++++ .../sd/symengine_scalar_operations.cc | 55 ++ 2 files changed, 539 insertions(+) diff --git a/include/deal.II/differentiation/sd/symengine_scalar_operations.h b/include/deal.II/differentiation/sd/symengine_scalar_operations.h index a66280d2a4..51417f7d4f 100644 --- a/include/deal.II/differentiation/sd/symengine_scalar_operations.h +++ b/include/deal.II/differentiation/sd/symengine_scalar_operations.h @@ -25,6 +25,9 @@ # include # include +# include +# include + # include # include # include @@ -36,6 +39,8 @@ namespace Differentiation { namespace SD { + namespace SE = ::SymEngine; + /** * @name Symbolic variable creation */ @@ -130,9 +135,488 @@ namespace Differentiation //@} + /** + * @name Symbolic map creation + */ + //@{ + + namespace internal + { + /** + * Return whether or not an @p entry is a valid symbol that + * we can expect to perform substitution on. + */ + bool + is_valid_substitution_symbol(const SE::Basic &entry); + } // namespace internal + + //@} + + /** + * @name Symbolic substitution map enlargement + */ + //@{ + + namespace internal + { + /** + * A convenience function to add an entry to the @p substitution_map. + * The new entry will have the key given by @p symbol with its paired + * @p value. Such maps are required to perform substitution of symbolic + * expressions, where all entries given as keys in the map are substituted + * by their value counterparts. + * + * There are cases where it is convenient to simply ignore the fact that + * we may be trying to add invalid symbols. For example, one may wish to + * add a tensor where only some entries are symbols and the others are + * zero'd. In this case we ensure that the user knows what they're doing + * by forcing them to override this safety mechanism with a template + * parameter. These types of functions are typically called in a manner + * indicating that the user knows exactly what they're passing into it. + * + * @tparam ignore_invalid_symbols A template parameter that enforces whether + * or not the @p symbol has to be a valid one or not. In the + * overwhelming majority of cases, the default value of + * false should be selected, with the result that an + * exception will be thrown if the input @p symbolic is, in fact, + * not a symbolic value or expression. + */ + template + void + add_to_substitution_map(types::substitution_map & substitution_map, + const SE::RCP &symbol, + const SE::RCP &value); + } // namespace internal + + /** + * A convenience function to add an entry to the @p substitution_map. + * The new entry will have the key given by @p symbol with its paired + * @p value. Such maps are required to perform substitution of symbolic + * expressions, with key entries being exchanges with their pair values. + * + * This function is guaranteed to create an ordering that is identical + * to the typical add_to_symbol_map() call that may be used when initially + * configuring a BatchOptimizer. + * This helps one conform to the requirement that the arguments sent into + * lambda and LLVM JIT compiled functions (created by optimizing symbolic + * expressions) (i) be the same, and (ii) have a constant ordering. + * + * @tparam ignore_invalid_symbols A template parameter that enforces whether + * or not the @p symbol has to be a valid one or not. In the + * overwhelming majority of cases, the default value of + * false should be selected, with the result that an + * exception will be thrown if the input @p symbolic is, in fact, + * not a symbolic value or expression. + * An exceptional case is, for example, when performing symbolic + * assembly on a finite element level. When extracting the symbolic + * equivalent of the shape function gradients using FEExtractors, + * the returned tensor will have some a priori determined + * zero-valued components. These trivial components are not valid + * symbols (as they are not symbolic expressions), and we would + * typically wish to guard against their (erroneous) inclusion. In + * this scenario, for convenience, one could set + * @p ignore_invalid_symbols to true and these zero-valued + * entries would be skipped over and ignored. + * + * @note In this function, the @p ValueType is somewhat arbitrary as + * it is only used to create dummy values. + */ + template + void + add_to_substitution_map(types::substitution_map &substitution_map, + const Expression & symbol, + const Expression & value); + + /** + * A convenience function to add an entry to the @p substitution_map. + * The new entry will have the key given by @p symbol with its paired + * @p value. + * + * The @p ExpressionType will be used to convert the @p value to a + * compatible SymEngine number type. + * It is therefore required that the @p ExpressionType + * 1. can be constructed from a @p ValueType, and that + * 2. it is convertible to a `const SE::RCP &`. + * + * For more context which this function is used, see the other + * \ref add_to_substitution_map(types::substitution_map &,const + * Expression &, const Expression &) function. + * + * @tparam ignore_invalid_symbols See the other + * \ref add_to_substitution_map(types::substitution_map &,const + * Expression &,const Expression &) function for a detailed + * discussion on the role of this template argument. + */ + template &>::value && + std::is_constructible::value>::type> + void + add_to_substitution_map(types::substitution_map &substitution_map, + const ExpressionType & symbol, + const ValueType & value); + + /** + * A convenience function for adding multiple entries to the + * @p substitution_map. The new entries will have the keys given in the + * @p symbols vector, each of which will be paired index-wise with its + * corresponding element in the @p values vector. + * + * The class represented by the @p ExpressionType template parameter will + * be used to convert the p value to a compatible SymEngine number type. + * It is therefore required that the @p ExpressionType + * 1. can be constructed from a @p ValueType, and that + * 2. it is convertible to a `const SE::RCP &`. + * + * For more context which this function is used, see the other + * \ref add_to_substitution_map(types::substitution_map &,const + * Expression &, const Expression &) function. + * + * @tparam ignore_invalid_symbols See the other + * \ref add_to_substitution_map(types::substitution_map &,const + * Expression &,const Expression &) function for a detailed + * discussion on the role of this template argument. + */ + template &>::value && + std::is_constructible::value>::type> + void + add_to_substitution_map(types::substitution_map & substitution_map, + const types::symbol_vector & symbols, + const std::vector &values); + + /** + * A convenience function for adding multiple entries to the + * @p substitution_map. The new entries will have the keys given in the + * @p symbols vector, each of which will be paired index-wise with its + * corresponding element in the @p values vector. It is expected that there + * are no duplicate entries between the two maps. + * + * For more context which this function is used, see the other + * \ref add_to_substitution_map(types::substitution_map &,const + * Expression &, const Expression &) function. + * + * @tparam ignore_invalid_symbols See the other + * \ref add_to_substitution_map(types::substitution_map &,const + * Expression &,const Expression &) function for a detailed + * discussion on the role of this template argument. + */ + template + void + add_to_substitution_map(types::substitution_map & substitution_map, + const types::substitution_map &symbol_values); + + /** + * A convenience function to add an entry to the @p substitution_map. + * The new entry will have the key given by the first element of + * @p symbol_value and the value given by its second element. + * It is expected that the key entry be a valid symbol or symbolic + * expression, and that the paired @p symbol_value elements are compatible + * with the other add_to_substitution_map() functions. + * + * The @p SymbolicType and its associated @ValueType need not be scalar + * types. So, for example, this function could be used to add tensor-valued + * data to the map in the following way: + * + * @code + * types::substitution_map symbol_value_map = ...; + * add_to_substitution_map( + * symbol_value_map, + * std::make_pair(Tensor<1,dim,Expression>(...), + * Tensor<1,dim,double>(...)) + * ); + * @endcode + * + * For more context which this function is used, see the other + * \ref add_to_substitution_map(types::substitution_map &,const + * Expression &, const Expression &) function. + * + * @tparam ignore_invalid_symbols See the other + * \ref add_to_substitution_map(types::substitution_map &,const + * Expression &,const Expression &) function for a detailed + * discussion on the role of this template argument. + */ + template + void + add_to_substitution_map( + types::substitution_map & substitution_map, + const std::pair &symbol_value); + + /** + * A convenience function for adding multiple entries to the + * @p substitution_map. The new entries will have the keys given by first + * entry of each element of @p symbol_values, and the values given its + * second entry. It is expected that the key entry be a valid symbols or + * symbolic expressions, and that the paired @p symbol_value elements are + * compatible with the other add_to_substitution_map() functions. + * + * The @p SymbolicType and its associated @ValueType need not be scalar + * types. So, for example, this function could be used to add tensor-valued + * data to the map in the following way: + * + * @code + * types::substitution_map symbol_value_map = ...; + * using vector_entry_t = std::vector, Tensor<1,dim,double> + * >>; + * add_to_substitution_map( + * symbol_value_map, + * vector_entry_t{ + * {Tensor<1,dim,Expression>(...), Tensor<1,dim,double>(...)}, + * {Tensor<1,dim,Expression>(...), Tensor<1,dim,double>(...)} + * }); + * @endcode + * + * For more context which this function is used, see the other + * \ref add_to_substitution_map(types::substitution_map &,const + * Expression &, const Expression &) function. + * + * @tparam ignore_invalid_symbols See the other + * \ref add_to_substitution_map(types::substitution_map &,const + * Expression &,const Expression &) function for a detailed + * discussion on the role of this template argument. + */ + template + void + add_to_substitution_map( + types::substitution_map & substitution_map, + const std::vector> &symbol_values); + + /** + * A convenience function for adding multiple entries to the + * @p substitution_map. The new entries will have the keys given by first + * entry of each element of @p symbol_values, and the values given its + * second entry, along with the addition of the @p other_symbol_values. + * It is expected that the key entry be a valid symbols or symbolic + * expressions, and that the paired @p symbol_value elements are compatible + * with the other add_to_substitution_map() functions. + * + * With this function it is possible to construct a symbolic substitution + * map from different types, so long as there exists a + * add_to_substitution_map() function with the signature corresponding to + * the pair types. An example may be as follows: + * + * @code + * types::substitution_map symbol_value_map = ...; + * add_to_substitution_map( + * symbol_value_map, + * std::make_pair(Expression(...), 3), + * std::make_pair(Tensor<1,dim,Expression>(...), + * Tensor<1,dim,float>(...)), + * std::make_pair(SymmetricTensor<2,dim,Expression>(...), + * SymmetricTensor<2,dim,double>(...))); + * @endcode + * + * It is possible to map symbolic types to other symbolic types + * using this function. For more details on this, see the other + * \ref make_symbol_value_map(const Expression &,const ValueType &) + * function. + */ + template + void + add_to_substitution_map( + types::substitution_map & substitution_map, + const std::pair &symbol_value, + const Args &... other_symbol_values); + + //@} + } // namespace SD } // namespace Differentiation + +/* -------------------- inline and template functions ------------------ */ + + +# ifndef DOXYGEN + + +namespace Differentiation +{ + namespace SD + { + /* ---------------- Symbolic substitution map enlargement --------------*/ + + namespace internal + { + /** + * There are cases where it is convenient to simply ignore + * the fact that we may be trying to add invalid symbols. + * For example, one may wish to add a tensor where only + * some entries are symbols and the others are zero'd. + * In this case we ensure that the user knows what they're + * doing by forcing them to override this safety mechanism + * with a template parameter. These types of functions are + * typically called in a manner indicating that the user + * knows exactly what they're passing into it. + */ + template + void + add_to_substitution_map(types::substitution_map & substitution_map, + const SE::RCP &symbol, + const SE::RCP &value) + { + if (ignore_invalid_symbols == false) + { + Assert( + internal::is_valid_substitution_symbol(*symbol), + ExcMessage( + "Substitution with a number that does not represent a symbol or symbolic derivative")); + Assert(substitution_map.find(symbol) == substitution_map.end(), + ExcMessage("This symbol is already in the map.")); + substitution_map[Expression(symbol)] = Expression(value); + } + else + { + if (internal::is_valid_substitution_symbol(*symbol)) + { + Assert(substitution_map.find(symbol) == substitution_map.end(), + ExcMessage("This symbol is already in the map.")); + substitution_map[Expression(symbol)] = Expression(value); + } + } + } + + } // namespace internal + + + template + void + add_to_substitution_map(types::substitution_map &substitution_map, + const Expression & symbol, + const Expression & value) + { + internal::add_to_substitution_map( + substitution_map, symbol.get_RCP(), value.get_RCP()); + } + + + template + void + add_to_substitution_map(types::substitution_map &substitution_map, + const ExpressionType & symbol, + const ValueType & value) + { + using SE_RCP_Basic = const SE::RCP &; + internal::add_to_substitution_map( + substitution_map, + static_cast(symbol), + static_cast(ExpressionType(value))); + } + + template + void + add_to_substitution_map(types::substitution_map & substitution_map, + const types::symbol_vector & symbols, + const std::vector &values) + { + Assert(symbols.size() == values.size(), + ExcMessage( + "Vector of symbols and values must be of equal length.")); + + typename types::symbol_vector::const_iterator it_s = symbols.begin(); + typename std::vector::const_iterator it_v = values.begin(); + using SE_RCP_Basic = const SE::RCP &; + for (; it_s != symbols.end(); ++it_s, ++it_v) + { + Assert(it_v != values.end(), ExcInternalError()); + internal::add_to_substitution_map( + substitution_map, + static_cast(*it_s), + static_cast(ExpressionType(*it_v))); + } + } + + + template + void + add_to_substitution_map( + types::substitution_map & substitution_map, + const std::pair &symbol_value) + { + add_to_substitution_map(substitution_map, + symbol_value.first, + symbol_value.second); + } + + + template + void + add_to_substitution_map( + types::substitution_map & substitution_map, + const std::vector> &symbol_values) + { + for (const auto &entry : symbol_values) + { + add_to_substitution_map(substitution_map, + entry); + } + } + + + template + void + add_to_substitution_map(types::substitution_map & substitution_map, + const types::substitution_map &symbol_values) + { + for (const auto &entry : symbol_values) + { + const SE::RCP &symbol = entry.first; + const SE::RCP &value = entry.second; + internal::add_to_substitution_map( + substitution_map, symbol, value); + } + } + + + template + void + add_to_substitution_map( + types::substitution_map & substitution_map, + const std::pair &symbol_value, + const Args &... other_symbol_values) + { + add_to_substitution_map(substitution_map, + symbol_value); + add_to_substitution_map(substitution_map, + other_symbol_values...); + } + + } // namespace SD +} // namespace Differentiation + + +# endif // DOXYGEN + DEAL_II_NAMESPACE_CLOSE #endif // DEAL_II_WITH_SYMENGINE diff --git a/source/differentiation/sd/symengine_scalar_operations.cc b/source/differentiation/sd/symengine_scalar_operations.cc index 7df4e22cc0..b9e898b5ac 100644 --- a/source/differentiation/sd/symengine_scalar_operations.cc +++ b/source/differentiation/sd/symengine_scalar_operations.cc @@ -23,6 +23,8 @@ # include # include +# include + DEAL_II_NAMESPACE_OPEN namespace Differentiation @@ -65,6 +67,59 @@ namespace Differentiation return func.differentiate(op); } + + /* ------------------------ Symbolic map creation ----------------------*/ + + + namespace internal + { + bool + is_valid_substitution_symbol(const SE::Basic &entry) + { + // Allow substitution of a symbol + // It is pretty clear as to why this is wanted... + if (SE::is_a(entry)) + return true; + + // Allow substitution of a function symbol + // If desired, we can transform general but undefined functional + // relationships to an explicit form that is concrete. This is + // required for a symbolic expression to be parsed by a Lambda or LLVM + // optimizer. + if (SE::is_a(entry)) + return true; + + // Allow substitution of the explicit expression of the derivative of + // an implicitly defined symbol (i.e. the result of the derivative of + // a FunctionSymbol). + if (SE::is_a(entry)) + return true; + + // When performing tensor differentiation, one may end up with a + // coefficient of one half due to symmetry operations, e.g. + // 0.5*Derivative(Qi_00(C_11, C_00, C_01), C_01) + // So we explicitly check for this exact case + if (SE::is_a(entry)) + { + const SE::Mul &entry_mul = SE::down_cast(entry); + // Check that the factor is a half... + if (SE::eq(*(entry_mul.get_coef()), *SE::real_double(0.5))) + { + // ...and that there is only one entry and that its a + // Derivative type + const SE::map_basic_basic &entry_mul_dict = + entry_mul.get_dict(); + if (entry_mul_dict.size() == 1 && + SE::is_a(*(entry_mul_dict.begin()->first))) + return true; + } + } + + return false; + } + + } // namespace internal + } // namespace SD } // namespace Differentiation -- 2.39.5