From 75ae313fe1d4d4c84fb7f185df9ebf5328c99dd9 Mon Sep 17 00:00:00 2001 From: Luca Heltai Date: Sat, 14 Dec 2019 21:03:03 +0100 Subject: [PATCH] Symbolic function. --- doc/news/changes/major/20191214LucaHeltai | 6 + include/deal.II/base/function_parser.h | 5 + include/deal.II/base/symbolic_function.h | 483 ++++++++++++++++++ .../base/symbolic_function.templates.h | 301 +++++++++++ source/base/CMakeLists.txt | 2 + source/base/symbolic_function.cc | 26 + source/base/symbolic_function.inst.in | 25 + tests/base/symbolic_function_01.cc | 75 +++ ...lic_function_01.with_symengine=true.output | 22 + tests/base/symbolic_function_02.cc | 38 ++ ...lic_function_02.with_symengine=true.output | 6 + tests/base/symbolic_function_03.cc | 50 ++ ...lic_function_03.with_symengine=true.output | 11 + tests/base/symbolic_function_04.cc | 53 ++ ...lic_function_04.with_symengine=true.output | 11 + 15 files changed, 1114 insertions(+) create mode 100644 doc/news/changes/major/20191214LucaHeltai create mode 100644 include/deal.II/base/symbolic_function.h create mode 100644 include/deal.II/base/symbolic_function.templates.h create mode 100644 source/base/symbolic_function.cc create mode 100644 source/base/symbolic_function.inst.in create mode 100644 tests/base/symbolic_function_01.cc create mode 100644 tests/base/symbolic_function_01.with_symengine=true.output create mode 100644 tests/base/symbolic_function_02.cc create mode 100644 tests/base/symbolic_function_02.with_symengine=true.output create mode 100644 tests/base/symbolic_function_03.cc create mode 100644 tests/base/symbolic_function_03.with_symengine=true.output create mode 100644 tests/base/symbolic_function_04.cc create mode 100644 tests/base/symbolic_function_04.with_symengine=true.output diff --git a/doc/news/changes/major/20191214LucaHeltai b/doc/news/changes/major/20191214LucaHeltai new file mode 100644 index 0000000000..5f835f3b88 --- /dev/null +++ b/doc/news/changes/major/20191214LucaHeltai @@ -0,0 +1,6 @@ +New: The SymbolicFunction class allows one to leaverage the SymEngine library to generate +dealii::Function objects where the gradients, Laplacians, and Hessians are computed symbolically +providing also the possibility to extract the time derivative of the SymbolicFunction object +as another SymbolicFunction object. +
+(Luca Heltai, 2019/12/14) diff --git a/include/deal.II/base/function_parser.h b/include/deal.II/base/function_parser.h index 87a6a5b1d4..1dc24a9cfe 100644 --- a/include/deal.II/base/function_parser.h +++ b/include/deal.II/base/function_parser.h @@ -211,6 +211,11 @@ class Vector; * constants); * @endcode * + * @note The difference between this class and the SymbolicFunction class is + * that the SymbolicFunction class allows to compute first and second order + * derivatives (in a symbolic way), while this class computes first order + * derivatives only, using finite differences. For complicated expressions, + * this class is generally faster than SymbolicFunction. * * @ingroup functions * @author Luca Heltai, Timo Heister 2005, 2014 diff --git a/include/deal.II/base/symbolic_function.h b/include/deal.II/base/symbolic_function.h new file mode 100644 index 0000000000..c188d68464 --- /dev/null +++ b/include/deal.II/base/symbolic_function.h @@ -0,0 +1,483 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2019 by the deal.II authors +// +// This file is part of the deal.II library. +// +// The deal.II library is free software; you can use it, redistribute +// it, and/or modify it under the terms of the GNU Lesser General +// Public License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// The full text of the license can be found in the file LICENSE.md at +// the top level directory of deal.II. +// +// --------------------------------------------------------------------- + +#ifndef dealii_symbolic_function_h +#define dealii_symbolic_function_h + +#include + +#include +#include +#include +#include +#include +#include + +#include + +#include +#include +#include + +DEAL_II_NAMESPACE_OPEN + +// Forward declarations +#ifndef DOXYGEN +template +class Vector; +namespace Functions +{ + template + class SymbolicFunction; +} +#endif + +namespace Functions +{ +#ifdef DEAL_II_WITH_SYMENGINE + /** + * A Function class that leverages symbolic differentiation to compute + * gradients, Laplacians, Hessians, and time derivatives. + * + * This class can be used to define functions using methods provided by the + * Differentiation::SD namespace. In particular, one can define a symbolic + * evaluation point (the argument of the function), as well as a symbolic + * expression. + * + * The symbolic gradients and the symbolic Hessians are computed at + * construction time, and when a substitution in the symbolic functions is + * requested by the user using the method update_user_substitution_map(). + * + * Whenever one of the evaluation methods is called, a substitution is + * attempted with the coordinate symbols argument replaced by the evaluation + * point and the symbolic time replaced by the current time, as returned by + * the get_time() method. The user has to make sure that at evaluation time + * argument substitution provides a fully evaluated expression (i.e., no other + * symbols are contained in the function expression, except numerical values), + * or an exception will be thrown. Additional symbols can be partially + * evaluated or substituted by storing them in a user supplied substitution + * maps, that can be updated by calling update_user_substitution_map() or the + * set_additional_function_arguments() methods. + * + * The simplest use case of this class is given in the following example: + * @code + * SymbolicFunction<2> fun("x^2+y; t*x*y"); + * fun.set_time(3.0); + * Point<2> p(1.0, 2.0); + * + * auto a = fun.value(p, / * component * / 0); // a = 3.0 + * auto b = fun.value(p, / * component * / 1); // b = 6.0 + * + * auto df_dt = fun.time_derivative(); + * + * auto c = df_dt.value(p, / * component * / 0); // c = 0.0 + * auto d = df_dt.value(p, / * component * / 1); // d = 2.0 + * @endcode + * where a Function with two components is defined using a string containing + * their expressions separated by semicolons. + * + * A more involved example, that explicitly uses + * Differentiation::SD::Expression objects, is given by + * @code + * using namespace Differentiation::SD; + * // Create a position Tensor<1,2,Differentiation::SD::Expression> + * // with symbols "x" and "y", and the symbol "t" + * const auto x = SymbolicFunction<2>::get_default_coordinate_symbols(); + * const auto t = make_symbol("t"); + * + * // Use directly x[0] (the symbol "x"), x[1] (the symbol "y"), and t + * // (the symbol "t"). + * Expression f = std::sin(x[0])*std::cos(x[1])*std::sin(t); + * // Alternatively, you can achieve the same result parsing a string: + * // Expression f("sin(x)*cos(y)*sin(t)", true); + * SymbolicFunction<2> function({f}, x); + * + * // Evaluate the function, its gradient, and its Laplacian + * Point<2> p(1.0, 2.0); + * auto fp = function.value(p); + * auto gradfp = function.gradient(p); + * auto lapfp = function.laplacian(p); + * + * // Evaluate the time derivative of the function, its gradient, and its + * // Laplacian + * auto time_derivative = function.time_derivative(); + * auto dt_fp = time_derivative.value(p); + * auto dt_gradfp = time_derivative.gradient(p); + * auto dt_lapfp = time_derivative.laplacian(p); + * @endcode + * + * Partial substitution is possible (i.e., you can define the function using + * additional symbols). However, as soon as you evaluate the function, you + * have to make sure that all extraneous symbols (i.e., those not referring + * to the spacial @p coordinate_symbols or to the @p time_symbol variable) + * have been substituted with numerical values, or expressions of the spatial + * or temporal argument, by calling the update_user_substitution_map() or the + * set_additional_function_arguments() methods. + * + * If your function requires additional arguments to be evaluated, you can + * specify them by calling the set_additional_function_arguments() method. + * + * If you call update_user_substitution_map() and + * set_additional_function_arguments() with the same argument, the effect on + * the function evaluation will be the same, however, the internal behaviour + * and function derivatives will be different. The method + * update_user_substitution_map() performs the substitution once (the first + * time it is required), and then stores internally a copy of the resulting + * expression, together with its derivatives (if required). These are then + * used in all subsequent evaluations. Calling + * set_additional_function_arguments() will evaluate the passed + * subsitution map on the fly during evaluation time, *after* all + * derivatives have been computed. + * + * @note The difference between this class and the FunctionParser class is + * that this class allows to compute first and second order derivatives (in a + * symbolic way), while the FunctionParser class computes first order + * derivatives only, using finite differences. For complicated expressions, + * this class may be slower than the FunctionParser class. + * + * @ingroup functions + * @author Luca Heltai 2019 + */ + template + class SymbolicFunction : public Function + { + public: + /** + * Constructor. + * + * The resulting Function object will have as many components as there + * are entries in the vector of symbolic expressions @p function. + * + * The vector @p function should contain a list of symbolic expression + * involving the coordinate symbols argument @p coordinate_symbols and + * possibly the symbolic time argument @p time_symbol. It is possible to + * define it in terms of other symbols, as long as the optional parameter + * @p user_substitution_map replaces all symbols except + * @p coordinate_symbols and @p time_symbol. + * This is useful if, for example, you want to express formulas in terms of + * material parameters that you want to name symbolically, rather than + * through their numeric values when defining the formula, or when you want + * to express your formula in terms of polar coordinates rather than + * cartesian ones, and you want the symbolic engine to compute the + * derivatives for you. + * You may later update the symbol map contained in @p user_substitution_map + * by calling update_user_substitution_map(). + * + * @param function A vector of symbolic expressions of type + * Differentiation::SD::Expression, representing the components of this + * Function. + * + * @param coordinate_symbols A tensor of symbols representing coordinates, + * used as input argument in the symbolic expressions contained in the + * @p function vector. The default @p coordinate_symbols is a + * Tensor<1,dim,Differentiation::SD::Expression> + * containing the symbols "x" for `dim` equal to one, "x", "y" for `dim` + * equal to two, and "x", "y", "z" for `dim` equal to three. + * + * @param time_symbol A symbolic variable representing time. It defaults + * to a symbolic variable named "t". + * + * @param user_substitution_map Any other symbol that may be contained in + * the symbolic function needs to be specified in this map. The map may be + * empty, and the functions may still contain unevaluated symbols, provided + * that you call update_user_substitution_map() and provide a replacement of + * all symbols except @p coordinate_symbols and @p time_symbol before any + * evaluation occurs. + */ + SymbolicFunction( + const std::vector &function, + const Tensor<1, dim, Differentiation::SD::Expression> + &coordinate_symbols = get_default_coordinate_symbols(), + const Differentiation::SD::Expression &time_symbol = + Differentiation::SD::make_symbol("t"), + const Differentiation::SD::types::substitution_map + &user_substitution_map = {}); + + /** + * Constructor that takes a single string that describes the function + * expression as a semicolon separated list of expressions. + * + * The symbolic expression can use the default argument and the default + * symbolic time variable, plus any additional symbols that you may + * need, provided that you update the user substitution map that substitutes + * all of them before you try to evaluate the function or its derivatives, + * by calling update_user_substitution_map(), and that you provide all the + * additional function arguments of your function using the method + * set_additional_function_arguments(). + */ + SymbolicFunction(const std::string &expressions); + + /** + * Store and apply the substitution map @p substitutions to each symbolic + * component of this Function object. + * + * Notice that this method will trigger a recomputation of the + * gradients, Hessians, and Laplacians of each component. + */ + void + update_user_substitution_map( + const Differentiation::SD::types::substitution_map &substitutions); + + /** + * Set the additional @p arguments to be substituted in next evaluation + * step. + * + * Notice that the @p arguments are substituted *after* evaluating the + * @p permanent_user_substitution_map, and after all derivatives are + * computed. If the additional arguments you pass still depend on the + * coordinate or time symbols, then evaluation of derivatives will result in + * a partial derivative evaluation. + * + * This method provides a way to evaluate functions that depend on more + * arguments than simply the coordinates and time. If you want to compute + * the total derivative w.r.t. to complicated symbolic expressions, you + * should call update_user_substitution_map() instead. + */ + void + set_additional_function_arguments( + const Differentiation::SD::types::substitution_map &arguments); + + /** + * Return a tensor of coordinate symbols that can be used to define the + * expressions of this symbolic function object. + * + * The default argument is a Tensor<1,dim,Differentiation::SD::Expression> + * containing the symbols "x" for `dim` equal to one, "x", "y" for `dim` + * equal to two, and "x", "y", "z" for `dim` equal to three. + */ + static Tensor<1, dim, Differentiation::SD::Expression> + get_default_coordinate_symbols(); + + /** + * Get the actual arguments used for the coordinates in the symbolic + * function. This object does not include any user-defined arguments. + */ + const Tensor<1, dim, Differentiation::SD::Expression> & + get_coordinate_symbols() const; + + /** + * Get the actual symbolic time in use in this symbolic function. + */ + const Differentiation::SD::Expression & + get_time_symbol() const; + + /** + * Get the actual symbolic expressions used in this symbolic function. + */ + const std::vector & + get_symbolic_function_expressions() const; + + /** + * Get the currently stored @p user_substitution_map. + */ + const Differentiation::SD::types::substitution_map & + get_user_substitution_map() const; + + /** + * Return a SymbolicFunction object that represents the time derivative of + * this function. The spatial argument, the symbolic time, and the currently + * stored user substitution map are forwarded to the new function. + */ + SymbolicFunction + time_derivative() const; + + // documentation inherited from the base class + virtual RangeNumberType + value(const Point &p, const unsigned int component = 0) const override; + + // documentation inherited from the base class + virtual Tensor<1, dim, RangeNumberType> + gradient(const Point & p, + const unsigned int component = 0) const override; + + // documentation inherited from the base class + virtual RangeNumberType + laplacian(const Point & p, + const unsigned int component = 0) const override; + + // documentation inherited from the base class + virtual SymmetricTensor<2, dim, RangeNumberType> + hessian(const Point & p, + const unsigned int component = 0) const override; + + /** + * Print the stored arguments and function expression, as it would be + * evaluated when calling the method value(). + */ + template + StreamType & + print(StreamType &out) const; + + private: + /** + * Return a substitution map that replaces the argument with the values of + * @p point, the symbolic time with the value of this->get_time(), and any + * additional arguments with the substitution map given by + * @p additional_function_arguments. + */ + Differentiation::SD::types::substitution_map + create_evaluation_substitution_map(const Point &point) const; + + /** + * Recompute the symbolic value of the function, applying the user + * substitution map. This may be an expensive computation, and it is called + * only if necessary. + */ + void + update_values() const; + + /** + * Recompute the symbolic gradient of the function, applying the user + * substitution map. This may be an expensive computation, and it is called + * only if necessary. + */ + void + update_first_derivatives() const; + + /** + * Recompute the symbolic Hessian and the symbolic Lapalacian of the + * function. This may be an expensive computation, and it is called + * only if necessary. + */ + void + update_second_derivatives() const; + + /** + * The components of this symbolic function, before any subustitution took + * place. This is immutable, and generated at construction time. + * + * Before any evaluation takes place, the @p user_substitution_map is + * applied to this object, and the result is stored in the internal variable + * function. + * + * During evaluation, the @p symbolic_coordinate, the @p symbolic_time, and + * any remaining symbols are substituted with the input evaluation point, + * the current time, and the content of @p additional_function_arguments. + */ + const std::vector user_function; + + /** + * Store the user substitution map used for expression substitutions. This + * may be updated with a call to update_user_substitution_map(). Notice that + * the function may still have unresolved symbols, provided that they are + * resolved by a call to set_additional_function_arguments(). + */ + Differentiation::SD::types::substitution_map user_substitution_map; + + /** + * Store a user substitution map used for additional argument + * substitutions. This will be updated by a call to + * set_additional_function_arguments(). + */ + Differentiation::SD::types::substitution_map additional_function_arguments; + + /** + * The actual components of this symbolic function. This is obtained from + * the @p user_function, after applying the @p user_substitution_map. + */ + mutable std::vector function; + + /** + * The gradients of each component of this symbolic function. This is + * obtained by computing the symbolic gradient of the object @p function, + * that is, after applying the @p user_substitution_map to @p user_function. + */ + mutable std::vector> + function_gradient; + + /** + * The Hessians of each component of this symbolic function. This is + * obtained by computing the symbolic Hessian of the object @p function, + * that is, after applying the @p user_substitution_map to @p user_function. + */ + mutable std::vector> + function_hessian; + + /** + * The Laplacians of each component of this symbolic function. This is + * obtained by computing the symbolic Laplacian of the object @p function, + * that is, after applying the @p user_substitution_map to @p user_function. + */ + mutable std::vector function_laplacian; + + /** + * The coordinate symbols argument of the function. + */ + Tensor<1, dim, Differentiation::SD::Expression> coordinate_symbols; + + /** + * The symbolic time argument of the function. + */ + mutable Differentiation::SD::Expression time_symbol; + }; + + /** + * Allow output using the bitwise left shift operator. + */ + template + inline std::ostream & + operator<<(std::ostream &out, const SymbolicFunction &f) + { + return f.print(out); + } + + + + // Inline and template functions + template + template + StreamType & + SymbolicFunction::print(StreamType &out) const + { + for (unsigned int i = 0; i < dim; ++i) + out << coordinate_symbols[i] << ", "; + for (const auto it : additional_function_arguments) + out << it.first << ", "; + out << time_symbol << " -> " << user_function[0]; + for (unsigned int i = 1; i < user_function.size(); ++i) + out << "; " << user_function[i]; + if (!user_substitution_map.empty()) + { + out << " # ( "; + std::string sep = ""; + for (const auto it : user_substitution_map) + { + out << sep << it.first << " = " << it.second; + sep = ", "; + } + out << " )"; + } + return out; + } +#else + template + class SymbolicFunction + { + public: + SymbolicFunction() + { + AssertThrow( + false, + ExcMessage( + "This class is not available if you did not enable SymEngine " + "when compiling deal.II.")); + } + }; +#endif +} // namespace Functions + +DEAL_II_NAMESPACE_CLOSE + +#endif diff --git a/include/deal.II/base/symbolic_function.templates.h b/include/deal.II/base/symbolic_function.templates.h new file mode 100644 index 0000000000..6f30bd7441 --- /dev/null +++ b/include/deal.II/base/symbolic_function.templates.h @@ -0,0 +1,301 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2019 by the deal.II authors +// +// This file is part of the deal.II library. +// +// The deal.II library is free software; you can use it, redistribute +// it, and/or modify it under the terms of the GNU Lesser General +// Public License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// The full text of the license can be found in the file LICENSE.md at +// the top level directory of deal.II. +// +// --------------------------------------------------------------------- + +#ifndef dealii_symbolic_function_templates_h +#define dealii_symbolic_function_templates_h + +#include + +#include + +DEAL_II_NAMESPACE_OPEN + +namespace Functions +{ +#ifdef DEAL_II_WITH_SYMENGINE + + namespace + { + inline std::vector + to_expressions(const std::vector &strings) + { + std::vector ret(strings.size()); + for (unsigned int i = 0; i < strings.size(); ++i) + ret[i] = Differentiation::SD::Expression(strings[i], true); + return ret; + } + } // namespace + + + + template + SymbolicFunction::SymbolicFunction( + const std::vector & functions, + const Tensor<1, dim, Differentiation::SD::Expression> &argument, + const Differentiation::SD::Expression & time, + const Differentiation::SD::types::substitution_map &user_substitution_map) + : Function(functions.size(), 0.0) + , user_function(functions) + , user_substitution_map(user_substitution_map) + , coordinate_symbols(argument) + , time_symbol(time) + {} + + + + template + SymbolicFunction::SymbolicFunction( + const std::string &expressions) + : Function( + Utilities::split_string_list(expressions, ';').size(), + 0.0) + , user_function( + to_expressions(Utilities::split_string_list(expressions, ';'))) + , coordinate_symbols(get_default_coordinate_symbols()) + , time_symbol(Differentiation::SD::make_symbol("t")) + {} + + + + template + void + SymbolicFunction::update_values() const + { + // This is only necessary if the function vector is empty. A call to + // update_user_substitution_map() will clear all internal vectors. + if (function.empty()) + { + function.resize(user_function.size()); + for (unsigned int i = 0; i < user_function.size(); ++i) + { + function[i] = user_function[i].substitute(user_substitution_map); + } + } + } + + + + template + void + SymbolicFunction::update_first_derivatives() const + { + // This is only necessary if the gradient vector is empty. A call to + // update_user_substitution_map() will clear all internal vectors. + if (function_gradient.empty()) + { + update_values(); + function_gradient.resize(user_function.size()); + for (unsigned int i = 0; i < user_function.size(); ++i) + { + function_gradient[i] = + Differentiation::SD::differentiate(function[i], + coordinate_symbols); + } + } + } + + + + template + void + SymbolicFunction::update_second_derivatives() const + { + // This is only necessary if the gradient vector is empty. A call to + // update_user_substitution_map() will clear all internal vectors. + if (function_hessian.empty()) + { + update_first_derivatives(); + function_hessian.resize(user_function.size()); + function_laplacian.resize(user_function.size()); + + for (unsigned int i = 0; i < user_function.size(); ++i) + { + function_hessian[i] = + Differentiation::SD::differentiate(function_gradient[i], + coordinate_symbols); + function_laplacian[i] = trace(function_hessian[i]); + } + } + } + + + template + void + SymbolicFunction::update_user_substitution_map( + const Differentiation::SD::types::substitution_map &substitutions) + { + user_substitution_map = substitutions; + + // Now reset all internal vectors. + function.resize(0); + function_gradient.resize(0); + function_hessian.resize(0); + function_laplacian.resize(0); + } + + + + template + void + SymbolicFunction::set_additional_function_arguments( + const Differentiation::SD::types::substitution_map &arguments) + { + additional_function_arguments = arguments; + } + + + + template + Tensor<1, dim, Differentiation::SD::Expression> + SymbolicFunction::get_default_coordinate_symbols() + { + static_assert(dim <= 3, "Not implemented yet."); + const std::vector names = {"x", "y", "z"}; + + Tensor<1, dim, Differentiation::SD::Expression> x; + for (unsigned int i = 0; i < dim; ++i) + x[i] = Differentiation::SD::make_symbol(names[i]); + return x; + } + + + + template + const Tensor<1, dim, Differentiation::SD::Expression> & + SymbolicFunction::get_coordinate_symbols() const + { + return coordinate_symbols; + } + + + + template + const Differentiation::SD::Expression & + SymbolicFunction::get_time_symbol() const + { + return time_symbol; + } + + + + template + const std::vector & + SymbolicFunction::get_symbolic_function_expressions() + const + { + return user_function; + } + + + + template + const Differentiation::SD::types::substitution_map & + SymbolicFunction::get_user_substitution_map() const + { + return user_substitution_map; + } + + + + template + SymbolicFunction + SymbolicFunction::time_derivative() const + { + std::vector df_dt(user_function.size()); + for (unsigned int i = 0; i < user_function.size(); ++i) + df_dt[i] = user_function[i].differentiate(time_symbol); + return SymbolicFunction(df_dt, + coordinate_symbols, + time_symbol, + user_substitution_map); + } + + + + template + RangeNumberType + SymbolicFunction::value( + const Point & p, + const unsigned int component) const + { + update_values(); + AssertIndexRange(component, function.size()); + return Differentiation::SD::substitute_and_evaluate( + function[component], create_evaluation_substitution_map(p)); + } + + + + template + Tensor<1, dim, RangeNumberType> + SymbolicFunction::gradient( + const Point & p, + const unsigned int component) const + { + update_first_derivatives(); + AssertIndexRange(component, function.size()); + return Differentiation::SD::substitute_and_evaluate( + function_gradient[component], create_evaluation_substitution_map(p)); + } + + + + template + RangeNumberType + SymbolicFunction::laplacian( + const Point & p, + const unsigned int component) const + { + update_second_derivatives(); + AssertIndexRange(component, function.size()); + return Differentiation::SD::substitute_and_evaluate( + function_laplacian[component], create_evaluation_substitution_map(p)); + } + + + + template + SymmetricTensor<2, dim, RangeNumberType> + SymbolicFunction::hessian( + const Point & p, + const unsigned int component) const + { + update_second_derivatives(); + AssertIndexRange(component, function.size()); + return SymmetricTensor<2, dim, RangeNumberType>( + Differentiation::SD::substitute_and_evaluate( + function_hessian[component], create_evaluation_substitution_map(p))); + } + + + + template + Differentiation::SD::types::substitution_map + SymbolicFunction::create_evaluation_substitution_map( + const Point &p) const + { + auto map = additional_function_arguments; + Differentiation::SD::add_to_substitution_map( + map, + std::make_pair(time_symbol, this->get_time()), + std::make_pair(coordinate_symbols, p)); + return map; + } +#endif +} // namespace Functions + +DEAL_II_NAMESPACE_CLOSE + +#endif diff --git a/source/base/CMakeLists.txt b/source/base/CMakeLists.txt index 2bc37ce8b4..be79de63f7 100644 --- a/source/base/CMakeLists.txt +++ b/source/base/CMakeLists.txt @@ -78,6 +78,7 @@ SET(_unity_include_src quadrature_selector.cc scalar_polynomials_base.cc subscriptor.cc + symbolic_function.cc table_handler.cc tensor_function.cc tensor_function_parser.cc @@ -127,6 +128,7 @@ SET(_inst partitioner.cuda.inst.in polynomials_rannacher_turek.inst.in symmetric_tensor.inst.in + symbolic_function.inst.in tensor_function.inst.in tensor_function_parser.inst.in time_stepping.inst.in diff --git a/source/base/symbolic_function.cc b/source/base/symbolic_function.cc new file mode 100644 index 0000000000..251050d9ba --- /dev/null +++ b/source/base/symbolic_function.cc @@ -0,0 +1,26 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 1998 - 2015 by the deal.II authors +// +// This file is part of the deal.II library. +// +// The deal.II library is free software; you can use it, redistribute +// it, and/or modify it under the terms of the GNU Lesser General +// Public License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// The full text of the license can be found in the file LICENSE.md at +// the top level directory of deal.II. +// +// --------------------------------------------------------------------- + +#include + +DEAL_II_NAMESPACE_OPEN +#ifdef DEAL_II_WITH_SYMENGINE +namespace Functions +{ + // explicit instantiations +# include "symbolic_function.inst" +} // namespace Functions +#endif +DEAL_II_NAMESPACE_CLOSE diff --git a/source/base/symbolic_function.inst.in b/source/base/symbolic_function.inst.in new file mode 100644 index 0000000000..2afaa63c37 --- /dev/null +++ b/source/base/symbolic_function.inst.in @@ -0,0 +1,25 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2015 - 2019 by the deal.II authors +// +// This file is part of the deal.II library. +// +// The deal.II library is free software; you can use it, redistribute +// it, and/or modify it under the terms of the GNU Lesser General +// Public License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// The full text of the license can be found in the file LICENSE.md at +// the top level directory of deal.II. +// +// --------------------------------------------------------------------- + + +for (S : REAL_SCALARS; dim : SPACE_DIMENSIONS) + { + template class SymbolicFunction; + } + +for (S : COMPLEX_SCALARS; dim : SPACE_DIMENSIONS) + { + template class SymbolicFunction; + } diff --git a/tests/base/symbolic_function_01.cc b/tests/base/symbolic_function_01.cc new file mode 100644 index 0000000000..fd26ae62e1 --- /dev/null +++ b/tests/base/symbolic_function_01.cc @@ -0,0 +1,75 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2019 by the deal.II authors +// +// This file is part of the deal.II library. +// +// The deal.II library is free software; you can use it, redistribute +// it, and/or modify it under the terms of the GNU Lesser General +// Public License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// The full text of the license can be found in the file LICENSE.md at +// the top level directory of deal.II. +// +// --------------------------------------------------------------------- + +// Test the SymbolicFunction class + +#include + +#include + +#include "../tests.h" + +using namespace Differentiation::SD; + +template +Functions::SymbolicFunction +test(const std::string &expression) +{ + auto x = Functions::SymbolicFunction::get_default_coordinate_symbols(); + auto t = make_symbol("t"); + auto f = Expression(expression, true); + auto df = differentiate(f, x); + auto Hf = differentiate(df, x); + + // Create a Function with a single component and an evaluation point + Functions::SymbolicFunction fun({f}, x); + Point p; + for (unsigned int i = 0; i < dim; ++i) + p[i] = i + 1.0; + + // Output all symbolic stuff + deallog << "==========================================================" + << std::endl + << "dim = " << dim << ", Symengine" << std::endl + << "x: " << x << ", f: " << f << ", df: " << df << ", H:" << Hf + << std::endl; + + // Output the function and its evaluation + deallog << "SymbolicFunction: " << fun << std::endl + << "p: " << p << ", f(p): " << fun.value(p) + << ", grad(f)(p): " << fun.gradient(p) + << ", laplacian(f)(p): " << fun.laplacian(p) + << ", H(f)(p): " << fun.hessian(p) << std::endl; + + // Output its time derivative + auto fun_t = fun.time_derivative(); + deallog << "Time derivative of SymbolicFunction: " << fun_t << std::endl + << "p: " << p << ", f_t(p): " << fun_t.value(p) + << ", grad(f_t)(p): " << fun_t.gradient(p) + << ", laplacian(f_t)(p): " << fun_t.laplacian(p) + << ", H(f_t)(p): " << fun_t.hessian(p) << std::endl; + + return fun; +} + +int +main() +{ + initlog(); + + test<1>("x**2 + t*x"); + test<2>("x+y+cos(exp(y*x*t))"); + test<3>("tan(z*t**2) + atan(y/x)"); +} diff --git a/tests/base/symbolic_function_01.with_symengine=true.output b/tests/base/symbolic_function_01.with_symengine=true.output new file mode 100644 index 0000000000..7523d0d1ec --- /dev/null +++ b/tests/base/symbolic_function_01.with_symengine=true.output @@ -0,0 +1,22 @@ + +DEAL::========================================================== +DEAL::dim = 1, Symengine +DEAL::x: x, f: t*x + x**2, df: t + 2*x, H:2 +DEAL::SymbolicFunction: x, t -> t*x + x**2 +DEAL::p: 1.00000, f(p): 1.00000, grad(f)(p): 2.00000, laplacian(f)(p): 2.00000, H(f)(p): 2.00000 +DEAL::Time derivative of SymbolicFunction: x, t -> x +DEAL::p: 1.00000, f_t(p): 1.00000, grad(f_t)(p): 1.00000, laplacian(f_t)(p): 0.00000, H(f_t)(p): 0.00000 +DEAL::========================================================== +DEAL::dim = 2, Symengine +DEAL::x: x y, f: x + y + cos(exp(t*x*y)), df: 1 - t*y*exp(t*x*y)*sin(exp(t*x*y)) 1 - t*x*exp(t*x*y)*sin(exp(t*x*y)), H:-t**2*y**2*exp(t*x*y)*sin(exp(t*x*y)) - t**2*y**2*exp(2*t*x*y)*cos(exp(t*x*y)) -t*exp(t*x*y)*sin(exp(t*x*y)) - t**2*x*y*exp(t*x*y)*sin(exp(t*x*y)) - t**2*x*y*exp(2*t*x*y)*cos(exp(t*x*y)) -t*exp(t*x*y)*sin(exp(t*x*y)) - t**2*x*y*exp(t*x*y)*sin(exp(t*x*y)) - t**2*x*y*exp(2*t*x*y)*cos(exp(t*x*y)) -t**2*x**2*exp(t*x*y)*sin(exp(t*x*y)) - t**2*x**2*exp(2*t*x*y)*cos(exp(t*x*y)) +DEAL::SymbolicFunction: x, y, t -> x + y + cos(exp(t*x*y)) +DEAL::p: 1.00000 2.00000, f(p): 3.54030, grad(f)(p): 1.00000 1.00000, laplacian(f)(p): 0.00000, H(f)(p): 0.00000 0.00000 0.00000 0.00000 +DEAL::Time derivative of SymbolicFunction: x, y, t -> -x*y*exp(t*x*y)*sin(exp(t*x*y)) +DEAL::p: 1.00000 2.00000, f_t(p): -1.68294, grad(f_t)(p): -1.68294 -0.841471, laplacian(f_t)(p): 0.00000, H(f_t)(p): 0.00000 -0.841471 -0.841471 0.00000 +DEAL::========================================================== +DEAL::dim = 3, Symengine +DEAL::x: x y z, f: tan(t**2*z) + atan(y/x), df: -y/(x**2*(1 + y**2/x**2)) 1/(x*(1 + y**2/x**2)) t**2*(1 + tan(t**2*z)**2), H:-2*y**3/(x**5*(1 + y**2/x**2)**2) + 2*y/(x**3*(1 + y**2/x**2)) -1/(x**2*(1 + y**2/x**2)) + 2*y**2/(x**4*(1 + y**2/x**2)**2) 0 -1/(x**2*(1 + y**2/x**2)) + 2*y**2/(x**4*(1 + y**2/x**2)**2) -2*y/(x**3*(1 + y**2/x**2)**2) 0 0 0 2*t**4*tan(t**2*z)*(1 + tan(t**2*z)**2) +DEAL::SymbolicFunction: x, y, z, t -> tan(t**2*z) + atan(y/x) +DEAL::p: 1.00000 2.00000 3.00000, f(p): 1.10715, grad(f)(p): -0.400000 0.200000 0.00000, laplacian(f)(p): 2.77556e-17, H(f)(p): 0.160000 0.120000 0.00000 0.120000 -0.160000 0.00000 0.00000 0.00000 0.00000 +DEAL::Time derivative of SymbolicFunction: x, y, z, t -> 2*t*z*(1 + tan(t**2*z)**2) +DEAL::p: 1.00000 2.00000 3.00000, f_t(p): 0.00000, grad(f_t)(p): 0.00000 0.00000 0.00000, laplacian(f_t)(p): 0.00000, H(f_t)(p): 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 diff --git a/tests/base/symbolic_function_02.cc b/tests/base/symbolic_function_02.cc new file mode 100644 index 0000000000..3cfebae06a --- /dev/null +++ b/tests/base/symbolic_function_02.cc @@ -0,0 +1,38 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2019 by the deal.II authors +// +// This file is part of the deal.II library. +// +// The deal.II library is free software; you can use it, redistribute +// it, and/or modify it under the terms of the GNU Lesser General +// Public License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// The full text of the license can be found in the file LICENSE.md at +// the top level directory of deal.II. +// +// --------------------------------------------------------------------- + +// Test the SymbolicFunction class: string constructor, and multiple components. + +#include + +#include + +#include "../tests.h" + +using namespace Differentiation::SD; + +int +main() +{ + initlog(); + + Functions::SymbolicFunction<2> fun("x; y; t"); + deallog << fun << std::endl; + Point<2> p(1, 2); + deallog << "p: " << p << std::endl + << "f_0(p): " << fun.value(p, 0) << std::endl + << "f_1(p): " << fun.value(p, 1) << std::endl + << "f_2(p): " << fun.value(p, 2) << std::endl; +} diff --git a/tests/base/symbolic_function_02.with_symengine=true.output b/tests/base/symbolic_function_02.with_symengine=true.output new file mode 100644 index 0000000000..7c07137dd5 --- /dev/null +++ b/tests/base/symbolic_function_02.with_symengine=true.output @@ -0,0 +1,6 @@ + +DEAL::x, y, t -> x; y; t +DEAL::p: 1.00000 2.00000 +DEAL::f_0(p): 1.00000 +DEAL::f_1(p): 2.00000 +DEAL::f_2(p): 0.00000 diff --git a/tests/base/symbolic_function_03.cc b/tests/base/symbolic_function_03.cc new file mode 100644 index 0000000000..76c2987883 --- /dev/null +++ b/tests/base/symbolic_function_03.cc @@ -0,0 +1,50 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2019 by the deal.II authors +// +// This file is part of the deal.II library. +// +// The deal.II library is free software; you can use it, redistribute +// it, and/or modify it under the terms of the GNU Lesser General +// Public License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// The full text of the license can be found in the file LICENSE.md at +// the top level directory of deal.II. +// +// --------------------------------------------------------------------- + +// Test the SymbolicFunction class: string constructor, multiple components, +// and updating constants + +#include + +#include "../tests.h" + +using namespace Differentiation::SD; + +int +main() +{ + initlog(0); + + Functions::SymbolicFunction<2> fun("x; alpha*y; t*alpha^2"); + + auto sub = make_substitution_map("alpha", Expression("1.0", true)); + fun.update_user_substitution_map(sub); + + deallog << fun << std::endl; + Point<2> p(1, 2); + deallog << "p: " << p << std::endl + << "f_0(p): " << fun.value(p, 0) << std::endl + << "f_1(p): " << fun.value(p, 1) << std::endl + << "f_2(p): " << fun.value(p, 2) << std::endl; + + sub = make_substitution_map(std::make_pair(make_symbol("alpha"), 2.0)); + fun.update_user_substitution_map(sub); + + deallog << fun << std::endl + << "p: " << p << std::endl + << "f_0(p): " << fun.value(p, 0) << std::endl + << "f_1(p): " << fun.value(p, 1) << std::endl + << "f_2(p): " << fun.value(p, 2) << std::endl; +} diff --git a/tests/base/symbolic_function_03.with_symengine=true.output b/tests/base/symbolic_function_03.with_symengine=true.output new file mode 100644 index 0000000000..d0f24a78f2 --- /dev/null +++ b/tests/base/symbolic_function_03.with_symengine=true.output @@ -0,0 +1,11 @@ + +DEAL::x, y, t -> x; y*alpha; t*alpha**2 # ( alpha = 1.0 ) +DEAL::p: 1.00000 2.00000 +DEAL::f_0(p): 1.00000 +DEAL::f_1(p): 2.00000 +DEAL::f_2(p): 0.00000 +DEAL::x, y, t -> x; y*alpha; t*alpha**2 # ( alpha = 2.0 ) +DEAL::p: 1.00000 2.00000 +DEAL::f_0(p): 1.00000 +DEAL::f_1(p): 4.00000 +DEAL::f_2(p): 0.00000 diff --git a/tests/base/symbolic_function_04.cc b/tests/base/symbolic_function_04.cc new file mode 100644 index 0000000000..0611d39764 --- /dev/null +++ b/tests/base/symbolic_function_04.cc @@ -0,0 +1,53 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2019 by the deal.II authors +// +// This file is part of the deal.II library. +// +// The deal.II library is free software; you can use it, redistribute +// it, and/or modify it under the terms of the GNU Lesser General +// Public License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// The full text of the license can be found in the file LICENSE.md at +// the top level directory of deal.II. +// +// --------------------------------------------------------------------- + +// Test the SymbolicFunction class: string constructor, multiple components, +// updating constants, and providing additional arguments. + +#include + +#include "../tests.h" + +using namespace Differentiation::SD; + +int +main() +{ + initlog(0); + + Functions::SymbolicFunction<2> fun("x+beta; alpha*y; t*alpha^2"); + + auto sub = make_substitution_map("alpha", Expression("1.0", true)); + fun.update_user_substitution_map(sub); + + auto args = make_substitution_map(make_symbol("beta"), 0.0); + fun.set_additional_function_arguments(args); + + deallog << fun << std::endl; + Point<2> p(1, 2); + deallog << "p: " << p << std::endl + << "f_0(p): " << fun.value(p, 0) << std::endl + << "f_1(p): " << fun.value(p, 1) << std::endl + << "f_2(p): " << fun.value(p, 2) << std::endl; + + sub = make_substitution_map(std::make_pair(make_symbol("alpha"), 2.0)); + fun.update_user_substitution_map(sub); + + deallog << fun << std::endl + << "p: " << p << std::endl + << "f_0(p): " << fun.value(p, 0) << std::endl + << "f_1(p): " << fun.value(p, 1) << std::endl + << "f_2(p): " << fun.value(p, 2) << std::endl; +} diff --git a/tests/base/symbolic_function_04.with_symengine=true.output b/tests/base/symbolic_function_04.with_symengine=true.output new file mode 100644 index 0000000000..fcded2f6a4 --- /dev/null +++ b/tests/base/symbolic_function_04.with_symengine=true.output @@ -0,0 +1,11 @@ + +DEAL::x, y, beta, t -> beta + x; y*alpha; t*alpha**2 # ( alpha = 1.0 ) +DEAL::p: 1.00000 2.00000 +DEAL::f_0(p): 1.00000 +DEAL::f_1(p): 2.00000 +DEAL::f_2(p): 0.00000 +DEAL::x, y, beta, t -> beta + x; y*alpha; t*alpha**2 # ( alpha = 2.0 ) +DEAL::p: 1.00000 2.00000 +DEAL::f_0(p): 1.00000 +DEAL::f_1(p): 4.00000 +DEAL::f_2(p): 0.00000 -- 2.39.5