From 7d20d6a798f620243592b2a124115c241500911c Mon Sep 17 00:00:00 2001 From: Jean-Paul Pelteret Date: Mon, 22 Apr 2019 15:57:20 +0200 Subject: [PATCH] Add a test for SD tensor differentiation --- .../symengine_tensor_operations_03.cc | 783 ++++++++++++++++++ .../symengine_tensor_operations_03.output | 4 + 2 files changed, 787 insertions(+) create mode 100644 tests/symengine/symengine_tensor_operations_03.cc create mode 100644 tests/symengine/symengine_tensor_operations_03.output diff --git a/tests/symengine/symengine_tensor_operations_03.cc b/tests/symengine/symengine_tensor_operations_03.cc new file mode 100644 index 0000000000..76cc90acdc --- /dev/null +++ b/tests/symengine/symengine_tensor_operations_03.cc @@ -0,0 +1,783 @@ +// --------------------------------------------------------------------- +// +// 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 at +// the top level of the deal.II distribution. +// +// --------------------------------------------------------------------- + + +// Test for symbolic tensor differentiation: +// Coupled double-differentiation example using SD-SD +// This test relates to ad_common_tests/helper_scalar_coupled_4_components_01.h + +#include +#include +#include +#include +#include +#include +#include + +#include + +#include + +#include +#include + +#include +#include + +#include "../tests.h" + +using namespace dealii; +namespace SD = Differentiation::SD; + +// Function and its derivatives +template +struct CoupledFunction +{ + static NumberType + det_t(const Tensor<2, dim, NumberType> &t) + { + return determinant(t); + } + + static Tensor<2, dim, NumberType> + ddet_t_dt(const Tensor<2, dim, NumberType> &t) + { + return det_t(t) * transpose(invert(t)); + } + + static Tensor<4, dim, NumberType> + d2det_t_dt_dt(const Tensor<2, dim, NumberType> &t) + { + const Tensor<2, dim, NumberType> t_inv = invert(t); + Tensor<4, dim, NumberType> dt_inv_trans_dt; + // https://en.wikiversity.org/wiki/Introduction_to_Elasticity/Tensors#Derivative_of_the_determinant_of_a_tensor + for (unsigned int i = 0; i < dim; ++i) + for (unsigned int j = 0; j < dim; ++j) + for (unsigned int k = 0; k < dim; ++k) + for (unsigned int l = 0; l < dim; ++l) + dt_inv_trans_dt[i][j][k][l] = -t_inv[l][i] * t_inv[j][k]; + + return det_t(t) * outer_product(transpose(t_inv), transpose(t_inv)) + + det_t(t) * dt_inv_trans_dt; + } + + static NumberType + det_t(const SymmetricTensor<2, dim, NumberType> &t) + { + return determinant(t); + } + + static SymmetricTensor<2, dim, NumberType> + ddet_t_dt(const SymmetricTensor<2, dim, NumberType> &t) + { + const SymmetricTensor<2, dim, NumberType> t_inv = invert(t); + return det_t(t) * t_inv; + } + + static SymmetricTensor<4, dim, NumberType> + d2det_t_dt_dt(const SymmetricTensor<2, dim, NumberType> &t) + { + const SymmetricTensor<2, dim, NumberType> t_inv = invert(t); + SymmetricTensor<4, dim, NumberType> dt_inv_dt; + // https://en.wikiversity.org/wiki/Introduction_to_Elasticity/Tensors#Derivative_of_the_determinant_of_a_tensor + for (unsigned int i = 0; i < dim; ++i) + for (unsigned int j = i; j < dim; ++j) + for (unsigned int k = 0; k < dim; ++k) + for (unsigned int l = k; l < dim; ++l) + dt_inv_dt[i][j][k][l] = + -0.5 * (t_inv[i][k] * t_inv[j][l] + t_inv[i][l] * t_inv[j][k]); + + return det_t(t) * outer_product(t_inv, t_inv) + det_t(t) * dt_inv_dt; + } + + static NumberType + v_squ(const Tensor<1, dim, NumberType> &v) + { + return v * v; + } + + static Tensor<1, dim, NumberType> + dv_squ_dv(const Tensor<1, dim, NumberType> &v) + { + return 2.0 * v; + } + + static Tensor<2, dim, NumberType> + d2v_squ_dv_dv(const Tensor<1, dim, NumberType> &v) + { + static const Tensor<2, dim, NumberType> I( + unit_symmetric_tensor()); + return 2.0 * I; + } + + // -------------------------------------------------------------------------- + + // Here we define a pseudo-energy function that we wish to differentiate twice + // with respect to all variables + + static const double sf; + + static NumberType + psi(const SymmetricTensor<2, dim, NumberType> &st, + const Tensor<2, dim, NumberType> & t, + const Tensor<1, dim, NumberType> & v, + const NumberType & s) + { + return std::pow(det_t(st), 2) * std::pow(det_t(t), 2) * + std::pow(v_squ(v), 3) * std::pow(s, sf); + }; + + // -------------------------------------------------------------------------- + + // Here's the hand-calculated expansion of this energy functions first and + // second derivatives. Lots of fun... + + static SymmetricTensor<2, dim, NumberType> + dpsi_dst(const SymmetricTensor<2, dim, NumberType> &st, + const Tensor<2, dim, NumberType> & t, + const Tensor<1, dim, NumberType> & v, + const NumberType & s) + { + return 2.0 * std::pow(det_t(st), 1) * ddet_t_dt(st) * + std::pow(det_t(t), 2) * std::pow(v_squ(v), 3) * std::pow(s, sf); + }; + + static Tensor<2, dim, NumberType> + dpsi_dt(const SymmetricTensor<2, dim, NumberType> &st, + const Tensor<2, dim, NumberType> & t, + const Tensor<1, dim, NumberType> & v, + const NumberType & s) + { + return 2.0 * std::pow(det_t(st), 2) * std::pow(det_t(t), 1) * ddet_t_dt(t) * + std::pow(v_squ(v), 3) * std::pow(s, sf); + }; + + static Tensor<1, dim, NumberType> + dpsi_dv(const SymmetricTensor<2, dim, NumberType> &st, + const Tensor<2, dim, NumberType> & t, + const Tensor<1, dim, NumberType> & v, + const NumberType & s) + { + return std::pow(det_t(st), 2) * std::pow(det_t(t), 2) * 3.0 * + std::pow(v_squ(v), 2) * dv_squ_dv(v) * std::pow(s, sf); + }; + + static NumberType + dpsi_ds(const SymmetricTensor<2, dim, NumberType> &st, + const Tensor<2, dim, NumberType> & t, + const Tensor<1, dim, NumberType> & v, + const NumberType & s) + { + return std::pow(det_t(st), 2) * std::pow(det_t(t), 2) * + std::pow(v_squ(v), 3) * sf * std::pow(s, sf - 1.0); + }; + + static SymmetricTensor<4, dim, NumberType> + d2psi_dst_x_dst(const SymmetricTensor<2, dim, NumberType> &st, + const Tensor<2, dim, NumberType> & t, + const Tensor<1, dim, NumberType> & v, + const NumberType & s) + { + return 2.0 * std::pow(det_t(t), 2) * std::pow(v_squ(v), 3) * + (std::pow(det_t(st), 0) * + outer_product(ddet_t_dt(st), ddet_t_dt(st)) + + std::pow(det_t(st), 1) * d2det_t_dt_dt(st)) * + std::pow(s, sf); + }; + + static Tensor<4, dim, NumberType> + d2psi_dst_x_dt(const SymmetricTensor<2, dim, NumberType> &st, + const Tensor<2, dim, NumberType> & t, + const Tensor<1, dim, NumberType> & v, + const NumberType & s) + { + return 4.0 * std::pow(det_t(st), 1) * std::pow(det_t(t), 1) * + outer_product(Tensor<2, dim, NumberType>(ddet_t_dt(st)), + ddet_t_dt(t)) * + std::pow(v_squ(v), 3) * std::pow(s, sf); + }; + + static Tensor<3, dim, NumberType> + d2psi_dst_x_dv(const SymmetricTensor<2, dim, NumberType> &st, + const Tensor<2, dim, NumberType> & t, + const Tensor<1, dim, NumberType> & v, + const NumberType & s) + { + return 2.0 * std::pow(det_t(t), 2) * std::pow(det_t(st), 1) * 3.0 * + std::pow(v_squ(v), 2) * + outer_product(Tensor<2, dim, NumberType>(ddet_t_dt(st)), + dv_squ_dv(v)) * + std::pow(s, sf); + }; + + static SymmetricTensor<2, dim, NumberType> + d2psi_dst_x_ds(const SymmetricTensor<2, dim, NumberType> &st, + const Tensor<2, dim, NumberType> & t, + const Tensor<1, dim, NumberType> & v, + const NumberType & s) + { + return 2.0 * std::pow(det_t(st), 1) * ddet_t_dt(st) * + std::pow(det_t(t), 2) * std::pow(v_squ(v), 3) * sf * + std::pow(s, sf - 1.0); + }; + + static Tensor<4, dim, NumberType> + d2psi_dt_x_dst(const SymmetricTensor<2, dim, NumberType> &st, + const Tensor<2, dim, NumberType> & t, + const Tensor<1, dim, NumberType> & v, + const NumberType & s) + { + return 4.0 * std::pow(det_t(st), 1) * std::pow(det_t(t), 1) * + outer_product(ddet_t_dt(t), + Tensor<2, dim, NumberType>(ddet_t_dt(st))) * + std::pow(v_squ(v), 3) * std::pow(s, sf); + }; + + static Tensor<4, dim, NumberType> + d2psi_dt_x_dt(const SymmetricTensor<2, dim, NumberType> &st, + const Tensor<2, dim, NumberType> & t, + const Tensor<1, dim, NumberType> & v, + const NumberType & s) + { + return 2.0 * std::pow(det_t(st), 2) * std::pow(v_squ(v), 3) * + (pow(det_t(t), 0) * outer_product(ddet_t_dt(t), ddet_t_dt(t)) + + std::pow(det_t(t), 1) * d2det_t_dt_dt(t)) * + std::pow(s, sf); + }; + + static Tensor<3, dim, NumberType> + d2psi_dt_x_dv(const SymmetricTensor<2, dim, NumberType> &st, + const Tensor<2, dim, NumberType> & t, + const Tensor<1, dim, NumberType> & v, + const NumberType & s) + { + return 2.0 * std::pow(det_t(st), 2) * std::pow(det_t(t), 1) * 3.0 * + std::pow(v_squ(v), 2) * outer_product(ddet_t_dt(t), dv_squ_dv(v)) * + std::pow(s, sf); + }; + + static Tensor<2, dim, NumberType> + d2psi_dt_x_ds(const SymmetricTensor<2, dim, NumberType> &st, + const Tensor<2, dim, NumberType> & t, + const Tensor<1, dim, NumberType> & v, + const NumberType & s) + { + return 2.0 * std::pow(det_t(st), 2) * std::pow(det_t(t), 1) * ddet_t_dt(t) * + std::pow(v_squ(v), 3) * sf * std::pow(s, sf - 1.0); + }; + + static Tensor<3, dim, NumberType> + d2psi_dv_x_dst(const SymmetricTensor<2, dim, NumberType> &st, + const Tensor<2, dim, NumberType> & t, + const Tensor<1, dim, NumberType> & v, + const NumberType & s) + { + return 2.0 * std::pow(det_t(t), 2) * std::pow(det_t(st), 1) * 3.0 * + std::pow(v_squ(v), 2) * + outer_product(dv_squ_dv(v), + Tensor<2, dim, NumberType>(ddet_t_dt(st))) * + std::pow(s, sf); + }; + + static Tensor<3, dim, NumberType> + d2psi_dv_x_dt(const SymmetricTensor<2, dim, NumberType> &st, + const Tensor<2, dim, NumberType> & t, + const Tensor<1, dim, NumberType> & v, + const NumberType & s) + { + return 2.0 * std::pow(det_t(st), 2) * std::pow(det_t(t), 1) * 3.0 * + std::pow(v_squ(v), 2) * outer_product(dv_squ_dv(v), ddet_t_dt(t)) * + std::pow(s, sf); + }; + + static Tensor<2, dim, NumberType> + d2psi_dv_x_dv(const SymmetricTensor<2, dim, NumberType> &st, + const Tensor<2, dim, NumberType> & t, + const Tensor<1, dim, NumberType> & v, + const NumberType & s) + { + return std::pow(det_t(st), 2) * std::pow(det_t(t), 2) * 3.0 * + (2.0 * std::pow(v_squ(v), 1) * + outer_product(dv_squ_dv(v), dv_squ_dv(v)) + + std::pow(v_squ(v), 2) * d2v_squ_dv_dv(v)) * + std::pow(s, sf); + }; + + static Tensor<1, dim, NumberType> + d2psi_dv_x_ds(const SymmetricTensor<2, dim, NumberType> &st, + const Tensor<2, dim, NumberType> & t, + const Tensor<1, dim, NumberType> & v, + const NumberType & s) + { + return std::pow(det_t(st), 2) * std::pow(det_t(t), 2) * 3.0 * + std::pow(v_squ(v), 2) * dv_squ_dv(v) * sf * std::pow(s, sf - 1.0); + }; + + static SymmetricTensor<2, dim, NumberType> + d2psi_ds_x_dst(const SymmetricTensor<2, dim, NumberType> &st, + const Tensor<2, dim, NumberType> & t, + const Tensor<1, dim, NumberType> & v, + const NumberType & s) + { + return 2.0 * ddet_t_dt(st) * std::pow(det_t(t), 2) * + std::pow(det_t(st), 1) * std::pow(v_squ(v), 3) * sf * + std::pow(s, sf - 1.0); + }; + + static Tensor<2, dim, NumberType> + d2psi_ds_x_dt(const SymmetricTensor<2, dim, NumberType> &st, + const Tensor<2, dim, NumberType> & t, + const Tensor<1, dim, NumberType> & v, + const NumberType & s) + { + return 2.0 * std::pow(det_t(st), 2) * std::pow(det_t(t), 1) * ddet_t_dt(t) * + std::pow(v_squ(v), 3) * sf * std::pow(s, sf - 1.0); + }; + + static Tensor<1, dim, NumberType> + d2psi_ds_x_dv(const SymmetricTensor<2, dim, NumberType> &st, + const Tensor<2, dim, NumberType> & t, + const Tensor<1, dim, NumberType> & v, + const NumberType & s) + { + return std::pow(det_t(st), 2) * std::pow(det_t(t), 2) * 3.0 * + std::pow(v_squ(v), 2) * dv_squ_dv(v) * sf * std::pow(s, sf - 1.0); + }; + + static NumberType + d2psi_ds_x_ds(const SymmetricTensor<2, dim, NumberType> &st, + const Tensor<2, dim, NumberType> & t, + const Tensor<1, dim, NumberType> & v, + const NumberType & s) + { + return std::pow(det_t(st), 2) * std::pow(det_t(t), 2) * + std::pow(v_squ(v), 3) * sf * (sf - 1.0) * std::pow(s, sf - 2.0); + }; +}; + +template +const double CoupledFunction::sf(2.2); + +template +void +print(Stream &stream, const std::string &name, const NumberType &val) +{ + stream << name << ": " << val << std::endl; +} +template +void +print(Stream & stream, + const std::string & name, + const Tensor<0, dim, NumberType> &val) +{ + stream << name << ": " << val << std::endl; +} +template +void +print(Stream & stream, + const std::string & name, + const Tensor<1, dim, NumberType> &t) +{ + for (unsigned int i = 0; i < dim; ++i) + stream << name << "[" << i << "]: " << t[i] << std::endl; +} +template +void +print(Stream & stream, + const std::string & name, + const Tensor<2, dim, NumberType> &t) +{ + for (unsigned int i = 0; i < dim; ++i) + for (unsigned int j = 0; j < dim; ++j) + stream << name << "[" << i << "][" << j << "]: " << t[i][j] << std::endl; +} +template +void +print(Stream & stream, + const std::string & name, + const Tensor<3, dim, NumberType> &t) +{ + for (unsigned int i = 0; i < dim; ++i) + for (unsigned int j = 0; j < dim; ++j) + for (unsigned int k = 0; k < dim; ++k) + stream << name << "[" << i << "][" << j << "][" << k + << "]: " << t[i][j][k] << std::endl; +} +template +void +print(Stream & stream, + const std::string & name, + const Tensor<4, dim, NumberType> &t) +{ + for (unsigned int i = 0; i < dim; ++i) + for (unsigned int j = 0; j < dim; ++j) + for (unsigned int k = 0; k < dim; ++k) + for (unsigned int l = 0; l < dim; ++l) + stream << name << "[" << i << "][" << j << "][" << k << "][" << l + << "]: " << t[i][j][k][l] << std::endl; +} + +template +void +test_symmetric_tensor_tensor_vector_scalar_coupled( + const SymmetricTensor<2, dim, number_t> &st, + const Tensor<2, dim, number_t> & t, + const Tensor<1, dim, number_t> & v, + const number_t & s) +{ + using SDNumberType = SD::Expression; + + const SymmetricTensor<2, dim, SDNumberType> symb_st( + SD::make_symmetric_tensor_of_symbols<2, dim>("st")); + const Tensor<2, dim, SDNumberType> symb_t( + SD::make_tensor_of_symbols<2, dim>("t")); + const Tensor<1, dim, SDNumberType> symb_v( + SD::make_vector_of_symbols("v")); + const SDNumberType symb_s(SD::make_symbol("s")); + + // Value + SDNumberType symb_psi; + + // First derivatives + SymmetricTensor<2, dim, SDNumberType> symb_dpsi_dst; + Tensor<2, dim, SDNumberType> symb_dpsi_dt; + Tensor<1, dim, SDNumberType> symb_dpsi_dv; + SDNumberType symb_dpsi_ds; + + // Second derivatives + SymmetricTensor<4, dim, SDNumberType> symb_d2psi_dst_x_dst; + Tensor<4, dim, SDNumberType> symb_d2psi_dst_x_dt; + Tensor<3, dim, SDNumberType> symb_d2psi_dst_x_dv; + SymmetricTensor<2, dim, SDNumberType> symb_d2psi_dst_x_ds; + + Tensor<4, dim, SDNumberType> symb_d2psi_dt_x_dst; + Tensor<4, dim, SDNumberType> symb_d2psi_dt_x_dt; + Tensor<3, dim, SDNumberType> symb_d2psi_dt_x_dv; + Tensor<2, dim, SDNumberType> symb_d2psi_dt_x_ds; + + Tensor<3, dim, SDNumberType> symb_d2psi_dv_x_dst; + Tensor<3, dim, SDNumberType> symb_d2psi_dv_x_dt; + Tensor<2, dim, SDNumberType> symb_d2psi_dv_x_dv; + Tensor<1, dim, SDNumberType> symb_d2psi_dv_x_ds; + + SymmetricTensor<2, dim, SDNumberType> symb_d2psi_ds_x_dst; + Tensor<2, dim, SDNumberType> symb_d2psi_ds_x_dt; + Tensor<1, dim, SDNumberType> symb_d2psi_ds_x_dv; + SDNumberType symb_d2psi_ds_x_ds; + + // Function and its derivatives + using func_sd = CoupledFunction; + symb_psi = func_sd::psi(symb_st, symb_t, symb_v, symb_s); + + // First derivative + symb_dpsi_dst = SD::differentiate(symb_psi, symb_st); + symb_dpsi_dt = SD::differentiate(symb_psi, symb_t); + symb_dpsi_dv = SD::differentiate(symb_psi, symb_v); + symb_dpsi_ds = SD::differentiate(symb_psi, symb_s); + + // Second derivatives + symb_d2psi_dst_x_dst = SD::differentiate(symb_dpsi_dst, symb_st); + symb_d2psi_dst_x_dt = SD::differentiate(symb_dpsi_dst, symb_t); + symb_d2psi_dst_x_dv = SD::differentiate(symb_dpsi_dst, symb_v); + symb_d2psi_dst_x_ds = SD::differentiate(symb_dpsi_dst, symb_s); + + symb_d2psi_dt_x_dst = SD::differentiate(symb_dpsi_dt, symb_st); + symb_d2psi_dt_x_dt = SD::differentiate(symb_dpsi_dt, symb_t); + symb_d2psi_dt_x_dv = SD::differentiate(symb_dpsi_dt, symb_v); + symb_d2psi_dt_x_ds = SD::differentiate(symb_dpsi_dt, symb_s); + + symb_d2psi_dv_x_dst = SD::differentiate(symb_dpsi_dv, symb_st); + symb_d2psi_dv_x_dt = SD::differentiate(symb_dpsi_dv, symb_t); + symb_d2psi_dv_x_dv = SD::differentiate(symb_dpsi_dv, symb_v); + symb_d2psi_dv_x_ds = SD::differentiate(symb_dpsi_dv, symb_s); + + symb_d2psi_ds_x_dst = SD::differentiate(symb_dpsi_ds, symb_st); + symb_d2psi_ds_x_dt = SD::differentiate(symb_dpsi_ds, symb_t); + symb_d2psi_ds_x_dv = SD::differentiate(symb_dpsi_ds, symb_v); + symb_d2psi_ds_x_ds = SD::differentiate(symb_dpsi_ds, symb_s); + +#if DEBUG + print(std::cout, "symb_st", symb_st); + print(std::cout, "symb_t", symb_t); + print(std::cout, "symb_v", symb_v); + print(std::cout, "symb_s", symb_s); + + print(std::cout, "symb_psi", symb_psi); + + print(std::cout, "symb_dpsi_dst", symb_dpsi_dst); + print(std::cout, "symb_dpsi_dt", symb_dpsi_dt); + print(std::cout, "symb_dpsi_dv", symb_dpsi_dv); + print(std::cout, "symb_dpsi_ds", symb_dpsi_ds); + + print(std::cout, "symb_d2psi_dst_x_dst", symb_d2psi_dst_x_dst); + print(std::cout, "symb_d2psi_dst_x_dt", symb_d2psi_dst_x_dt); + print(std::cout, "symb_d2psi_dst_x_dv", symb_d2psi_dst_x_dv); + print(std::cout, "symb_d2psi_dst_x_ds", symb_d2psi_dst_x_ds); + + print(std::cout, "symb_d2psi_dt_x_dst", symb_d2psi_dt_x_dst); + print(std::cout, "symb_d2psi_dt_x_dt", symb_d2psi_dt_x_dt); + print(std::cout, "symb_d2psi_dt_x_dv", symb_d2psi_dt_x_dv); + print(std::cout, "symb_d2psi_dt_x_ds", symb_d2psi_dt_x_ds); + + print(std::cout, "symb_d2psi_dv_x_dst", symb_d2psi_dv_x_dst); + print(std::cout, "symb_d2psi_dv_x_dt", symb_d2psi_dv_x_dt); + print(std::cout, "symb_d2psi_dv_x_dv", symb_d2psi_dv_x_dv); + print(std::cout, "symb_d2psi_dv_x_ds", symb_d2psi_dv_x_ds); + + print(std::cout, "symb_d2psi_ds_x_dst", symb_d2psi_ds_x_dst); + print(std::cout, "symb_d2psi_ds_x_dt", symb_d2psi_ds_x_dt); + print(std::cout, "symb_d2psi_ds_x_dv", symb_d2psi_ds_x_dv); + print(std::cout, "symb_d2psi_ds_x_ds", symb_d2psi_ds_x_ds); +#endif + + // It pains me to do this manually... + // This is why we have utility functions to create substitution maps + // and do evaluation + typename SD::types::substitution_map sub_vals; + sub_vals[symb_s] = SDNumberType(s); + for (unsigned int i = 0; i < dim; ++i) + { + sub_vals[symb_v[i]] = SDNumberType(v[i]); + for (unsigned int j = 0; j < dim; ++j) + { + sub_vals[symb_t[i][j]] = SDNumberType(t[i][j]); + sub_vals[symb_st[i][j]] = SDNumberType(st[i][j]); + } + } + + // Perform substitution and evaluation + // Value + number_t psi; + + // First derivatives + SymmetricTensor<2, dim, number_t> dpsi_dst; + Tensor<2, dim, number_t> dpsi_dt; + Tensor<1, dim, number_t> dpsi_dv; + number_t dpsi_ds; + + // Second derivatives + SymmetricTensor<4, dim, number_t> d2psi_dst_x_dst; + Tensor<4, dim, number_t> d2psi_dst_x_dt; + Tensor<3, dim, number_t> d2psi_dst_x_dv; + SymmetricTensor<2, dim, number_t> d2psi_dst_x_ds; + + Tensor<4, dim, number_t> d2psi_dt_x_dst; + Tensor<4, dim, number_t> d2psi_dt_x_dt; + Tensor<3, dim, number_t> d2psi_dt_x_dv; + Tensor<2, dim, number_t> d2psi_dt_x_ds; + + Tensor<3, dim, number_t> d2psi_dv_x_dst; + Tensor<3, dim, number_t> d2psi_dv_x_dt; + Tensor<2, dim, number_t> d2psi_dv_x_dv; + Tensor<1, dim, number_t> d2psi_dv_x_ds; + + SymmetricTensor<2, dim, number_t> d2psi_ds_x_dst; + Tensor<2, dim, number_t> d2psi_ds_x_dt; + Tensor<1, dim, number_t> d2psi_ds_x_dv; + number_t d2psi_ds_x_ds; + + psi = symb_psi.substitute_and_evaluate(sub_vals); + dpsi_ds = symb_dpsi_ds.substitute_and_evaluate(sub_vals); + d2psi_ds_x_ds = + symb_d2psi_ds_x_ds.substitute_and_evaluate(sub_vals); + for (unsigned int i = 0; i < dim; ++i) + { + dpsi_dv[i] = + symb_dpsi_dv[i].template substitute_and_evaluate(sub_vals); + d2psi_dv_x_ds[i] = + symb_d2psi_dv_x_ds[i].template substitute_and_evaluate( + sub_vals); + d2psi_ds_x_dv[i] = + symb_d2psi_ds_x_dv[i].template substitute_and_evaluate( + sub_vals); + + for (unsigned int j = 0; j < dim; ++j) + { + dpsi_dst[i][j] = + symb_dpsi_dst[i][j].template substitute_and_evaluate( + sub_vals); + dpsi_dt[i][j] = + symb_dpsi_dt[i][j].template substitute_and_evaluate( + sub_vals); + d2psi_dst_x_ds[i][j] = + symb_d2psi_dst_x_ds[i][j] + .template substitute_and_evaluate(sub_vals); + d2psi_dt_x_ds[i][j] = + symb_d2psi_dt_x_ds[i][j].template substitute_and_evaluate( + sub_vals); + d2psi_dv_x_dv[i][j] = + symb_d2psi_dv_x_dv[i][j].template substitute_and_evaluate( + sub_vals); + d2psi_ds_x_dst[i][j] = + symb_d2psi_ds_x_dst[i][j] + .template substitute_and_evaluate(sub_vals); + d2psi_ds_x_dt[i][j] = + symb_d2psi_ds_x_dt[i][j].template substitute_and_evaluate( + sub_vals); + + for (unsigned int k = 0; k < dim; ++k) + { + d2psi_dst_x_dv[i][j][k] = + symb_d2psi_dst_x_dv[i][j][k] + .template substitute_and_evaluate(sub_vals); + d2psi_dt_x_dv[i][j][k] = + symb_d2psi_dt_x_dv[i][j][k] + .template substitute_and_evaluate(sub_vals); + d2psi_dv_x_dst[i][j][k] = + symb_d2psi_dv_x_dst[i][j][k] + .template substitute_and_evaluate(sub_vals); + d2psi_dv_x_dt[i][j][k] = + symb_d2psi_dv_x_dt[i][j][k] + .template substitute_and_evaluate(sub_vals); + + for (unsigned int l = 0; l < dim; ++l) + { + d2psi_dst_x_dst[i][j][k][l] = + symb_d2psi_dst_x_dst[i][j][k][l] + .template substitute_and_evaluate(sub_vals); + d2psi_dst_x_dt[i][j][k][l] = + symb_d2psi_dst_x_dt[i][j][k][l] + .template substitute_and_evaluate(sub_vals); + d2psi_dt_x_dst[i][j][k][l] = + symb_d2psi_dt_x_dst[i][j][k][l] + .template substitute_and_evaluate(sub_vals); + d2psi_dt_x_dt[i][j][k][l] = + symb_d2psi_dt_x_dt[i][j][k][l] + .template substitute_and_evaluate(sub_vals); + } + } + } + } + + // Verify the result + using func = CoupledFunction; + static const double tol = 1e-9; + Assert(std::abs(psi - func::psi(st, t, v, s)) < tol, + ExcMessage("No match for function value.")); + + Assert(std::abs((dpsi_dst - func::dpsi_dst(st, t, v, s)).norm()) < tol, + ExcMessage("No match for first derivative.")); + Assert(std::abs((dpsi_dt - func::dpsi_dt(st, t, v, s)).norm()) < tol, + ExcMessage("No match for first derivative.")); + Assert(std::abs((dpsi_dv - func::dpsi_dv(st, t, v, s)).norm()) < tol, + ExcMessage("No match for first derivative.")); + Assert(std::abs(dpsi_ds - func::dpsi_ds(st, t, v, s)) < tol, + ExcMessage("No match for first derivative.")); + + Assert(std::abs( + (d2psi_dst_x_dst - func::d2psi_dst_x_dst(st, t, v, s)).norm()) < tol, + ExcMessage("No match for second derivative.")); + Assert(std::abs((d2psi_dst_x_dt - func::d2psi_dst_x_dt(st, t, v, s)).norm()) < + tol, + ExcMessage("No match for second derivative.")); + Assert(std::abs((d2psi_dst_x_dv - func::d2psi_dst_x_dv(st, t, v, s)).norm()) < + tol, + ExcMessage("No match for second derivative.")); + Assert(std::abs((d2psi_dst_x_ds - func::d2psi_dst_x_ds(st, t, v, s)).norm()) < + tol, + ExcMessage("No match for second derivative.")); + Assert(std::abs((d2psi_dt_x_dst - func::d2psi_dt_x_dst(st, t, v, s)).norm()) < + tol, + ExcMessage("No match for second derivative.")); + Assert(std::abs((d2psi_dt_x_dt - func::d2psi_dt_x_dt(st, t, v, s)).norm()) < + tol, + ExcMessage("No match for second derivative.")); + Assert(std::abs((d2psi_dt_x_dv - func::d2psi_dt_x_dv(st, t, v, s)).norm()) < + tol, + ExcMessage("No match for second derivative.")); + Assert(std::abs((d2psi_dt_x_ds - func::d2psi_dt_x_ds(st, t, v, s)).norm()) < + tol, + ExcMessage("No match for second derivative.")); + Assert(std::abs((d2psi_dv_x_dst - func::d2psi_dv_x_dst(st, t, v, s)).norm()) < + tol, + ExcMessage("No match for second derivative.")); + Assert(std::abs((d2psi_dv_x_dt - func::d2psi_dv_x_dt(st, t, v, s)).norm()) < + tol, + ExcMessage("No match for second derivative.")); + Assert(std::abs((d2psi_dv_x_dv - func::d2psi_dv_x_dv(st, t, v, s)).norm()) < + tol, + ExcMessage("No match for second derivative.")); + Assert(std::abs((d2psi_dv_x_ds - func::d2psi_dv_x_ds(st, t, v, s)).norm()) < + tol, + ExcMessage("No match for second derivative.")); + Assert(std::abs((d2psi_ds_x_dst - func::d2psi_ds_x_dst(st, t, v, s)).norm()) < + tol, + ExcMessage("No match for second derivative.")); + Assert(std::abs((d2psi_ds_x_dt - func::d2psi_ds_x_dt(st, t, v, s)).norm()) < + tol, + ExcMessage("No match for second derivative.")); + Assert(std::abs((d2psi_ds_x_dv - func::d2psi_ds_x_dv(st, t, v, s)).norm()) < + tol, + ExcMessage("No match for second derivative.")); + Assert(std::abs(d2psi_ds_x_ds - func::d2psi_ds_x_ds(st, t, v, s)) < tol, + ExcMessage("No match for second derivative.")); +} + +template +void +run() +{ + deallog.push("dim " + Utilities::to_string(dim)); + + SymmetricTensor<2, dim, number_t> st = + SymmetricTensor<2, dim, number_t>(unit_symmetric_tensor()); + Tensor<2, dim, number_t> t = + Tensor<2, dim, number_t>(unit_symmetric_tensor()); + Tensor<1, dim, number_t> v; + double s = 0.57; + for (unsigned int i = 0; i < dim; ++i) + { + v[i] = 0.275 * (1.0 + i); + for (unsigned int j = 0; j < dim; ++j) + { + t[i][j] = 1.03 + 1.15 * (((i + 1) * j) % dim); + st[i][j] = 1.01 + 1.25 * (((i + 2) * j) % dim); + } + } + st *= 0.85; + t *= 0.9; + v *= 0.63; + s *= 1.21; + + test_symmetric_tensor_tensor_vector_scalar_coupled(st, + t, + v, + s); + + deallog << "OK" << std::endl; + deallog.pop(); +} + +int +main(int argc, char **argv) +{ + initlog(); + + deallog.push("Double"); + { + run<2, double>(); + run<3, double>(); + } + deallog.pop(); + + // deallog.push("Float"); + // { + // run<2, float>(); + // run<3, float>(); + // } + // deallog.pop(); + + deallog << "OK" << std::endl; +} diff --git a/tests/symengine/symengine_tensor_operations_03.output b/tests/symengine/symengine_tensor_operations_03.output new file mode 100644 index 0000000000..18800bce74 --- /dev/null +++ b/tests/symengine/symengine_tensor_operations_03.output @@ -0,0 +1,4 @@ + +DEAL:Double:dim 2::OK +DEAL:Double:dim 3::OK +DEAL::OK -- 2.39.5