From: Jean-Paul Pelteret Date: Thu, 14 Mar 2019 17:27:51 +0000 (+0100) Subject: Add some initial tests for SymEngine X-Git-Tag: v9.1.0-rc1~271^2~3 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=8f713e8f0f889cea3269404fb65742ddefd17d3a;p=dealii.git Add some initial tests for SymEngine --- diff --git a/tests/symengine/CMakeLists.txt b/tests/symengine/CMakeLists.txt new file mode 100644 index 0000000000..46145a0466 --- /dev/null +++ b/tests/symengine/CMakeLists.txt @@ -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 index 0000000000..1022f6aa9e --- /dev/null +++ b/tests/symengine/basic_01.cc @@ -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 +#include +#include +#include +#include +#include +#pragma clang diagnostic pop + +#include +#include + +using namespace dealii; +namespace SE = SymEngine; + +SE::RCP +make_symengine_rcp(const int &val) +{ + return SE::integer(val); +} + +SE::RCP +make_symengine_rcp(const float &val) +{ + return SE::real_double(val); +} + +SE::RCP +make_symengine_rcp(const double &val) +{ + return SE::real_double(val); +} + +template +SE::RCP +make_symengine_rcp(const std::complex &val) +{ + // Build complex from two SymEngine numbers + return SE::Complex::from_two_nums(*make_symengine_rcp(val.real()), + *make_symengine_rcp(val.imag())); +} + +template +void +test_number() +{ + SE::RCP a = make_symengine_rcp(NumberType(4.2)); + deallog << "a: " << *a << std::endl; + + SE::RCP b(make_symengine_rcp(NumberType(2.1))); + deallog << "b: " << *b << std::endl; + + SE::RCP 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(); + deallog.pop(); + + deallog.push("Float"); + test_number(); + deallog.pop(); + + deallog.push("Double"); + test_number(); + deallog.pop(); + + // Not available yet + // SymEngine::SymEngineException: Invalid Format: Expected Integer or Rational + // deallog << "Complex double" << std::endl; + // test_number>(); + // deallog.pop(); + + deallog << "OK" << std::endl; +} diff --git a/tests/symengine/basic_01.output b/tests/symengine/basic_01.output new file mode 100644 index 0000000000..f6307bc99a --- /dev/null +++ b/tests/symengine/basic_01.output @@ -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 index 0000000000..7284e26f9c --- /dev/null +++ b/tests/symengine/basic_02.cc @@ -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 +#include +#include +#include +#include +#include +#include +#include +#include +#pragma clang diagnostic pop + +#include +#include + +using namespace dealii; +namespace SE = SymEngine; + +SE::RCP +make_symengine_rcp(const int &val) +{ + return SE::integer(val); +} + +SE::RCP +make_symengine_rcp(const float &val) +{ + return SE::real_double(val); +} + +SE::RCP +make_symengine_rcp(const double &val) +{ + return SE::real_double(val); +} + +template +SE::RCP +make_symengine_rcp(const std::complex &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 +make_symengine_rcp(const std::string &name) +{ + return SE::symbol(name); +} + +template +void +test_number() +{ + SE::RCP a = make_symengine_rcp(NumberType(3.1)); + SE::RCP b = make_symengine_rcp(NumberType(7.5)); + deallog << "a: " << *a << std::endl; + deallog << "b: " << *b << std::endl; + + SE::RCP x(make_symengine_rcp("x")); + SE::RCP y(make_symengine_rcp("y")); + deallog << "x: " << *x << std::endl; + deallog << "y: " << *y << std::endl; + + SE::RCP 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 dc_dx = c->diff(x); + // SE::RCP dc_dy = diff(c,SE::implicit_cast &>(y)); + SE::RCP dc_dy = + diff(c, SE::rcp_static_cast(y)); + SE::RCP 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(); + deallog.pop(); + + deallog.push("Float"); + test_number(); + deallog.pop(); + + deallog.push("Double"); + test_number(); + deallog.pop(); + + // Not available yet + // SymEngine::SymEngineException: Invalid Format: Expected Integer or Rational + // deallog << "Complex double" << std::endl; + // test_number>(); + // deallog.pop(); + + deallog << "OK" << std::endl; +} diff --git a/tests/symengine/basic_02.output b/tests/symengine/basic_02.output new file mode 100644 index 0000000000..a023eec9d8 --- /dev/null +++ b/tests/symengine/basic_02.output @@ -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 index 0000000000..3a11eb2be9 --- /dev/null +++ b/tests/symengine/basic_03.cc @@ -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 +#include +#include +#include +#include +#include +#include +#include +#include +#pragma clang diagnostic pop + +#include +#include + +using namespace dealii; +namespace SE = SymEngine; + +SE::RCP +make_symengine_rcp(const int &val) +{ + return SE::integer(val); +} + +SE::RCP +make_symengine_rcp(const float &val) +{ + return SE::real_double(val); +} + +SE::RCP +make_symengine_rcp(const double &val) +{ + return SE::real_double(val); +} + +template +SE::RCP +make_symengine_rcp(const std::complex &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 +make_symengine_rcp(const std::string &name) +{ + return SE::symbol(name); +} + +SE::RCP +make_symengine_rcp(const std::string & name, + const std::vector> &args) +{ + return SE::function_symbol(name, args); +} + +SE::RCP +make_symengine_rcp(const std::string &name, const SE::RCP &arg) +{ + return SE::function_symbol(name, arg); +} + +template +void +test_number() +{ + // Modified snippet from SymEngine test + // { + // SE::RCP x = SE::symbol("x"); + // SE::RCP _xi_1 = SE::symbol("_xi_1"); + // SE::RCP f = SE::function_symbol("f", x); + // SE::RCP 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 x(make_symengine_rcp("x")); + SE::RCP y(make_symengine_rcp("y")); + deallog << "x: " << *x << std::endl; + deallog << "y: " << *y << std::endl; + + // Function symbols + SE::RCP f(make_symengine_rcp("f", x)); + SE::RCP g(make_symengine_rcp("g", {x, y})); + SE::RCP 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 df_dx = f->diff(x); + SE::RCP dg_dx = g->diff(x); + SE::RCP 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 df_dy = + f->diff(SE::rcp_static_cast(y)); + SE::RCP dg_dy = + g->diff(SE::rcp_static_cast(y)); + SE::RCP dh_dy = + h->diff(SE::rcp_static_cast(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(); + deallog.pop(); + + deallog.push("Float"); + test_number(); + deallog.pop(); + + deallog.push("Double"); + test_number(); + deallog.pop(); + + // Not available yet + // SymEngine::SymEngineException: Invalid Format: Expected Integer or Rational + // deallog << "Complex double" << std::endl; + // test_number>(); + // deallog.pop(); + + deallog << "OK" << std::endl; +} diff --git a/tests/symengine/basic_03.output b/tests/symengine/basic_03.output new file mode 100644 index 0000000000..d65a0e705c --- /dev/null +++ b/tests/symengine/basic_03.output @@ -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 index 0000000000..9ffad422cf --- /dev/null +++ b/tests/symengine/basic_04.cc @@ -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 +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#pragma GCC diagnostic pop + +#include +#include +#include + +using namespace dealii; +namespace SE = SymEngine; + +int +main(int argc, char *argv[]) +{ + initlog(); + + const unsigned int n_runs = 1000; + + SE::RCP t_11 = SE::symbol("t_11"); + SE::RCP t_12 = SE::symbol("t_12"); + SE::RCP s = SE::symbol("s"); + SE::RCP t_02 = SE::symbol("t_02"); + SE::RCP v_0 = SE::symbol("v_0"); + SE::RCP v_1 = SE::symbol("v_1"); + SE::RCP v_2 = SE::symbol("v_2"); + SE::RCP t_10 = SE::symbol("t_10"); + SE::RCP t_00 = SE::symbol("t_00"); + SE::RCP t_01 = SE::symbol("t_01"); + SE::RCP t_20 = SE::symbol("t_20"); + SE::RCP t_21 = SE::symbol("t_21"); + SE::RCP 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 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 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(*h->subs(dict)).i; + } + auto t2 = std::chrono::high_resolution_clock::now(); + const double diff_symm_subs = + std::chrono::duration_cast(t2 - t1).count(); + std::cout + << "Subs " << n_runs << " calls :" + << std::chrono::duration_cast(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(t2 - t1).count(); + std::cout + << "LambdaDoubleVisitor setup-time :" + << std::chrono::duration_cast(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(t2 - t1).count(); + std::cout + << "LambdaDoubleVisitor run-time :" + << std::chrono::duration_cast(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(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(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 index 0000000000..7eb32c0d55 --- /dev/null +++ b/tests/symengine/basic_04.with_symengine_with_llvm=false.output @@ -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 index 0000000000..e230bb2c34 --- /dev/null +++ b/tests/symengine/basic_04.with_symengine_with_llvm=true.output @@ -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 index 0000000000..41d59d174d --- /dev/null +++ b/tests/symengine/basic_05.cc @@ -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 +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#pragma GCC diagnostic pop + +#include +#include +#include + +using namespace dealii; +namespace SE = SymEngine; + +int +main(int argc, char *argv[]) +{ + initlog(); + + const unsigned int n_runs = 1000; + + SE::RCP t_11 = SE::symbol("t_11"); + SE::RCP t_12 = SE::symbol("t_12"); + SE::RCP s = SE::symbol("s"); + SE::RCP t_02 = SE::symbol("t_02"); + SE::RCP v_0 = SE::symbol("v_0"); + SE::RCP v_1 = SE::symbol("v_1"); + SE::RCP v_2 = SE::symbol("v_2"); + SE::RCP t_10 = SE::symbol("t_10"); + SE::RCP t_00 = SE::symbol("t_00"); + SE::RCP t_01 = SE::symbol("t_01"); + SE::RCP t_20 = SE::symbol("t_20"); + SE::RCP t_21 = SE::symbol("t_21"); + SE::RCP 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 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 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(*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(t2 - t1).count(); + std::cout + << "Subs " << n_runs << " calls :" + << std::chrono::duration_cast(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(t2 - t1).count(); + std::cout + << "LambdaDoubleVisitor setup-time :" + << std::chrono::duration_cast(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(t2 - t1).count(); + std::cout + << "LambdaDoubleVisitor run-time :" + << std::chrono::duration_cast(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(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(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 index 0000000000..b68d7bd3b6 --- /dev/null +++ b/tests/symengine/basic_05.with_symengine_with_llvm=false.output @@ -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 index 0000000000..14de91cc39 --- /dev/null +++ b/tests/symengine/basic_05.with_symengine_with_llvm=true.output @@ -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 index 0000000000..a664af2b46 --- /dev/null +++ b/tests/symengine/basic_06.cc @@ -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 +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#pragma GCC diagnostic pop + +#include +#include + +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 C = SE::symbol("C"); // Primary variable + SE::RCP Qi = SE::symbol("Qi"); // Some internal variable + SE::RCP 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 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 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 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 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 f_CQ_subs = f_CQ_symb->subs(all_var_dict); + deallog << "f_CQ_subs: " << *f_CQ_subs << std::endl; + + SE::RCP 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 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 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 index 0000000000..ad5bc7ba11 --- /dev/null +++ b/tests/symengine/basic_06.output @@ -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 index 0000000000..2963db1902 --- /dev/null +++ b/tests/symengine/basic_07.cc @@ -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 +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#pragma GCC diagnostic pop + +#include +#include + +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 C = SE::symbol("C"); // Primary variable + SE::RCP Qi = SE::symbol("Qi"); // Some internal variable + SE::RCP 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 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 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 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 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 f_CQ_subs = f_CQ_symb->subs(all_var_dict); + deallog << "f_CQ_subs: " << *f_CQ_subs << std::endl; + + SE::RCP 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 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 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 index 0000000000..6b1e0b525c --- /dev/null +++ b/tests/symengine/basic_07.output @@ -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