From: Jean-Paul Pelteret Date: Tue, 24 Oct 2017 15:06:34 +0000 (+0200) Subject: Add some modified Sacado::Fad::DFad examples to the testsuite X-Git-Tag: v9.0.0-rc1~883^2 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=dfcce049a8b6b56718f6c1d11e279acc1a142814;p=dealii.git Add some modified Sacado::Fad::DFad examples to the testsuite --- diff --git a/tests/sacado/basic_01.cc b/tests/sacado/basic_01.cc new file mode 100644 index 0000000000..f2155b9c86 --- /dev/null +++ b/tests/sacado/basic_01.cc @@ -0,0 +1,92 @@ +// --------------------------------------------------------------------- +// +// 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 +// forward mode AD. The Sacado::Fad::DFad class, which utilizes dynamic +// memory allocation for the number of derivative components, is selected +// as the number type with which to perform the derivative calculations. +// +// A related example that is shipped with Trilinos can be found at +// https://github.com/trilinos/Trilinos/blob/master/packages/sacado/example/dfad_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; + + // Number of independent variables + const int num_deriv = 2; + + // FAD objects: Independent variables + const Sacado::Fad::DFad x_ad(num_deriv, 0, x); + const Sacado::Fad::DFad y_ad(num_deriv, 1, y); + // FAD objects: Passive variables + const Sacado::Fad::DFad z_ad(z); + + deallog << "x_ad: " << x_ad << std::endl; + deallog << "y_ad: " << y_ad << std::endl; + deallog << "z_ad: " << z_ad << 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::Fad::DFad f_fad = ::f(x_ad, y_ad, z_ad); + + deallog << "f_fad: " << f_fad << std::endl; + + // Extract value and derivatives + const double f_ad = f_fad.val(); // f + const double df_dx_ad = f_fad.dx(0); // df/dx + const double df_dy_ad = f_fad.dx(1); // 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_01.output b/tests/sacado/basic_01.output new file mode 100644 index 0000000000..79d1438967 --- /dev/null +++ b/tests/sacado/basic_01.output @@ -0,0 +1,6 @@ + +DEAL::x_ad: 5.00000 [ 1.00000 0.00000 ] +DEAL::y_ad: 10.0000 [ 0.00000 1.00000 ] +DEAL::z_ad: 4.00000 [ ] +DEAL::f_fad: 380.000 [ 44.0000 36.0000 ] +DEAL::OK diff --git a/tests/sacado/basic_02.cc b/tests/sacado/basic_02.cc new file mode 100644 index 0000000000..7eeab0b687 --- /dev/null +++ b/tests/sacado/basic_02.cc @@ -0,0 +1,121 @@ +// --------------------------------------------------------------------- +// +// 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 forward mode AD. The Sacado::Fad::DFad class, which utilizes +// dynamic memory allocation for the number of derivative components, is +// used for both tiers of the derivative calculations. +// +// A related example that is shipped with Trilinos can be found at +// https://github.com/trilinos/Trilinos/blob/master/packages/sacado/example/dfad_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::Fad::DFad< Sacado::Fad::DFad > x_ad(num_deriv, 0, x); + Sacado::Fad::DFad< Sacado::Fad::DFad > y_ad(num_deriv, 1, y); + // FAD objects: Passive variables + const Sacado::Fad::DFad z_ad(z); + + // Initialize the internal data of variables from which + // second derivatives will be computed + x_ad.val() = Sacado::Fad::DFad(num_deriv, 0, x); + y_ad.val() = Sacado::Fad::DFad(num_deriv, 1, y); + + deallog << "x_ad: " << x_ad << std::endl; + deallog << "y_ad: " << y_ad << std::endl; + deallog << "z_ad: " << z_ad << 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::Fad::DFad< Sacado::Fad::DFad > f_fad = ::f(x_ad, y_ad, z_ad); + + deallog << "f_fad: " << f_fad << std::endl; + + // Extract value and derivatives + const double f_ad = f_fad.val().val(); // f + const double df_dx_ad = f_fad.dx(0).val(); // df/dx + const double df_dy_ad = f_fad.dx(1).val(); // df/dy + const double d2f_dx_dx_ad = f_fad.dx(0).dx(0); // d^2f/dx^2 + const double d2f_dy_dx_ad = f_fad.dx(0).dx(1); // d^2f/dy_dx + const double d2f_dx_dy_ad = f_fad.dx(1).dx(0); // d^2f/dx_dy + const double d2f_dy_dy_ad = f_fad.dx(1).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_02.output b/tests/sacado/basic_02.output new file mode 100644 index 0000000000..1572d95eea --- /dev/null +++ b/tests/sacado/basic_02.output @@ -0,0 +1,6 @@ + +DEAL::x_ad: -3.00000 [ 1.00000 0.00000 ] [ 1.00000 [ ] 0.00000 [ ] ] +DEAL::y_ad: 2.00000 [ 0.00000 1.00000 ] [ 0.00000 [ ] 1.00000 [ ] ] +DEAL::z_ad: 7.00000 [ ] +DEAL::f_fad: -35.0000 [ 203.000 154.000 ] [ 203.000 [ -126.000 14.0000 ] 154.000 [ 14.0000 77.0000 ] ] +DEAL::OK