From bb4f404c97025938052e1378fd2b338a4a20a34c Mon Sep 17 00:00:00 2001 From: Jean-Paul Pelteret Date: Thu, 26 Oct 2017 16:00:45 +0200 Subject: [PATCH] Add some Sacado::RFad and nested RFad-DFad examples to the testsuite --- tests/sacado/basic_01a.cc | 88 +++++++++++++++++++++++++ tests/sacado/basic_01a.output | 6 ++ tests/sacado/basic_02a.cc | 118 ++++++++++++++++++++++++++++++++++ tests/sacado/basic_02a.output | 6 ++ 4 files changed, 218 insertions(+) create mode 100644 tests/sacado/basic_01a.cc create mode 100644 tests/sacado/basic_01a.output create mode 100644 tests/sacado/basic_02a.cc create mode 100644 tests/sacado/basic_02a.output diff --git a/tests/sacado/basic_01a.cc b/tests/sacado/basic_01a.cc new file mode 100644 index 0000000000..0fd88bc73c --- /dev/null +++ b/tests/sacado/basic_01a.cc @@ -0,0 +1,88 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2017 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. +// +// --------------------------------------------------------------------- + +// Tests the computation of the first derivatives of a function using +// reverse mode AD with the Sacado::Rad::ADvar class. +// +// A related example that is shipped with Trilinos can be found at +// https://github.com/trilinos/Trilinos/blob/master/packages/sacado/example/trad_example.cpp + + +#include "../tests.h" + +#include + +// The function to differentiate +template +NumberType f(const NumberType &x, const NumberType &y, const NumberType &z) +{ + return z*(x + z*y + x*y); +} + +// The analytic derivative of f(x,y,z) with respect to x and y +void +df(const double &x, const double &y, const double &z, + double &df_dx, double &df_dy) +{ + df_dx = z*(1.0 + y); + df_dy = z*(z + x); +} + +int main() +{ + initlog(); + + // Values of function arguments + const double x = 5.0; + const double y = 10.0; + const double z = 4.0; + + // RAD objects: Independent variables + const Sacado::Rad::ADvar x_ad (x); + const Sacado::Rad::ADvar y_ad (y); + // RAD objects: Passive variables + const Sacado::Rad::ADvar z_ad (z); + + deallog << "x_ad: " << x_ad.val() << std::endl; + deallog << "y_ad: " << y_ad.val() << std::endl; + deallog << "z_ad: " << z_ad.val() << std::endl; + + // Compute function + const double f = ::f(x, y, z); + + // Compute derivative analytically + double df_dx = 0.0, df_dy = 0.0; + df(x, y, z, df_dx, df_dy); + + // Compute function and derivative with AD + const Sacado::Rad::ADvar f_rad = ::f(x_ad, y_ad, z_ad); + Sacado::Rad::ADvar::Gradcomp(); + + deallog << "f_rad: " << f_rad.val() << std::endl; + + // Extract value and derivatives + const double f_ad = f_rad.val(); // f + const double df_dx_ad = x_ad.adj(); // df/dx + const double df_dy_ad = y_ad.adj(); // df/dy + + const double tol = 1.0e-14; + Assert(std::fabs(f - f_ad) < tol, + ExcMessage("Computation incorrect: Value")); + Assert(std::fabs(df_dx - df_dx_ad) < tol && + std::fabs(df_dy - df_dy_ad) < tol, + ExcMessage("Computation incorrect: First derivative")); + + deallog << "OK" << std::endl; +} diff --git a/tests/sacado/basic_01a.output b/tests/sacado/basic_01a.output new file mode 100644 index 0000000000..99ded411f2 --- /dev/null +++ b/tests/sacado/basic_01a.output @@ -0,0 +1,6 @@ + +DEAL::x_ad: 5.00000 +DEAL::y_ad: 10.0000 +DEAL::z_ad: 4.00000 +DEAL::f_rad: 380.000 +DEAL::OK diff --git a/tests/sacado/basic_02a.cc b/tests/sacado/basic_02a.cc new file mode 100644 index 0000000000..c3144644b7 --- /dev/null +++ b/tests/sacado/basic_02a.cc @@ -0,0 +1,118 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2017 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. +// +// --------------------------------------------------------------------- + +// Tests the computation of the second derivatives of a function using +// nested reverse-forward mode AD. The Sacado::Rad::ADvar class is used to +// compute the first derivatives while the Sacado::Fad::DFad class, which utilizes +// dynamic memory allocation for the number of derivative components, is +// used for the second derivative calculations. +// +// A related example that is shipped with Trilinos can be found at +// https://github.com/trilinos/Trilinos/blob/master/packages/sacado/example/trad_dfad_example.cpp + + +#include "../tests.h" + +#include + +// The function to differentiate +template +NumberType +f(const NumberType &x, const NumberType &y, const NumberType2 &z) +{ + return z*(x*x*x + z*y*y + 0.5*x*y*y); +} + +// The analytic derivative of f(x,y,z) with respect to x and y +void +df(const double &x, const double &y, const double &z, + double &df_dx, double &df_dy) +{ + df_dx = z*(3.0*x*x + 0.5*y*y); + df_dy = z*(2.0*z*y + x*y); +} + +// The analytic second derivatives of f(x,y,z) with respect to x and y +void +d2f(const double &x, const double &y, const double &z, + double &d2f_dx_dx, double &d2f_dy_dy, + double &d2f_dy_dx) +{ + d2f_dx_dx = z*(6.0*x); + d2f_dy_dx = z*y; + d2f_dy_dy = z*(2.0*z + x); +} + +int main() +{ + initlog(); + + // Values of function arguments + const double x = -3.0; + const double y = 2.0; + const double z = 7.0; + + // Number of independent variables + const int num_deriv = 2; + + // FAD objects: Independent variables + Sacado::Rad::ADvar< Sacado::Fad::DFad > x_ad(Sacado::Fad::DFad(num_deriv, 0, x)); + Sacado::Rad::ADvar< Sacado::Fad::DFad > y_ad(Sacado::Fad::DFad(num_deriv, 1, y)); + // FAD objects: Passive variables + const Sacado::Rad::ADvar< Sacado::Fad::DFad > z_ad(z); + + deallog << "x_ad: " << x_ad.val() << std::endl; + deallog << "y_ad: " << y_ad.val() << std::endl; + deallog << "z_ad: " << z_ad.val() << std::endl; + + // Compute function + const double f = ::f(x, y, z); + + // Compute derivative analytically + double df_dx = 0.0, df_dy = 0.0; + df(x, y, z, df_dx, df_dy); + + // Compute second derivative analytically + double d2f_dx_dx = 0.0, d2f_dy_dy = 0.0, d2f_dy_dx = 0.0; + d2f(x, y, z, d2f_dx_dx, d2f_dy_dy, d2f_dy_dx); + + // Compute function and derivative with AD + const Sacado::Rad::ADvar< Sacado::Fad::DFad > f_rfad = ::f(x_ad, y_ad, z_ad); + Sacado::Rad::ADvar< Sacado::Fad::DFad >::Gradcomp(); + + deallog << "f_rad: " << f_rfad.val() << std::endl; + + // Extract value and derivatives + const double f_ad = f_rfad.val().val(); // f + const double df_dx_ad = x_ad.adj().val(); // df/dx + const double df_dy_ad = y_ad.adj().val(); // df/dy + const double d2f_dx_dx_ad = x_ad.adj().dx(0); // d^2f/dx^2 + const double d2f_dy_dx_ad = x_ad.adj().dx(1); // d^2f/dy_dx + const double d2f_dx_dy_ad = y_ad.adj().dx(0); // d^2f/dx_dy + const double d2f_dy_dy_ad = y_ad.adj().dx(1); // d^2f/dy^2 + + const double tol = 1.0e-14; + Assert(std::fabs(f - f_ad) < tol, + ExcMessage("Computation incorrect: Value")); + Assert(std::fabs(df_dx - df_dx_ad) < tol && + std::fabs(df_dy - df_dy_ad) < tol, + ExcMessage("Computation incorrect: First derivative")); + Assert(std::fabs(d2f_dx_dx - d2f_dx_dx_ad) < tol && + std::fabs(d2f_dy_dy - d2f_dy_dy_ad) < tol && + std::fabs(d2f_dy_dx - d2f_dy_dx_ad) < tol, + ExcMessage("Computation incorrect: Second derivative")); + + deallog << "OK" << std::endl; +} diff --git a/tests/sacado/basic_02a.output b/tests/sacado/basic_02a.output new file mode 100644 index 0000000000..f748d6e824 --- /dev/null +++ b/tests/sacado/basic_02a.output @@ -0,0 +1,6 @@ + +DEAL::x_ad: -3.00000 [ 1.00000 0.00000 ] +DEAL::y_ad: 2.00000 [ 0.00000 1.00000 ] +DEAL::z_ad: 7.00000 [ ] +DEAL::f_rad: -35.0000 [ 203.000 154.000 ] +DEAL::OK -- 2.39.5