]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Add some initial tests for SymEngine
authorJean-Paul Pelteret <jppelteret@gmail.com>
Thu, 14 Mar 2019 17:27:51 +0000 (18:27 +0100)
committerJean-Paul Pelteret <jppelteret@gmail.com>
Fri, 15 Mar 2019 10:04:37 +0000 (11:04 +0100)
17 files changed:
tests/symengine/CMakeLists.txt [new file with mode: 0644]
tests/symengine/basic_01.cc [new file with mode: 0644]
tests/symengine/basic_01.output [new file with mode: 0644]
tests/symengine/basic_02.cc [new file with mode: 0644]
tests/symengine/basic_02.output [new file with mode: 0644]
tests/symengine/basic_03.cc [new file with mode: 0644]
tests/symengine/basic_03.output [new file with mode: 0644]
tests/symengine/basic_04.cc [new file with mode: 0644]
tests/symengine/basic_04.with_symengine_with_llvm=false.output [new file with mode: 0644]
tests/symengine/basic_04.with_symengine_with_llvm=true.output [new file with mode: 0644]
tests/symengine/basic_05.cc [new file with mode: 0644]
tests/symengine/basic_05.with_symengine_with_llvm=false.output [new file with mode: 0644]
tests/symengine/basic_05.with_symengine_with_llvm=true.output [new file with mode: 0644]
tests/symengine/basic_06.cc [new file with mode: 0644]
tests/symengine/basic_06.output [new file with mode: 0644]
tests/symengine/basic_07.cc [new file with mode: 0644]
tests/symengine/basic_07.output [new file with mode: 0644]

diff --git a/tests/symengine/CMakeLists.txt b/tests/symengine/CMakeLists.txt
new file mode 100644 (file)
index 0000000..46145a0
--- /dev/null
@@ -0,0 +1,6 @@
+CMAKE_MINIMUM_REQUIRED(VERSION 2.8.9)
+INCLUDE(../setup_testsubproject.cmake)
+PROJECT(testsuite CXX)
+IF(DEAL_II_WITH_SYMENGINE AND DEAL_II_WITH_CXX11)
+  DEAL_II_PICKUP_TESTS()
+ENDIF()
diff --git a/tests/symengine/basic_01.cc b/tests/symengine/basic_01.cc
new file mode 100644 (file)
index 0000000..1022f6a
--- /dev/null
@@ -0,0 +1,134 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2016 - 2018 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 at
+// the top level of the deal.II distribution.
+//
+// ---------------------------------------------------------------------
+
+
+// Check that all fundamental operations with primitive SymEngine number
+// types work as expected
+
+// References:
+// https://github.com/symengine/symengine/blob/master/symengine/tests/basic/test_basic.cpp
+// https://github.com/symengine/symengine/blob/master/symengine/tests/basic/test_number.cpp
+
+#include "../tests.h"
+
+#pragma clang diagnostic push
+#pragma clang diagnostic ignored "-Wextra-semi"
+#include <symengine/complex_double.h>
+#include <symengine/eval_double.h>
+#include <symengine/integer.h>
+#include <symengine/rational.h>
+#include <symengine/real_double.h>
+#include <symengine/visitor.h>
+#pragma clang diagnostic pop
+
+#include <fstream>
+#include <iomanip>
+
+using namespace dealii;
+namespace SE = SymEngine;
+
+SE::RCP<const SE::Number>
+make_symengine_rcp(const int &val)
+{
+  return SE::integer(val);
+}
+
+SE::RCP<const SE::Number>
+make_symengine_rcp(const float &val)
+{
+  return SE::real_double(val);
+}
+
+SE::RCP<const SE::Number>
+make_symengine_rcp(const double &val)
+{
+  return SE::real_double(val);
+}
+
+template <typename NumberType>
+SE::RCP<const SE::Number>
+make_symengine_rcp(const std::complex<NumberType> &val)
+{
+  // Build complex from two SymEngine numbers
+  return SE::Complex::from_two_nums(*make_symengine_rcp(val.real()),
+                                    *make_symengine_rcp(val.imag()));
+}
+
+template <typename NumberType>
+void
+test_number()
+{
+  SE::RCP<const SE::Number> a = make_symengine_rcp(NumberType(4.2));
+  deallog << "a: " << *a << std::endl;
+
+  SE::RCP<const SE::Basic> b(make_symengine_rcp(NumberType(2.1)));
+  deallog << "b: " << *b << std::endl;
+
+  SE::RCP<const SE::Basic> c;
+  // deallog << "c (no constructor): " << *c << std::endl; // Fails
+
+  c = SE::add(a, b);
+  deallog << "c = a+b: " << *c << std::endl;
+
+  c = SE::sub(a, b);
+  deallog << "c = a-b: " << *c << std::endl;
+
+  c = SE::mul(a, b);
+  deallog << "c = a*b: " << *c << std::endl;
+
+  c = SE::div(a, b);
+  deallog << "c = a/b: " << *c << std::endl;
+
+  c = a;
+  deallog << "c = a: " << *c << std::endl;
+
+  c = SE::add(c, a);
+  deallog << "c += a: " << *c << std::endl;
+
+  c = SE::sub(c, a);
+  deallog << "c -= a: " << *c << std::endl;
+
+  c = SE::mul(c, a);
+  deallog << "c *= a: " << *c << std::endl;
+
+  c = SE::div(c, a);
+  deallog << "c /= a: " << *c << std::endl;
+}
+
+int
+main()
+{
+  initlog();
+
+  deallog.push("Integer");
+  test_number<int>();
+  deallog.pop();
+
+  deallog.push("Float");
+  test_number<float>();
+  deallog.pop();
+
+  deallog.push("Double");
+  test_number<double>();
+  deallog.pop();
+
+  // Not available yet
+  // SymEngine::SymEngineException: Invalid Format: Expected Integer or Rational
+  // deallog << "Complex double" << std::endl;
+  // test_number<std::complex<double>>();
+  // deallog.pop();
+
+  deallog << "OK" << std::endl;
+}
diff --git a/tests/symengine/basic_01.output b/tests/symengine/basic_01.output
new file mode 100644 (file)
index 0000000..f6307bc
--- /dev/null
@@ -0,0 +1,35 @@
+
+DEAL:Integer::a: 4
+DEAL:Integer::b: 2
+DEAL:Integer::c = a+b: 6
+DEAL:Integer::c = a-b: 2
+DEAL:Integer::c = a*b: 8
+DEAL:Integer::c = a/b: 2
+DEAL:Integer::c = a: 4
+DEAL:Integer::c += a: 8
+DEAL:Integer::c -= a: 4
+DEAL:Integer::c *= a: 16
+DEAL:Integer::c /= a: 4
+DEAL:Float::a: 4.19999980926514
+DEAL:Float::b: 2.09999990463257
+DEAL:Float::c = a+b: 6.29999971389771
+DEAL:Float::c = a-b: 2.09999990463257
+DEAL:Float::c = a*b: 8.81999919891359
+DEAL:Float::c = a/b: 2.
+DEAL:Float::c = a: 4.19999980926514
+DEAL:Float::c += a: 8.39999961853027
+DEAL:Float::c -= a: 4.19999980926514
+DEAL:Float::c *= a: 17.6399983978272
+DEAL:Float::c /= a: 4.19999980926514
+DEAL:Double::a: 4.2
+DEAL:Double::b: 2.1
+DEAL:Double::c = a+b: 6.3
+DEAL:Double::c = a-b: 2.1
+DEAL:Double::c = a*b: 8.82
+DEAL:Double::c = a/b: 2.
+DEAL:Double::c = a: 4.2
+DEAL:Double::c += a: 8.4
+DEAL:Double::c -= a: 4.2
+DEAL:Double::c *= a: 17.64
+DEAL:Double::c /= a: 4.2
+DEAL::OK
diff --git a/tests/symengine/basic_02.cc b/tests/symengine/basic_02.cc
new file mode 100644 (file)
index 0000000..7284e26
--- /dev/null
@@ -0,0 +1,152 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2016 - 2018 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 at
+// the top level of the deal.II distribution.
+//
+// ---------------------------------------------------------------------
+
+
+// Check that SymEngine can perform symbolic differentiation correctly
+
+// References:
+// https://github.com/symengine/symengine/blob/master/symengine/tests/basic/test_basic.cpp
+// https://github.com/symengine/symengine/blob/master/symengine/tests/basic/test_number.cpp
+// https://github.com/symengine/symengine/blob/master/symengine/tests/basic/test_subs.cpp
+// https://github.com/symengine/symengine/blob/master/symengine/symengine_casts.h
+// https://github.com/symengine/symengine/blob/master/symengine/symengine_rcp.h
+// https://github.com/symengine/symengine/blob/master/symengine/derivative.h
+// https://github.com/symengine/symengine/blob/master/symengine/subs.h
+
+#include "../tests.h"
+
+#pragma clang diagnostic push
+#pragma clang diagnostic ignored "-Wextra-semi"
+#include <symengine/complex_double.h>
+#include <symengine/derivative.h>
+#include <symengine/eval_double.h>
+#include <symengine/integer.h>
+#include <symengine/rational.h>
+#include <symengine/real_double.h>
+#include <symengine/symbol.h>
+#include <symengine/symengine_casts.h>
+#include <symengine/visitor.h>
+#pragma clang diagnostic pop
+
+#include <fstream>
+#include <iomanip>
+
+using namespace dealii;
+namespace SE = SymEngine;
+
+SE::RCP<const SE::Number>
+make_symengine_rcp(const int &val)
+{
+  return SE::integer(val);
+}
+
+SE::RCP<const SE::Number>
+make_symengine_rcp(const float &val)
+{
+  return SE::real_double(val);
+}
+
+SE::RCP<const SE::Number>
+make_symengine_rcp(const double &val)
+{
+  return SE::real_double(val);
+}
+
+template <typename NumberType>
+SE::RCP<const SE::Number>
+make_symengine_rcp(const std::complex<NumberType> &val)
+{
+  // Build complex from two SymEngine numbers
+  return SE::Complex::from_two_nums(*make_symengine_rcp(val.real()),
+                                    *make_symengine_rcp(val.imag()));
+}
+
+SE::RCP<const SE::Symbol>
+make_symengine_rcp(const std::string &name)
+{
+  return SE::symbol(name);
+}
+
+template <typename NumberType>
+void
+test_number()
+{
+  SE::RCP<const SE::Number> a = make_symengine_rcp(NumberType(3.1));
+  SE::RCP<const SE::Basic>  b = make_symengine_rcp(NumberType(7.5));
+  deallog << "a: " << *a << std::endl;
+  deallog << "b: " << *b << std::endl;
+
+  SE::RCP<const SE::Symbol> x(make_symengine_rcp("x"));
+  SE::RCP<const SE::Basic>  y(make_symengine_rcp("y"));
+  deallog << "x: " << *x << std::endl;
+  deallog << "y: " << *y << std::endl;
+
+  SE::RCP<const SE::Basic> c;
+  // deallog << "c (no constructor): " << *c << std::endl; // Fails
+
+  // Construction of symbolic function
+  c = SE::mul(y, SE::mul(SE::sub(y, b), SE::add(a, x)));
+  deallog << "c = y*(y-b)*(a+x): " << *c << std::endl;
+
+  // Perform symbolic differentiation
+  SE::RCP<const SE::Basic> dc_dx = c->diff(x);
+  // SE::RCP<const SE::Basic> dc_dy = diff(c,SE::implicit_cast<const
+  // SE::RCP<const SE::Symbol> &>(y));
+  SE::RCP<const SE::Basic> dc_dy =
+    diff(c, SE::rcp_static_cast<const SE::Symbol>(y));
+  SE::RCP<const SE::Basic> dc_dy_2 = sdiff(c, y);
+
+  deallog << "dc_dx = a*y*(y-b): " << *dc_dx << std::endl;
+  deallog << "dc_dy = 2*y*(a+x): " << *dc_dy << std::endl;
+  deallog << "dc_dy = 2*y*(a+x): " << *dc_dy_2 << std::endl;
+
+  // Substitute values
+  SE::map_basic_basic sub_vals;
+  sub_vals[x] = make_symengine_rcp(NumberType(1));
+  sub_vals[y] = make_symengine_rcp(NumberType(2.2));
+  deallog << "dc_dx(x=1,y=2.2): " << *(dc_dx->subs(sub_vals)) << std::endl;
+
+  sub_vals[x] = make_symengine_rcp(NumberType(0.75));
+  sub_vals[y] = make_symengine_rcp(NumberType(-2.08));
+  deallog << "dc_dy(x=0.75,y=-1.08): " << *(dc_dy->subs(sub_vals)) << std::endl;
+
+  deallog << std::endl;
+}
+
+int
+main()
+{
+  initlog();
+
+  deallog.push("Integer");
+  test_number<int>();
+  deallog.pop();
+
+  deallog.push("Float");
+  test_number<float>();
+  deallog.pop();
+
+  deallog.push("Double");
+  test_number<double>();
+  deallog.pop();
+
+  // Not available yet
+  // SymEngine::SymEngineException: Invalid Format: Expected Integer or Rational
+  // deallog << "Complex double" << std::endl;
+  // test_number<std::complex<double>>();
+  // deallog.pop();
+
+  deallog << "OK" << std::endl;
+}
diff --git a/tests/symengine/basic_02.output b/tests/symengine/basic_02.output
new file mode 100644 (file)
index 0000000..a023eec
--- /dev/null
@@ -0,0 +1,35 @@
+
+DEAL:Integer::a: 3
+DEAL:Integer::b: 7
+DEAL:Integer::x: x
+DEAL:Integer::y: y
+DEAL:Integer::c = y*(y-b)*(a+x): y*(-7 + y)*(3 + x)
+DEAL:Integer::dc_dx = a*y*(y-b): y*(-7 + y)
+DEAL:Integer::dc_dy = 2*y*(a+x): y*(3 + x) + (-7 + y)*(3 + x)
+DEAL:Integer::dc_dy = 2*y*(a+x): y*(3 + x) + (-7 + y)*(3 + x)
+DEAL:Integer::dc_dx(x=1,y=2.2): -10
+DEAL:Integer::dc_dy(x=0.75,y=-1.08): -33
+DEAL:Integer::
+DEAL:Float::a: 3.09999990463257
+DEAL:Float::b: 7.5
+DEAL:Float::x: x
+DEAL:Float::y: y
+DEAL:Float::c = y*(y-b)*(a+x): y*(3.09999990463257 + x)*(-7.5 + y)
+DEAL:Float::dc_dx = a*y*(y-b): y*(-7.5 + y)
+DEAL:Float::dc_dy = 2*y*(a+x): y*(3.09999990463257 + x) + (3.09999990463257 + x)*(-7.5 + y)
+DEAL:Float::dc_dy = 2*y*(a+x): y*(3.09999990463257 + x) + (3.09999990463257 + x)*(-7.5 + y)
+DEAL:Float::dc_dx(x=1,y=2.2): -11.6600001478195
+DEAL:Float::dc_dy(x=0.75,y=-1.08): -44.8909983005524
+DEAL:Float::
+DEAL:Double::a: 3.1
+DEAL:Double::b: 7.5
+DEAL:Double::x: x
+DEAL:Double::y: y
+DEAL:Double::c = y*(y-b)*(a+x): y*(3.1 + x)*(-7.5 + y)
+DEAL:Double::dc_dx = a*y*(y-b): y*(-7.5 + y)
+DEAL:Double::dc_dy = 2*y*(a+x): y*(3.1 + x) + (3.1 + x)*(-7.5 + y)
+DEAL:Double::dc_dy = 2*y*(a+x): y*(3.1 + x) + (3.1 + x)*(-7.5 + y)
+DEAL:Double::dc_dx(x=1,y=2.2): -11.66
+DEAL:Double::dc_dy(x=0.75,y=-1.08): -44.891
+DEAL:Double::
+DEAL::OK
diff --git a/tests/symengine/basic_03.cc b/tests/symengine/basic_03.cc
new file mode 100644 (file)
index 0000000..3a11eb2
--- /dev/null
@@ -0,0 +1,211 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2016 - 2018 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 at
+// the top level of the deal.II distribution.
+//
+// ---------------------------------------------------------------------
+
+
+// Check that SymEngine can perform symbolic differentiation using function
+// symbols correctly
+
+// References:
+// https://github.com/symengine/symengine/blob/master/symengine/tests/basic/test_functions.cpp
+// https://github.com/symengine/symengine/blob/master/symengine/function.h
+// https://github.com/symengine/symengine/blob/master/symengine/tests/basic/test_basic.cpp
+// https://github.com/symengine/symengine/blob/master/symengine/tests/basic/test_number.cpp
+// https://github.com/symengine/symengine/blob/master/symengine/tests/basic/test_subs.cpp
+// https://github.com/symengine/symengine/blob/master/symengine/symengine_casts.h
+// https://github.com/symengine/symengine/blob/master/symengine/symengine_rcp.h
+// https://github.com/symengine/symengine/blob/master/symengine/derivative.h
+// https://github.com/symengine/symengine/blob/master/symengine/subs.h
+
+#include "../tests.h"
+
+#pragma clang diagnostic push
+#pragma clang diagnostic ignored "-Wextra-semi"
+#include <symengine/complex_double.h>
+#include <symengine/derivative.h>
+#include <symengine/eval_double.h>
+#include <symengine/integer.h>
+#include <symengine/rational.h>
+#include <symengine/real_double.h>
+#include <symengine/symbol.h>
+#include <symengine/symengine_casts.h>
+#include <symengine/visitor.h>
+#pragma clang diagnostic pop
+
+#include <fstream>
+#include <iomanip>
+
+using namespace dealii;
+namespace SE = SymEngine;
+
+SE::RCP<const SE::Number>
+make_symengine_rcp(const int &val)
+{
+  return SE::integer(val);
+}
+
+SE::RCP<const SE::Number>
+make_symengine_rcp(const float &val)
+{
+  return SE::real_double(val);
+}
+
+SE::RCP<const SE::Number>
+make_symengine_rcp(const double &val)
+{
+  return SE::real_double(val);
+}
+
+template <typename NumberType>
+SE::RCP<const SE::Number>
+make_symengine_rcp(const std::complex<NumberType> &val)
+{
+  // Build complex from two SymEngine numbers
+  return SE::Complex::from_two_nums(*make_symengine_rcp(val.real()),
+                                    *make_symengine_rcp(val.imag()));
+}
+
+SE::RCP<const SE::Symbol>
+make_symengine_rcp(const std::string &name)
+{
+  return SE::symbol(name);
+}
+
+SE::RCP<const SE::Basic>
+make_symengine_rcp(const std::string &                          name,
+                   const std::vector<SE::RCP<const SE::Basic>> &args)
+{
+  return SE::function_symbol(name, args);
+}
+
+SE::RCP<const SE::Basic>
+make_symengine_rcp(const std::string &name, const SE::RCP<const SE::Basic> &arg)
+{
+  return SE::function_symbol(name, arg);
+}
+
+template <typename NumberType>
+void
+test_number()
+{
+  // Modified snippet from SymEngine test
+  // {
+  //   SE::RCP<const SE::Symbol> x = SE::symbol("x");
+  //   SE::RCP<const SE::Symbol> _xi_1 = SE::symbol("_xi_1");
+  //   SE::RCP<const SE::Basic> f = SE::function_symbol("f", x);
+  //   SE::RCP<const SE::Basic> r1, r2;
+  //
+  //   f = function_symbol("f", SE::pow(x, SE::integer(2)));
+  //   deallog << "f: " << *f << std::endl;
+  //   r1 = f->diff(x);
+  //   deallog << "r1: " << *r1 << std::endl;
+  //   r2 = SE::Derivative::create(SE::function_symbol("f", _xi_1), {_xi_1});
+  //   deallog << "r2: " << *r2 << std::endl;
+  //   r2 = SE::Subs::create(r2, {{_xi_1, SE::pow(x, SE::integer(2))}});
+  //   deallog << "r2: " << *r2 << std::endl;
+  // }
+
+  // Normal symbols
+  SE::RCP<const SE::Symbol> x(make_symengine_rcp("x"));
+  SE::RCP<const SE::Basic>  y(make_symengine_rcp("y"));
+  deallog << "x: " << *x << std::endl;
+  deallog << "y: " << *y << std::endl;
+
+  // Function symbols
+  SE::RCP<const SE::Basic> f(make_symengine_rcp("f", x));
+  SE::RCP<const SE::Basic> g(make_symengine_rcp("g", {x, y}));
+  SE::RCP<const SE::Basic> h(
+    make_symengine_rcp("h", {SE::add(x, y), SE::mul(x, y)}));
+  deallog << "f: " << *f << std::endl;
+  deallog << "g: " << *g << std::endl;
+  deallog << "h: " << *h << std::endl;
+
+  // Perform symbolic differentiation
+  SE::RCP<const SE::Basic> df_dx = f->diff(x);
+  SE::RCP<const SE::Basic> dg_dx = g->diff(x);
+  SE::RCP<const SE::Basic> dh_dx = h->diff(x);
+  deallog << "df_dx: " << *df_dx << std::endl;
+  deallog << "dg_dx: " << *dg_dx << std::endl;
+  deallog << "dh_dx: " << *dh_dx << std::endl;
+  SE::RCP<const SE::Basic> df_dy =
+    f->diff(SE::rcp_static_cast<const SE::Symbol>(y));
+  SE::RCP<const SE::Basic> dg_dy =
+    g->diff(SE::rcp_static_cast<const SE::Symbol>(y));
+  SE::RCP<const SE::Basic> dh_dy =
+    h->diff(SE::rcp_static_cast<const SE::Symbol>(y));
+  deallog << "df_dy: " << *df_dy << std::endl;
+  deallog << "dg_dy: " << *dg_dy << std::endl;
+  deallog << "dh_dy: " << *dh_dy << std::endl;
+
+  // Substitute values
+  SE::map_basic_basic sub_vals;
+  sub_vals[x] = make_symengine_rcp(NumberType(1));
+  sub_vals[y] = make_symengine_rcp(NumberType(2.2));
+  sub_vals[f] = SE::pow(x, SE::integer(2));
+  sub_vals[g] = SE::add(x, y);
+  sub_vals[h] = SE::sub(x, y);
+
+  deallog << "f(x=1,y=2.2): " << *(f->subs(sub_vals))
+          << std::endl; // ->subs(sub_vals)
+  deallog << "g(x=1,y=2.2): " << *(g->subs(sub_vals))
+          << std::endl; // ->subs(sub_vals)
+  deallog << "h(x=1,y=2.2): " << *(h->subs(sub_vals))
+          << std::endl; // ->subs(sub_vals)
+  deallog << "Eval: f(x=1,y=2.2): "
+          << eval_double(*(f->subs(sub_vals)->subs(sub_vals))) << std::endl;
+  deallog << "Eval: g(x=1,y=2.2): "
+          << eval_double(*(g->subs(sub_vals)->subs(sub_vals))) << std::endl;
+  deallog << "Eval: h(x=1,y=2.2): "
+          << eval_double(*(h->subs(sub_vals)->subs(sub_vals))) << std::endl;
+
+  // Not yet implemented
+  deallog << "df_dx(x=1,y=2.2): " << *(df_dx->subs(sub_vals))
+          << std::endl; // ->subs(sub_vals)
+  deallog << "dg_dx(x=1,y=2.2): " << *(dg_dx->subs(sub_vals))
+          << std::endl; // ->subs(sub_vals)
+  deallog << "dh_dx(x=1,y=2.2): " << *(dh_dx->subs(sub_vals))
+          << std::endl; // ->subs(sub_vals)
+  // deallog << "Eval: df_dx(x=1,y=2.2): " <<
+  // eval_double(*(df_dx->subs(sub_vals)->subs(sub_vals))) << std::endl; deallog
+  // << "Eval: dg_dx(x=1,y=2.2): " <<
+  // eval_double(*(dg_dx->subs(sub_vals)->subs(sub_vals))) << std::endl; deallog
+  // << "Eval: dh_dx(x=1,y=2.2): " <<
+  // eval_double(*(dh_dx->subs(sub_vals)->subs(sub_vals))) << std::endl;
+}
+
+int
+main()
+{
+  initlog();
+
+  deallog.push("Integer");
+  test_number<int>();
+  deallog.pop();
+
+  deallog.push("Float");
+  test_number<float>();
+  deallog.pop();
+
+  deallog.push("Double");
+  test_number<double>();
+  deallog.pop();
+
+  // Not available yet
+  // SymEngine::SymEngineException: Invalid Format: Expected Integer or Rational
+  // deallog << "Complex double" << std::endl;
+  // test_number<std::complex<double>>();
+  // deallog.pop();
+
+  deallog << "OK" << std::endl;
+}
diff --git a/tests/symengine/basic_03.output b/tests/symengine/basic_03.output
new file mode 100644 (file)
index 0000000..d65a0e7
--- /dev/null
@@ -0,0 +1,62 @@
+
+DEAL:Integer::x: x
+DEAL:Integer::y: y
+DEAL:Integer::f: f(x)
+DEAL:Integer::g: g(x, y)
+DEAL:Integer::h: h(x + y, x*y)
+DEAL:Integer::df_dx: Derivative(f(x), x)
+DEAL:Integer::dg_dx: Derivative(g(x, y), x)
+DEAL:Integer::dh_dx: y*Subs(Derivative(h(x + y, _xi_2), _xi_2), (_xi_2), (x*y)) + Subs(Derivative(h(_xi_1, x*y), _xi_1), (_xi_1), (x + y))
+DEAL:Integer::df_dy: 0
+DEAL:Integer::dg_dy: Derivative(g(x, y), y)
+DEAL:Integer::dh_dy: x*Subs(Derivative(h(x + y, _xi_2), _xi_2), (_xi_2), (x*y)) + Subs(Derivative(h(_xi_1, x*y), _xi_1), (_xi_1), (x + y))
+DEAL:Integer::f(x=1,y=2.2): x**2
+DEAL:Integer::g(x=1,y=2.2): x + y
+DEAL:Integer::h(x=1,y=2.2): x - y
+DEAL:Integer::Eval: f(x=1,y=2.2): 1.00
+DEAL:Integer::Eval: g(x=1,y=2.2): 3.00
+DEAL:Integer::Eval: h(x=1,y=2.2): -1.00
+DEAL:Integer::df_dx(x=1,y=2.2): 2*x
+DEAL:Integer::dg_dx(x=1,y=2.2): 1
+DEAL:Integer::dh_dx(x=1,y=2.2): 2*Subs(Derivative(h(3, _xi_2), _xi_2), (_xi_2), (2)) + Subs(Derivative(h(_xi_1, 2), _xi_1), (_xi_1), (3))
+DEAL:Float::x: x
+DEAL:Float::y: y
+DEAL:Float::f: f(x)
+DEAL:Float::g: g(x, y)
+DEAL:Float::h: h(x + y, x*y)
+DEAL:Float::df_dx: Derivative(f(x), x)
+DEAL:Float::dg_dx: Derivative(g(x, y), x)
+DEAL:Float::dh_dx: y*Subs(Derivative(h(x + y, _xi_2), _xi_2), (_xi_2), (x*y)) + Subs(Derivative(h(_xi_1, x*y), _xi_1), (_xi_1), (x + y))
+DEAL:Float::df_dy: 0
+DEAL:Float::dg_dy: Derivative(g(x, y), y)
+DEAL:Float::dh_dy: x*Subs(Derivative(h(x + y, _xi_2), _xi_2), (_xi_2), (x*y)) + Subs(Derivative(h(_xi_1, x*y), _xi_1), (_xi_1), (x + y))
+DEAL:Float::f(x=1,y=2.2): x**2
+DEAL:Float::g(x=1,y=2.2): x + y
+DEAL:Float::h(x=1,y=2.2): x - y
+DEAL:Float::Eval: f(x=1,y=2.2): 1.00
+DEAL:Float::Eval: g(x=1,y=2.2): 3.20
+DEAL:Float::Eval: h(x=1,y=2.2): -1.20
+DEAL:Float::df_dx(x=1,y=2.2): 2*x
+DEAL:Float::dg_dx(x=1,y=2.2): 1
+DEAL:Float::dh_dx(x=1,y=2.2): 2.20000004768372*Subs(Derivative(h(3.20000004768372, _xi_2), _xi_2), (_xi_2), (2.20000004768372)) + Subs(Derivative(h(_xi_1, 2.20000004768372), _xi_1), (_xi_1), (3.20000004768372))
+DEAL:Double::x: x
+DEAL:Double::y: y
+DEAL:Double::f: f(x)
+DEAL:Double::g: g(x, y)
+DEAL:Double::h: h(x + y, x*y)
+DEAL:Double::df_dx: Derivative(f(x), x)
+DEAL:Double::dg_dx: Derivative(g(x, y), x)
+DEAL:Double::dh_dx: y*Subs(Derivative(h(x + y, _xi_2), _xi_2), (_xi_2), (x*y)) + Subs(Derivative(h(_xi_1, x*y), _xi_1), (_xi_1), (x + y))
+DEAL:Double::df_dy: 0
+DEAL:Double::dg_dy: Derivative(g(x, y), y)
+DEAL:Double::dh_dy: x*Subs(Derivative(h(x + y, _xi_2), _xi_2), (_xi_2), (x*y)) + Subs(Derivative(h(_xi_1, x*y), _xi_1), (_xi_1), (x + y))
+DEAL:Double::f(x=1,y=2.2): x**2
+DEAL:Double::g(x=1,y=2.2): x + y
+DEAL:Double::h(x=1,y=2.2): x - y
+DEAL:Double::Eval: f(x=1,y=2.2): 1.00
+DEAL:Double::Eval: g(x=1,y=2.2): 3.20
+DEAL:Double::Eval: h(x=1,y=2.2): -1.20
+DEAL:Double::df_dx(x=1,y=2.2): 2*x
+DEAL:Double::dg_dx(x=1,y=2.2): 1
+DEAL:Double::dh_dx(x=1,y=2.2): 2.2*Subs(Derivative(h(3.2, _xi_2), _xi_2), (_xi_2), (2.2)) + Subs(Derivative(h(_xi_1, 2.2), _xi_1), (_xi_1), (3.2))
+DEAL::OK
diff --git a/tests/symengine/basic_04.cc b/tests/symengine/basic_04.cc
new file mode 100644 (file)
index 0000000..9ffad42
--- /dev/null
@@ -0,0 +1,170 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2016 - 2018 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 at
+// the top level of the deal.II distribution.
+//
+// ---------------------------------------------------------------------
+
+
+// Check that SymEngine can do some optimsation using lambda functions
+
+// References:
+// https://github.com/symengine/symengine/blob/master/symengine/tests/eval/test_lambda_double.cpp
+
+#include "../tests.h"
+
+#pragma GCC diagnostic push
+#pragma GCC diagnostic ignored "-Wextra-semi"
+
+#include <symengine/add.h>
+#include <symengine/basic.h>
+#include <symengine/constants.h>
+#include <symengine/dict.h>
+#include <symengine/integer.h>
+#include <symengine/lambda_double.h>
+#include <symengine/llvm_double.h>
+#include <symengine/mul.h>
+#include <symengine/parser.h>
+#include <symengine/pow.h>
+#include <symengine/real_double.h>
+#include <symengine/symbol.h>
+
+#pragma GCC diagnostic pop
+
+#include <chrono>
+#include <iostream>
+#include <vector>
+
+using namespace dealii;
+namespace SE = SymEngine;
+
+int
+main(int argc, char *argv[])
+{
+  initlog();
+
+  const unsigned int n_runs = 1000;
+
+  SE::RCP<const SE::Symbol> t_11 = SE::symbol("t_11");
+  SE::RCP<const SE::Symbol> t_12 = SE::symbol("t_12");
+  SE::RCP<const SE::Symbol> s    = SE::symbol("s");
+  SE::RCP<const SE::Symbol> t_02 = SE::symbol("t_02");
+  SE::RCP<const SE::Symbol> v_0  = SE::symbol("v_0");
+  SE::RCP<const SE::Symbol> v_1  = SE::symbol("v_1");
+  SE::RCP<const SE::Symbol> v_2  = SE::symbol("v_2");
+  SE::RCP<const SE::Symbol> t_10 = SE::symbol("t_10");
+  SE::RCP<const SE::Symbol> t_00 = SE::symbol("t_00");
+  SE::RCP<const SE::Symbol> t_01 = SE::symbol("t_01");
+  SE::RCP<const SE::Symbol> t_20 = SE::symbol("t_20");
+  SE::RCP<const SE::Symbol> t_21 = SE::symbol("t_21");
+  SE::RCP<const SE::Symbol> t_22 = SE::symbol("t_22");
+
+  SE::vec_basic v = {
+    t_11, t_12, s, t_02, v_0, v_1, v_2, t_10, t_00, t_01, t_20, t_21, t_22};
+
+  SE::RCP<const SE::Basic> h =
+    SE::parse("2.*s**2.2*(v_0**2 + v_1**2 + v_2**2)**3*"
+              "(-1.*t_20*(1.*t_01*t_12 - 1.*t_02*t_11) "
+              "+ 1.*t_21*(1.*t_00*t_12 - 1.*t_02*t_10) "
+              "- 1.*t_22*(1.*t_00*t_11 - 1.*t_01*t_10))"
+              "*(1.*t_21*t_12 - 1.*t_22*t_11)");
+
+  SE::map_basic_basic dict;
+  std::vector<double> vals(v.size());
+  for (unsigned i = 0; i < v.size(); i++)
+    {
+      dict[v[i]] = SE::real_double(i);
+      vals[i]    = i;
+    }
+
+  double res = 0.0, res1 = 0.0, res2 = 0.0;
+
+
+  // Standard substitution
+
+  auto t1 = std::chrono::high_resolution_clock::now();
+  for (unsigned j = 0; j < n_runs; j++)
+    {
+      res += static_cast<const SE::RealDouble &>(*h->subs(dict)).i;
+    }
+  auto         t2 = std::chrono::high_resolution_clock::now();
+  const double diff_symm_subs =
+    std::chrono::duration_cast<std::chrono::microseconds>(t2 - t1).count();
+  std::cout
+    << "Subs " << n_runs << " calls                  :"
+    << std::chrono::duration_cast<std::chrono::microseconds>(t2 - t1).count()
+    << " us" << std::endl;
+
+  deallog << "Value (subs): " << res << std::endl;
+
+  // Optimisation via lambda functions
+
+  t1 = std::chrono::high_resolution_clock::now();
+  SE::LambdaRealDoubleVisitor l;
+  l.init(v, *h);
+  t2 = std::chrono::high_resolution_clock::now();
+  const double diff_lambda_real_double_setup =
+    std::chrono::duration_cast<std::chrono::microseconds>(t2 - t1).count();
+  std::cout
+    << "LambdaDoubleVisitor setup-time   :"
+    << std::chrono::duration_cast<std::chrono::microseconds>(t2 - t1).count()
+    << " us" << std::endl;
+
+  t1 = std::chrono::high_resolution_clock::now();
+  for (unsigned j = 0; j < n_runs; j++)
+    {
+      res1 += l.call(vals);
+    }
+  t2 = std::chrono::high_resolution_clock::now();
+  const double diff_lambda_real_double_call =
+    std::chrono::duration_cast<std::chrono::microseconds>(t2 - t1).count();
+  std::cout
+    << "LambdaDoubleVisitor run-time     :"
+    << std::chrono::duration_cast<std::chrono::microseconds>(t2 - t1).count()
+    << " us" << std::endl;
+
+  deallog << "Value (lambda): " << res1 << std::endl;
+
+  std::cout << "Timing ratio (subs/lambda_double): "
+            << (diff_symm_subs) /
+                 (diff_lambda_real_double_setup + diff_lambda_real_double_call)
+            << std::endl;
+
+#ifdef HAVE_SYMENGINE_LLVM
+  // Optimisation via LLVM JIT optimiser
+
+  t1 = std::chrono::high_resolution_clock::now();
+  SE::LLVMDoubleVisitor l2;
+  l2.init(v, *h);
+  t2 = std::chrono::high_resolution_clock::now();
+  std::cout
+    << "LLVMDoubleVisitor setup-time     :"
+    << std::chrono::duration_cast<std::chrono::microseconds>(t2 - t1).count()
+    << " us" << std::endl;
+
+  t1 = std::chrono::high_resolution_clock::now();
+  for (unsigned j = 0; j < n_runs; j++)
+    {
+      res2 += l2.call(vals);
+    }
+  t2 = std::chrono::high_resolution_clock::now();
+  std::cout
+    << "LLVMDoubleVisitor run-time       :"
+    << std::chrono::duration_cast<std::chrono::microseconds>(t2 - t1).count()
+    << " us : Value: " << res2 << std::endl;
+
+  deallog << "Value (llvm): " << res2 << std::endl;
+#endif
+
+  deallog << "OK" << std::endl;
+
+  return 0;
+}
diff --git a/tests/symengine/basic_04.with_symengine_with_llvm=false.output b/tests/symengine/basic_04.with_symengine_with_llvm=false.output
new file mode 100644 (file)
index 0000000..7eb32c0
--- /dev/null
@@ -0,0 +1,4 @@
+
+DEAL::Value (subs): 2.41358e+13
+DEAL::Value (lambda): 2.41358e+13
+DEAL::OK
diff --git a/tests/symengine/basic_04.with_symengine_with_llvm=true.output b/tests/symengine/basic_04.with_symengine_with_llvm=true.output
new file mode 100644 (file)
index 0000000..e230bb2
--- /dev/null
@@ -0,0 +1,5 @@
+
+DEAL::Value (subs): 2.41358e+13
+DEAL::Value (lambda): 2.41358e+13
+DEAL::Value (llvm): 2.41358e+13
+DEAL::OK
diff --git a/tests/symengine/basic_05.cc b/tests/symengine/basic_05.cc
new file mode 100644 (file)
index 0000000..41d59d1
--- /dev/null
@@ -0,0 +1,190 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2016 - 2018 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 at
+// the top level of the deal.II distribution.
+//
+// ---------------------------------------------------------------------
+
+// Check that SymEngine can do some optimsation on diffentiated expressions
+// using lambda functions
+// This is an extension of basic_04.cc
+
+// References:
+// https://github.com/symengine/symengine/blob/master/symengine/tests/eval/test_lambda_double.cpp
+
+#include "../tests.h"
+
+#pragma GCC diagnostic push
+#pragma GCC diagnostic ignored "-Wextra-semi"
+
+#include <symengine/add.h>
+#include <symengine/basic.h>
+#include <symengine/dict.h>
+#include <symengine/integer.h>
+#include <symengine/lambda_double.h>
+#include <symengine/llvm_double.h>
+#include <symengine/mul.h>
+#include <symengine/parser.h>
+#include <symengine/pow.h>
+#include <symengine/real_double.h>
+#include <symengine/symbol.h>
+
+#pragma GCC diagnostic pop
+
+#include <chrono>
+#include <iostream>
+#include <vector>
+
+using namespace dealii;
+namespace SE = SymEngine;
+
+int
+main(int argc, char *argv[])
+{
+  initlog();
+
+  const unsigned int n_runs = 1000;
+
+  SE::RCP<const SE::Symbol> t_11 = SE::symbol("t_11");
+  SE::RCP<const SE::Symbol> t_12 = SE::symbol("t_12");
+  SE::RCP<const SE::Symbol> s    = SE::symbol("s");
+  SE::RCP<const SE::Symbol> t_02 = SE::symbol("t_02");
+  SE::RCP<const SE::Symbol> v_0  = SE::symbol("v_0");
+  SE::RCP<const SE::Symbol> v_1  = SE::symbol("v_1");
+  SE::RCP<const SE::Symbol> v_2  = SE::symbol("v_2");
+  SE::RCP<const SE::Symbol> t_10 = SE::symbol("t_10");
+  SE::RCP<const SE::Symbol> t_00 = SE::symbol("t_00");
+  SE::RCP<const SE::Symbol> t_01 = SE::symbol("t_01");
+  SE::RCP<const SE::Symbol> t_20 = SE::symbol("t_20");
+  SE::RCP<const SE::Symbol> t_21 = SE::symbol("t_21");
+  SE::RCP<const SE::Symbol> t_22 = SE::symbol("t_22");
+
+  SE::vec_basic v = {
+    t_11, t_12, s, t_02, v_0, v_1, v_2, t_10, t_00, t_01, t_20, t_21, t_22};
+
+  SE::RCP<const SE::Basic> h =
+    SE::parse("2.*s**2.2*(v_0**2 + v_1**2 + v_2**2)**3*"
+              "(-1.*t_20*(1.*t_01*t_12 - 1.*t_02*t_11) "
+              "+ 1.*t_21*(1.*t_00*t_12 - 1.*t_02*t_10) "
+              "- 1.*t_22*(1.*t_00*t_11 - 1.*t_01*t_10))"
+              "*(1.*t_21*t_12 - 1.*t_22*t_11)");
+
+  SE::vec_basic diffs = {h->diff(t_20),
+                         h->diff(t_21),
+                         h->diff(t_22),
+                         h->diff(t_10),
+                         h->diff(t_11),
+                         h->diff(t_12),
+                         h->diff(t_00),
+                         h->diff(t_01),
+                         h->diff(t_02),
+                         h->diff(s)};
+
+  std::cout << *h << std::endl;
+
+  SE::map_basic_basic dict;
+  std::vector<double> vals(v.size());
+  for (unsigned i = 0; i < v.size(); i++)
+    {
+      dict[v[i]] = SE::real_double(i);
+      vals[i]    = i;
+    }
+
+  double res[10] = {}, res1[10] = {}, res2[10] = {};
+  double r = 0.0, r2 = 0.0, r3 = 0.0;
+
+  // Standard substitution
+
+  auto t1 = std::chrono::high_resolution_clock::now();
+  for (unsigned j = 0; j < n_runs; j++)
+    {
+      for (unsigned k = 0; k < diffs.size(); k++)
+        {
+          res[k] = static_cast<const SE::RealDouble &>(*diffs[k]->subs(dict)).i;
+          r += res[k];
+        }
+    }
+  auto         t2 = std::chrono::high_resolution_clock::now();
+  const double diff_symm_subs =
+    std::chrono::duration_cast<std::chrono::microseconds>(t2 - t1).count();
+  std::cout
+    << "Subs " << n_runs << " calls                  :"
+    << std::chrono::duration_cast<std::chrono::microseconds>(t2 - t1).count()
+    << " us" << std::endl;
+
+  deallog << "Value (subs): " << r << std::endl;
+
+  // Optimisation via lambda functions
+
+  t1 = std::chrono::high_resolution_clock::now();
+  SE::LambdaRealDoubleVisitor l;
+  l.init(v, diffs);
+  t2 = std::chrono::high_resolution_clock::now();
+  const double diff_lambda_real_double_setup =
+    std::chrono::duration_cast<std::chrono::microseconds>(t2 - t1).count();
+  std::cout
+    << "LambdaDoubleVisitor setup-time   :"
+    << std::chrono::duration_cast<std::chrono::microseconds>(t2 - t1).count()
+    << " us" << std::endl;
+
+  t1 = std::chrono::high_resolution_clock::now();
+  for (unsigned j = 0; j < n_runs; j++)
+    {
+      l.call(res1, vals.data());
+      for (unsigned k = 0; k < diffs.size(); k++)
+        {
+          r2 += res1[k];
+        }
+    }
+  t2 = std::chrono::high_resolution_clock::now();
+  const double diff_lambda_real_double_call =
+    std::chrono::duration_cast<std::chrono::microseconds>(t2 - t1).count();
+  std::cout
+    << "LambdaDoubleVisitor run-time     :"
+    << std::chrono::duration_cast<std::chrono::microseconds>(t2 - t1).count()
+    << " us" << std::endl;
+
+  deallog << "Value (lambda): " << r2 << std::endl;
+
+#ifdef HAVE_SYMENGINE_LLVM
+  // Optimisation via LLVM JIT optimiser
+
+  t1 = std::chrono::high_resolution_clock::now();
+  SE::LLVMDoubleVisitor l2;
+  l2.init(v, diffs);
+  t2 = std::chrono::high_resolution_clock::now();
+  std::cout
+    << "LLVMDoubleVisitor setup-time     :"
+    << std::chrono::duration_cast<std::chrono::microseconds>(t2 - t1).count()
+    << " us" << std::endl;
+
+  t1 = std::chrono::high_resolution_clock::now();
+  for (unsigned j = 0; j < 1000; j++)
+    {
+      l2.call(res2, vals.data());
+      for (unsigned k = 0; k < diffs.size(); k++)
+        {
+          r3 += res2[k];
+        }
+    }
+  t2 = std::chrono::high_resolution_clock::now();
+  std::cout
+    << "LLVMDoubleVisitor run-time       :"
+    << std::chrono::duration_cast<std::chrono::microseconds>(t2 - t1).count()
+    << " us : Value: " << r3 << std::endl;
+
+  deallog << "Value (llvm): " << r3 << std::endl;
+#endif
+
+  deallog << "OK" << std::endl;
+
+  return 0;
+}
diff --git a/tests/symengine/basic_05.with_symengine_with_llvm=false.output b/tests/symengine/basic_05.with_symengine_with_llvm=false.output
new file mode 100644 (file)
index 0000000..b68d7bd
--- /dev/null
@@ -0,0 +1,4 @@
+
+DEAL::Value (subs): 2.91338e+13
+DEAL::Value (lambda): 2.91338e+13
+DEAL::OK
diff --git a/tests/symengine/basic_05.with_symengine_with_llvm=true.output b/tests/symengine/basic_05.with_symengine_with_llvm=true.output
new file mode 100644 (file)
index 0000000..14de91c
--- /dev/null
@@ -0,0 +1,5 @@
+
+DEAL::Value (subs): 2.91338e+13
+DEAL::Value (lambda): 2.91338e+13
+DEAL::Value (llvm): 2.91338e+13
+DEAL::OK
diff --git a/tests/symengine/basic_06.cc b/tests/symengine/basic_06.cc
new file mode 100644 (file)
index 0000000..a664af2
--- /dev/null
@@ -0,0 +1,115 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2016 - 2018 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 at
+// the top level of the deal.II distribution.
+//
+// ---------------------------------------------------------------------
+
+// Check that SymEngine can do delayed substitution of explicitly dependent
+// symbolic functions and evaluate their derivatives as well.
+// This is inline with what is required to implement some material laws that
+// have internal variables
+
+#include "../tests.h"
+
+#pragma GCC diagnostic push
+#pragma GCC diagnostic ignored "-Wextra-semi"
+
+#include <symengine/add.h>
+#include <symengine/basic.h>
+#include <symengine/dict.h>
+#include <symengine/eval_double.h>
+#include <symengine/functions.h>
+#include <symengine/integer.h>
+#include <symengine/mul.h>
+#include <symengine/pow.h>
+#include <symengine/real_double.h>
+#include <symengine/symbol.h>
+
+#pragma GCC diagnostic pop
+
+#include <iostream>
+#include <vector>
+
+using namespace dealii;
+namespace SE = SymEngine;
+
+int
+main(int argc, char *argv[])
+{
+  initlog();
+
+  // 0. Define some independent variables
+  // Here we have two definitions of Q: One it is treated as a completely
+  // indepdent variable, and the other where its dependence on C is
+  // explicitly defined
+  SE::RCP<const SE::Symbol> C  = SE::symbol("C");  // Primary variable
+  SE::RCP<const SE::Symbol> Qi = SE::symbol("Qi"); // Some internal variable
+  SE::RCP<const SE::Basic>  Q =
+    SE::function_symbol("Q", C); // This sets up the dependence of Q on C
+
+  deallog << *C << ",  " << *Qi << ",  " << *Q << std::endl;
+
+  // 1. Define a function with Q being independent of C
+  SE::RCP<const SE::Basic> f_CQ_symb =
+    SE::mul(SE::real_double(0.5), SE::mul(Qi, SE::pow(C, SE::integer(2))));
+
+  deallog << "f_CQ_symb: " << *f_CQ_symb << std::endl;
+
+  // 2. Compute the partial derivative of f wrt C. This is a "partial
+  // derivative" as Q != Q(C).
+  SE::RCP<const SE::Basic> df_CQ_dC_symb = f_CQ_symb->diff(C);
+
+  deallog << "df_CQ_dC_symb: " << *df_CQ_dC_symb << std::endl;
+
+  // 3. Substitute Qi -> Q=Q(C) to now make SymEngine aware of Q's dependence
+  // on C.
+  SE::map_basic_basic int_var_dict;
+  int_var_dict[Qi] = Q;
+  SE::RCP<const SE::Basic> df_CQ_dC_symb_subs =
+    df_CQ_dC_symb->subs(int_var_dict);
+
+  deallog << "df_CQ_dC_symb_subs: " << *df_CQ_dC_symb_subs << std::endl;
+
+  // 4. Compute total derivative of df_dC wrt C.
+  SE::RCP<const SE::Basic> D2f_CQ_DC_dC_symb = df_CQ_dC_symb_subs->diff(C);
+
+  deallog << "D2f_CQ_DC_dC_symb: " << *D2f_CQ_DC_dC_symb << std::endl;
+
+  // 5. Perform some numerical substitutions
+  SE::map_basic_basic all_var_dict;
+  all_var_dict[C]  = SE::real_double(5);
+  all_var_dict[Qi] = SE::real_double(3);
+  all_var_dict[Q] =
+    SE::mul(SE::real_double(0.6),
+            C); // This explicitly defines the relationship between Q and C
+
+  SE::RCP<const SE::Basic> f_CQ_subs = f_CQ_symb->subs(all_var_dict);
+  deallog << "f_CQ_subs: " << *f_CQ_subs << std::endl;
+
+  SE::RCP<const SE::Basic> df_CQ_dC_subs = df_CQ_dC_symb->subs(all_var_dict);
+  deallog << "df_CQ_dC: subs: " << *df_CQ_dC_subs << std::endl
+          << "df_CQ_dC: eval: " << SE::eval_double(*df_CQ_dC_subs) << std::endl;
+
+  // This first substitution should presumably convert Q(C)->0.6C and thus
+  // dQ(C)_dC into 0.6
+  SE::RCP<const SE::Basic> D2f_CQ_DC_dC_subs_1 =
+    D2f_CQ_DC_dC_symb->subs(all_var_dict);
+  // This second substitution should fill out the remaining entries for C
+  SE::RCP<const SE::Basic> D2f_CQ_DC_dC_subs =
+    D2f_CQ_DC_dC_subs_1->subs(all_var_dict);
+  deallog << "D2f_CQ_DC_dC: subs 1: " << *D2f_CQ_DC_dC_subs_1 << std::endl
+          << "D2f_CQ_DC_dC: subs 2: " << *D2f_CQ_DC_dC_subs << std::endl
+          << "D2f_CQ_DC_dC: eval: " << SE::eval_double(*D2f_CQ_DC_dC_subs)
+          << std::endl;
+
+  return 0;
+}
diff --git a/tests/symengine/basic_06.output b/tests/symengine/basic_06.output
new file mode 100644 (file)
index 0000000..ad5bc7b
--- /dev/null
@@ -0,0 +1,12 @@
+
+DEAL::C,  Qi,  Q(C)
+DEAL::f_CQ_symb: 0.5*C**2*Qi
+DEAL::df_CQ_dC_symb: 1.0*C*Qi
+DEAL::df_CQ_dC_symb_subs: 1.0*C*Q(C)
+DEAL::D2f_CQ_DC_dC_symb: 1.0*C*Derivative(Q(C), C) + 1.0*Q(C)
+DEAL::f_CQ_subs: 37.5
+DEAL::df_CQ_dC: subs: 15.
+DEAL::df_CQ_dC: eval: 15.0000
+DEAL::D2f_CQ_DC_dC: subs 1: 3.0 + 0.6*C
+DEAL::D2f_CQ_DC_dC: subs 2: 6.0
+DEAL::D2f_CQ_DC_dC: eval: 6.00000
diff --git a/tests/symengine/basic_07.cc b/tests/symengine/basic_07.cc
new file mode 100644 (file)
index 0000000..2963db1
--- /dev/null
@@ -0,0 +1,123 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2016 - 2018 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 at
+// the top level of the deal.II distribution.
+//
+// ---------------------------------------------------------------------
+
+// Check that SymEngine can do delayed substitution of implicitly dependent
+// symbolic functions and evaluate their derivatives as well.
+// This is inline with what is required to implement some material laws that
+// have internal variables
+
+#include "../tests.h"
+
+#pragma GCC diagnostic push
+#pragma GCC diagnostic ignored "-Wextra-semi"
+
+#include <symengine/add.h>
+#include <symengine/basic.h>
+#include <symengine/dict.h>
+#include <symengine/eval_double.h>
+#include <symengine/functions.h>
+#include <symengine/integer.h>
+#include <symengine/mul.h>
+#include <symengine/pow.h>
+#include <symengine/real_double.h>
+#include <symengine/symbol.h>
+
+#pragma GCC diagnostic pop
+
+#include <iostream>
+#include <vector>
+
+using namespace dealii;
+namespace SE = SymEngine;
+
+int
+main(int argc, char *argv[])
+{
+  initlog();
+
+  // 0. Define some independent variables
+  // Here we have two definitions of Q: One it is treated as a completely
+  // independent variable, and the other where its dependence on C is
+  // implicitly defined
+  // That is to say, we'd have to compute Q=Q(C) using a nonlinear solution
+  // scheme, as might be the case for a viscoelastic material with a
+  // nonlinear evolution law.
+  SE::RCP<const SE::Symbol> C  = SE::symbol("C");  // Primary variable
+  SE::RCP<const SE::Symbol> Qi = SE::symbol("Qi"); // Some internal variable
+  SE::RCP<const SE::Basic>  Q =
+    SE::function_symbol("Q", C); // This sets up the dependence of Q on C
+
+  deallog << *C << ",  " << *Qi << ",  " << *Q << std::endl;
+
+  // 1. Define a function with Q being independent of C
+  SE::RCP<const SE::Basic> f_CQ_symb =
+    SE::mul(SE::real_double(0.5), SE::mul(Qi, SE::pow(C, SE::integer(2))));
+
+  deallog << "f_CQ_symb: " << *f_CQ_symb << std::endl;
+
+  // 2. Compute the partial derivative of f wrt C. This is a "partial
+  // derivative" as Q != Q(C).
+  SE::RCP<const SE::Basic> df_CQ_dC_symb = f_CQ_symb->diff(C);
+
+  deallog << "df_CQ_dC_symb: " << *df_CQ_dC_symb << std::endl;
+
+  // 3. Substitute Qi -> Q=Q(C) to now make SymEngine aware of Q's dependence
+  // on C.
+  SE::map_basic_basic int_var_dict;
+  int_var_dict[Qi] = Q;
+  SE::RCP<const SE::Basic> df_CQ_dC_symb_subs =
+    df_CQ_dC_symb->subs(int_var_dict);
+
+  deallog << "df_CQ_dC_symb_subs: " << *df_CQ_dC_symb_subs << std::endl;
+
+  // 4. Compute total derivative of df_dC wrt C.
+  SE::RCP<const SE::Basic> D2f_CQ_DC_dC_symb = df_CQ_dC_symb_subs->diff(C);
+
+  deallog << "D2f_CQ_DC_dC_symb: " << *D2f_CQ_DC_dC_symb << std::endl;
+
+  // 5. Perform some numerical substitutions
+  // Lets assume Q = 0.6*C*C
+  // Then dQ_dC = 1.2*C
+  SE::map_basic_basic all_var_dict;
+  all_var_dict[C] = SE::real_double(5);
+  all_var_dict[Qi] =
+    SE::real_double(0.6 * 5 * 5); // This would be the numerical end result of
+                                  // an internal nonlinear solver for Q
+  all_var_dict[Q] = Qi;           // These are numerically the same thing
+  all_var_dict[Q->diff(C)] =
+    SE::mul(SE::real_double(1.2),
+            C); // This resolves the implicit relationship between Q and C
+
+  SE::RCP<const SE::Basic> f_CQ_subs = f_CQ_symb->subs(all_var_dict);
+  deallog << "f_CQ_subs: " << *f_CQ_subs << std::endl;
+
+  SE::RCP<const SE::Basic> df_CQ_dC_subs = df_CQ_dC_symb->subs(all_var_dict);
+  deallog << "df_CQ_dC: subs: " << *df_CQ_dC_subs << std::endl
+          << "df_CQ_dC: eval: " << SE::eval_double(*df_CQ_dC_subs) << std::endl;
+
+  // This first substitution should presumably convert Q(C)->0.6C and thus
+  // dQ(C)_dC into 0.6
+  SE::RCP<const SE::Basic> D2f_CQ_DC_dC_subs_1 =
+    D2f_CQ_DC_dC_symb->subs(all_var_dict);
+  // This second substitution should fill out the remaining entries for C
+  SE::RCP<const SE::Basic> D2f_CQ_DC_dC_subs =
+    D2f_CQ_DC_dC_subs_1->subs(all_var_dict);
+  deallog << "D2f_CQ_DC_dC: subs 1: " << *D2f_CQ_DC_dC_subs_1 << std::endl
+          << "D2f_CQ_DC_dC: subs 2: " << *D2f_CQ_DC_dC_subs << std::endl
+          << "D2f_CQ_DC_dC: eval: " << SE::eval_double(*D2f_CQ_DC_dC_subs)
+          << std::endl;
+
+  return 0;
+}
diff --git a/tests/symengine/basic_07.output b/tests/symengine/basic_07.output
new file mode 100644 (file)
index 0000000..6b1e0b5
--- /dev/null
@@ -0,0 +1,12 @@
+
+DEAL::C,  Qi,  Q(C)
+DEAL::f_CQ_symb: 0.5*C**2*Qi
+DEAL::df_CQ_dC_symb: 1.0*C*Qi
+DEAL::df_CQ_dC_symb_subs: 1.0*C*Q(C)
+DEAL::D2f_CQ_DC_dC_symb: 1.0*C*Derivative(Q(C), C) + 1.0*Q(C)
+DEAL::f_CQ_subs: 187.5
+DEAL::df_CQ_dC: subs: 75.
+DEAL::df_CQ_dC: eval: 75.0000
+DEAL::D2f_CQ_DC_dC: subs 1: 6.0*C + 1.0*Qi
+DEAL::D2f_CQ_DC_dC: subs 2: 45.0
+DEAL::D2f_CQ_DC_dC: eval: 45.0000

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.