From 5c06092fa54191d2c19f2cae4615fa63a12b460a Mon Sep 17 00:00:00 2001 From: Ryan Lutz Date: Tue, 17 Jun 2025 10:36:31 -0700 Subject: [PATCH 1/4] added example and test --- src/examples/CMakeLists.txt | 1 + src/examples/neohookean_timing.cpp | 727 +++++++++++++++++++++++++++++ src/tests/CMakeLists.txt | 1 + src/tests/enzyme_neohookean.cpp | 624 +++++++++++++++++++++++++ 4 files changed, 1353 insertions(+) create mode 100644 src/examples/neohookean_timing.cpp create mode 100644 src/tests/enzyme_neohookean.cpp diff --git a/src/examples/CMakeLists.txt b/src/examples/CMakeLists.txt index d0345d3a..0c7a30d1 100644 --- a/src/examples/CMakeLists.txt +++ b/src/examples/CMakeLists.txt @@ -8,6 +8,7 @@ set( contact_examples common_plane_gpu.cpp common_plane.cpp mortar_lm_patch_test.cpp + neohookean_timing.cpp ) diff --git a/src/examples/neohookean_timing.cpp b/src/examples/neohookean_timing.cpp new file mode 100644 index 00000000..a694de55 --- /dev/null +++ b/src/examples/neohookean_timing.cpp @@ -0,0 +1,727 @@ +#include +#include +#include +#include "gtest/gtest.h" +#include +#include "axom/core.hpp" +#include "axom/slic/interface/slic.hpp" +#include +#include +#include "tribol/common/Enzyme.hpp" + + +void multiply3x3(const double F[9], const double F_T[9], double C[9]) { + for (int i = 0; i < 3; ++i) { // row of A + for (int j = 0; j < 3; ++j) { // column of B + C[i * 3 + j] = 0.0; + for (int k = 0; k < 3; ++k) { + C[i * 3 + j] += F_T[i * 3 + k] * F[k * 3 + j]; + } + } + } +} + +void calc_E_from_F(const double F[9], double E[9]) { + double I[9] = {1.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 1.0}; + double F_T[9] = {0.0}; + F_T[0] = F[0]; + F_T[1] = F[3]; + F_T[2] = F[6]; + F_T[3] = F[1]; + F_T[4] = F[4]; + F_T[5] = F[7]; + F_T[6] = F[2]; + F_T[7] = F[5]; + F_T[8] = F[8]; + double C[9] = {0.0}; + multiply3x3(F, F_T, C); + for(int i = 0; i < 9; ++i) { + E[i] = 0.5 * (C[i] - I[i]); + } +} + +//Calculates right cauchy stress tensor +void calc_cauchy_stress_tensor(const double E[9], double C[9]){ + double I[9] = {1, 0, 0, 0, 1, 0, 0, 0, 1}; + for(int i = 0; i < 9; ++i) { + C[i] = 2 * E[i]; + C[i] += I[i]; + } +} + + +//Calculates Trace +double calc_trace(const double C[9]) { + double Tr_C = C[0] + C[4] + C[8]; + return Tr_C; +} + + + +double calc_determinant(const double C[9]) { + double J = C[0] * (C[4] * C[8] - C[5] * C[7]) - C[1] * (C[3] * C[8] - C[5] * C[6]) + C[2] * (C[3] * C[7] - C[4] * C[6]); + J = std::sqrt(J); + return J; +} + + +bool invert3x3(const double F[9], double Finv[9]) +{ + // Compute the determinant + double det = + F[0]*(F[4]*F[8] - F[5]*F[7]) - + F[1]*(F[3]*F[8] - F[5]*F[6]) + + F[2]*(F[3]*F[7] - F[4]*F[6]); + + if (std::abs(det) < std::numeric_limits::epsilon()) + return false; // Singular matrix + + double invDet = 1.0 / det; + + // Compute the inverse using the formula for the inverse of a 3x3 matrix + Finv[0] = (F[4]*F[8] - F[5]*F[7]) * invDet; + Finv[1] = -(F[1]*F[8] - F[2]*F[7]) * invDet; + Finv[2] = (F[1]*F[5] - F[2]*F[4]) * invDet; + + Finv[3] = -(F[3]*F[8] - F[5]*F[6]) * invDet; + Finv[4] = (F[0]*F[8] - F[2]*F[6]) * invDet; + Finv[5] = -(F[0]*F[5] - F[2]*F[3]) * invDet; + + Finv[6] = (F[3]*F[7] - F[4]*F[6]) * invDet; + Finv[7] = -(F[0]*F[7] - F[1]*F[6]) * invDet; + Finv[8] = (F[0]*F[4] - F[1]*F[3]) * invDet; + + return true; +} + + +//build strain energy equation +void strain_energy(double* E, double mu, double lambda, double* W) { + double C[9] = {0.0}; + calc_cauchy_stress_tensor(E, C); + double Tr_C = calc_trace(C); + double J = calc_determinant(C); + *W = mu/2.0 * (Tr_C - 3.0) - mu * log(J) + lambda/2.0 * pow((log(J)), 2.0); +} + +//calc stress using enzyme fwddiff +void stress(double* E, double mu, double lambda, double* dW_dE) { + double W = 0.0; + for(int i = 0; i < 9; ++i) { + double dE[9] = {0.0}; + dE[i] = 1.0; + double dmu = 0.0; + double dlambda = 0.0; + double dw = 0.0; + __enzyme_fwddiff( (void*) strain_energy, E, dE, mu, dmu, lambda, dlambda, &W, &dw); + dW_dE[i] = dw; + } + +} + +//calc stress using enzyme autodiff +void stress_reverse(double* E, double mu, double lambda, double* dW_dE) { + double dE[9] = {0.0}; + double W = 0.0; + double dW = 1.0; + __enzyme_autodiff( strain_energy, enzyme_dup, E, dE, enzyme_const, mu, enzyme_const, lambda, enzyme_dup, &W, &dW); + + for (int i = 0; i < 9; ++i) { + dW_dE[i] = dE[i]; + } +} + + +//calc stress using finite Difference +void stress_FD(double* E, double mu, double lambda, double* dW_dE, double h = 1e-7) { + double E_plus[9] = {0.0}; + double E_minus[9] = {0.0}; + double W_plus; + double W_minus; + for(int i = 0; i < 9; ++i) { + for(int j = 0; j < 9; ++j) { + E_plus[j] = E[j]; + E_minus[j] = E[j]; + } + E_plus[i] = E[i] + h; + E_minus[i] = E[i] - h; + strain_energy(E_plus, mu, lambda, &W_plus); + strain_energy(E_minus, mu, lambda, &W_minus); + dW_dE[i] = (W_plus - W_minus) / (2 * h); + + } + +} + + +void hand_code_deriv(double* E, double mu, double lambda, double* S) { + double C[9] = {0.0}; + double Cinv[9]; + double I[9] = {1.0, 0.0, 0.0, 0, 1.0, 0.0, 0.0, 0.0, 1.0}; + calc_cauchy_stress_tensor(E, C); + double J = calc_determinant(C); + invert3x3(C, Cinv); + double first_term[9]; + for(int i = 0; i < 9; ++i) { + first_term[i] = lambda * std::log(J) * Cinv[i]; + } + double second_term[9]; + for(int i = 0; i < 9; ++i) { + second_term[i] = I[i] - Cinv[i]; + second_term[i] *= mu; + } + for(int i = 0; i < 9; ++i) { + S[i] = first_term[i] + second_term[i]; + } + +} + +void second_deriv_fwd_fwd(double* E, double mu, double lambda, double* d2W_d2E) { + double dW[9] = {0.0}; + double d2w[9] = {0.0}; + + for(int i = 0; i < 9; ++i) { + double d2E[9] = {0.0}; + d2E[i] = 1.0; + double d2mu = 0.0; + double d2lambda = 0.0; + __enzyme_fwddiff ( (void*) stress, E, d2E, mu, d2mu, lambda, d2lambda, &dW, &d2w ); + for(int j = 0; j < 9; ++j) { + d2W_d2E[9 * i + j] = d2w[j]; + } + } +} + +void second_deriv_rev_fwd(double* E, double mu, double lambda, double* d2W_d2E) { + double dW[9] = {0.0}; + double d2w[9] = {0.0}; + + for(int i = 0; i < 9; ++i) { + double d2E[9] = {0.0}; + d2E[i] = 1.0; + double d2mu = 0.0; + double d2lambda = 0.0; + __enzyme_fwddiff ( (void*) stress_reverse, E, d2E, mu, d2mu, lambda, d2lambda, &dW, &d2w ); + for(int j = 0; j < 9; ++j) { + d2W_d2E[9 * i + j] = d2w[j]; + } + } +} + +void second_deriv_rev_rev(double* E, double mu, double lambda, double* d2W_d2E) { + for (int i = 0; i < 9; ++i) { + double d2E[81] = {0.0}; + double W[9] = {0.0}; + double d2W[9] = {0.0}; + d2W[i] = 1.0; + __enzyme_autodiff( stress_reverse, enzyme_dup, E, d2E, enzyme_const, mu, enzyme_const, lambda, enzyme_dup, &W, &d2W); + for(int j = 0; j < 9; ++j) { + d2W_d2E[9 * j + i] = d2E[j]; + } + } +} + +void second_deriv_fwd_rev(double* E, double mu, double lambda, double* d2W_d2E) { + for (int i = 0; i < 9; ++i) { + double d2E[81] = {0.0}; + double W[9] = {0.0}; + double d2W[9] = {0.0}; + d2W[i] = 1.0; + __enzyme_autodiff( stress, enzyme_dup, E, d2E, enzyme_const, mu, enzyme_const, lambda, enzyme_dup, &W, &d2W); + for(int j = 0; j < 9; ++j) { + d2W_d2E[9 * i + j] = d2E[j]; + } + } +} + + +void second_deriv_fwd_FD(double* E, double mu, double lambda, double* d2W_d2E, double h = 1e-7){ + double E_plus[9] = {0.0}; + double E_minus[9] = {0.0}; + double dW_plus[9] = {0.0}; + double dW_minus[9] = {0.0}; + for(int i = 0; i < 9; ++i) { + for(int j = 0; j < 9; ++j) { + E_plus[j] = E[j]; + E_minus[j] = E[j]; + } + E_plus[i] = E[i] + h; + E_minus[i] = E[i] - h; + stress(E_plus, mu, lambda, dW_plus); + stress(E_minus, mu, lambda, dW_minus); + for(int j = 0; j < 9; ++j) { + d2W_d2E[9 * i + j] = (dW_plus[j] - dW_minus[j]) / (2 * h); + } + } +} + +void second_deriv_hand_fwd(double *E, double mu, double lambda, double* d2W_d2E) { + double dW[9] = {0.0}; + double d2w[9] = {0.0}; + + + for(int i = 0; i < 9; ++i) { + double d2E[9] = {0.0}; + d2E[i] = 1.0; + double d2mu = 0.0; + double d2lambda = 0.0; + __enzyme_fwddiff ( (void*) hand_code_deriv, E, d2E, mu, d2mu, lambda, d2lambda, &dW, &d2w ); + for(int j = 0; j < 9; ++j) { + d2W_d2E[9 * i + j] = d2w[j]; + } + } +} + +void second_deriv_hand_FD(double* E, double mu, double lambda, double* d2W_d2E, double h = 1e-7) { + double E_plus[9] = {0.0}; + double E_minus[9] = {0.0}; + double dw_minus[9] = {0.0}; + double dw_plus[9] = {0.0}; + for(int i = 0; i < 9; ++i){ + for (int j = 0; j < 9; ++j) { + E_plus[j] = E[j]; + E_minus[j] = E[j]; + } + E_plus[i] = E[i] + h; + E_minus[i] = E[i] - h; + hand_code_deriv(E_plus, mu, lambda, dw_plus); + hand_code_deriv(E_minus, mu, lambda, dw_minus); + for(int j = 0; j < 9; ++j) { + d2W_d2E[9 * i + j] = (dw_plus[j] - dw_minus[j]) / (2 * h); + } + } +} + + + + + +void run_fwd_mode(double* E, double mu, double lambda, double* dw_df, int N) { + axom::utilities::Timer timer{ false }; + stress(E, mu, lambda, dw_df); + double Dw[9] = {0.0}; + timer.start(); + for(int i = 0; i < N; ++i) { + stress(E, mu, lambda, dw_df); + for(int j = 0; j < 9; ++j) { + Dw[j] += dw_df[j]; + } + } + timer.stop(); + std::cout << axom::fmt::format( "Time to calc fwd_diff: {0:f}ms", timer.elapsedTimeInMilliSec() ) << std::endl; +} + +void run_bkwd_mode(double* E, double mu, double lambda, double* dw_dE, int N) { + axom::utilities::Timer timer{ false }; + + stress_reverse(E, mu, lambda, dw_dE); + double Dw[9] = {0.0}; + + timer.start(); + for(int i = 0; i < N; ++i) { + stress_reverse(E, mu, lambda, dw_dE); + for(int j = 0; j < 9; ++j) { + Dw[j] += dw_dE[j]; + } + } + timer.stop(); + std::cout << axom::fmt::format( "Time to calc backward_diff: {0:f}ms", timer.elapsedTimeInMilliSec() ) << std::endl; +} + +void run_hand_derivative(double* E, double mu, double lambda, double* S, int N) { + axom::utilities::Timer timer{ false }; + hand_code_deriv(E, mu, lambda, S); + double Dw[9] = {0.0}; + timer.start(); + for(int i = 0; i < N; ++i) { + hand_code_deriv(E, mu, lambda, S); + for(int j = 0; j < 9; ++j) { + Dw[j] += S[j]; + } + } + timer.stop(); + std::cout << axom::fmt::format( "Time to calc hand_diff: {0:f}ms", timer.elapsedTimeInMilliSec() ) << std::endl; +} + +void run_fwd_fwd(double* E, double mu, double lambda, double* S, int N) { + axom::utilities::Timer timer{ false }; + timer.start(); + for(int i = 0; i < N; ++i) { + second_deriv_fwd_fwd(E, mu, lambda, S); + } + timer.stop(); + std::cout << axom::fmt::format( "Time to calc fwd_fwd_diff: {0:f}ms", timer.elapsedTimeInMilliSec() ) << std::endl; +} + +void run_rev_fwd(double* E, double mu, double lambda, double* S, int N) { + axom::utilities::Timer timer{ false }; + timer.start(); + for(int i = 0; i < N; ++i) { + second_deriv_rev_fwd(E, mu, lambda, S); + } + timer.stop(); + std::cout << axom::fmt::format( "Time to calc rev_fwd_diff: {0:f}ms", timer.elapsedTimeInMilliSec() ) << std::endl; +} + +void run_fwd_rev(double* E, double mu, double lambda, double* S, int N) { + axom::utilities::Timer timer{ false }; + timer.start(); + for(int i = 0; i < N; ++i) { + second_deriv_fwd_rev(E, mu, lambda, S); + } + timer.stop(); + std::cout << axom::fmt::format( "Time to calc fwd_rev_diff: {0:f}ms", timer.elapsedTimeInMilliSec() ) << std::endl; +} + +void run_rev_rev(double* E, double mu, double lambda, double* S, int N) { + axom::utilities::Timer timer{ false }; + timer.start(); + for(int i = 0; i < N; ++i) { + second_deriv_rev_rev(E, mu, lambda, S); + } + timer.stop(); + std::cout << axom::fmt::format( "Time to calc rev_rev_diff: {0:f}ms", timer.elapsedTimeInMilliSec() ) << std::endl; +} + +void run_fwd_FD(double* E, double mu, double lambda, double* S, int N, double h) { + axom::utilities::Timer timer{ false }; + timer.start(); + for(int i = 0; i < N; ++i) { + second_deriv_fwd_FD(E, mu, lambda, S, h); + } + timer.stop(); + std::cout << axom::fmt::format( "Time to calc fwd_FD_diff: {0:f}ms", timer.elapsedTimeInMilliSec() ) << std::endl; +} + +void run_hand_fwd(double* E, double mu, double lambda, double* S, int N) { + axom::utilities::Timer timer{ false }; + timer.start(); + for(int i = 0; i < N; ++i) { + second_deriv_hand_fwd(E, mu, lambda, S); + } + timer.stop(); + std::cout << axom::fmt::format( "Time to calc hand_fwd_diff: {0:f}ms", timer.elapsedTimeInMilliSec() ) << std::endl; +} + + +int main() { + double E[9] = {0.1, 0.05, 0.02, 0.05, 0.2, 0.01, 0.02, 0.01, 0.15}; + double mu = 1.0; + double lambda = 1.0; + int N = 0; + std::cout << "Enter Number: "; + std::cin >> N; + double dw_dE[9] = {0.0}; + double d2W_d2E[81] = {0.0}; + double h = 1e-7; + + run_fwd_mode(E, mu, lambda, dw_dE, N); + run_bkwd_mode(E, mu, lambda, dw_dE, N); + run_hand_derivative(E, mu, lambda, dw_dE, N); + run_fwd_fwd(E, mu, lambda, d2W_d2E, N); + run_rev_fwd(E, mu, lambda, d2W_d2E, N); + run_fwd_rev(E, mu, lambda, d2W_d2E, N); + run_rev_rev(E, mu, lambda, d2W_d2E, N); + run_fwd_FD(E, mu, lambda, d2W_d2E, N, h); + run_hand_fwd(E, mu, lambda, d2W_d2E, N); + + + return 0; + +} + +// TEST(StrainEnergyTest, StressFiniteDifferenceVsAutodiff) { +// double F[9] = {1.0, 0.5, 0.0, 0.0, 1.2, 0.1, 0.0, 0.0, 1.0}; +// double E[9] = {0.0}; +// double mu = 1.0; +// double lambda = 1.0; +// double dW_dF_fd[9]; +// double dW_dF_ad[9]; +// double J = calc_determinant(F); +// calc_green_lagrange(F, E); +// std::cout << "Calling stress_FD..." << std::endl; +// stress_FD(E, J, mu, lambda, dW_dF_fd); +// std::cout << "Finite difference stress computed:" << std::endl; +// for (int i = 0; i < 9; ++i) { +// std::cout << "dW_dF_fd[" << i << "] = " << dW_dF_fd[i] << std::endl; +// } + +// std::cout << "Calling stress..." << std::endl; +// stress(E, J, mu, lambda, dW_dF_ad); +// std::cout << "Autodiff stress computed:" << std::endl; +// for (int i = 0; i < 9; ++i) { +// std::cout << "dW_dF_ad[" << i << "] = " << dW_dF_ad[i] << std::endl; +// } + +// for (int i = 0; i < 9; ++i) { +// std::cout << "Comparing index " << i << ": FD = " << dW_dF_fd[i] +// << ", AD = " << dW_dF_ad[i] << std::endl; +// EXPECT_NEAR(dW_dF_fd[i], dW_dF_ad[i], 1e-6) << "Mismatch at index " << i; +// } +// } + +// TEST(StrainEnergyTest, StressFiniteDifferenceVsReverseAutodiff) { +// double E[9] = {1.0, 0.5, 0.0, 0.0, 1.2, 0.1, 0.0, 0.0, 1.0}; +// double mu = 1.0; +// double lambda = 1.0; +// double dW_dF_fd[9]; +// double dW_dF_rev[9]; +// std::cout << "E: "; +// for (int i = 0; i < 9; ++i) std::cout << E[i] << " "; +// std::cout << std::endl; +// std::cout << "Calling stress_FD..." << std::endl; +// stress_FD(E, mu, lambda, dW_dF_fd); +// std::cout << "Finite difference stress computed:" << std::endl; +// for (int i = 0; i < 9; ++i) { +// std::cout << "dW_dF_fd[" << i << "] = " << dW_dF_fd[i] << std::endl; +// } + +// std::cout << "Calling stress_reverse..." << std::endl; +// stress_reverse(E, mu, lambda, dW_dF_rev); +// std::cout << "Reverse-mode autodiff stress computed:" << std::endl; +// for (int i = 0; i < 9; ++i) { +// std::cout << "dW_dF_rev[" << i << "] = " << dW_dF_rev[i] << std::endl; +// } + +// for (int i = 0; i < 9; ++i) { +// std::cout << "Comparing index " << i << ": FD = " << dW_dF_fd[i] +// << ", Reverse AD = " << dW_dF_rev[i] << std::endl; +// EXPECT_NEAR(dW_dF_fd[i], dW_dF_rev[i], 1e-6) << "Mismatch at index " << i; +// } +// } + +// TEST(StrainEnergyTest, StressFiniteDifferenceVsReverseAutodiff) { +// double E[9] = {0.1, 0.05, 0.02, 0.05, 0.2, 0.01, 0.02, 0.01, 0.15}; +// double mu = 1.0; +// double lambda = 1.0; +// double dW_dF_fd[9]; +// double dW_dF_hand[9]; +// std::cout << "E: "; +// for (int i = 0; i < 9; ++i) std::cout << E[i] << " "; +// std::cout << std::endl; +// std::cout << "Calling stress_FD..." << std::endl; +// stress_FD(E, mu, lambda, dW_dF_fd); +// std::cout << "Finite difference stress computed:" << std::endl; +// for (int i = 0; i < 9; ++i) { +// std::cout << "dW_dF_fd[" << i << "] = " << dW_dF_fd[i] << std::endl; +// } + +// std::cout << "Calling hand_code_deriv..." << std::endl; +// hand_code_deriv(E, mu, lambda, dW_dF_hand); +// std::cout << "Hand coded stress computed:" << std::endl; +// for (int i = 0; i < 9; ++i) { +// std::cout << "dW_dF_hand[" << i << "] = " << dW_dF_hand[i] << std::endl; +// } + +// for (int i = 0; i < 9; ++i) { +// std::cout << "Comparing index " << i << ": FD = " << dW_dF_fd[i] +// << ", Hand derivative = " << dW_dF_hand[i] << std::endl; +// EXPECT_NEAR(dW_dF_fd[i], dW_dF_hand[i], 1e-6) << "Mismatch at index " << i; +// } +// } + +// TEST(StrainEnergyTest, StressFiniteDifferenceVsReverseAutodiff) { +// double E[9] = {0.1, 0.05, 0.02, 0.05, 0.2, 0.01, 0.02, 0.01, 0.15}; +// double mu = 1.0; +// double lambda = 1.0; +// double dW_dF_fwd_fd[81]; +// double dW_dF_fwd_fwd[81]; +// std::cout << "E: "; +// for (int i = 0; i < 81; ++i) std::cout << E[i] << " "; +// std::cout << std::endl; +// std::cout << "Calling stress_FD..." << std::endl; +// second_deriv_fwd_FD(E, mu, lambda, dW_dF_fwd_fd); +// std::cout << "Finite difference stress computed:" << std::endl; +// for (int i = 0; i < 81; ++i) { +// std::cout << "dW_dF_fd[" << i << "] = " << dW_dF_fwd_fd[i] << std::endl; +// } + +// std::cout << "Calling hand_code_deriv..." << std::endl; +// second_deriv_fwd_fwd(E, mu, lambda, dW_dF_fwd_fwd); +// std::cout << "Hand coded stress computed:" << std::endl; +// for (int i = 0; i < 81; ++i) { +// std::cout << "dW_dF_hand[" << i << "] = " << dW_dF_fwd_fwd[i] << std::endl; +// } + +// for (int i = 0; i < 81; ++i) { +// std::cout << "Comparing index " << i << ": FD = " << dW_dF_fwd_fd[i] +// << ", Hand derivative = " << dW_dF_fwd_fwd[i] << std::endl; +// EXPECT_NEAR(dW_dF_fwd_fd[i], dW_dF_fwd_fwd[i], 1e-6) << "Mismatch at index " << i; +// } +// } + +// TEST(StrainEnergyTest, StressFiniteDifferenceVsReverseAutodiff) { +// double E[9] = {0.1, 0.05, 0.02, 0.05, 0.2, 0.01, 0.02, 0.01, 0.15}; +// double mu = 1.0; +// double lambda = 1.0; +// double dW_dF_fwd_fd[81]; +// double dW_dF_fwd_rev[81]; + +// std::cout << "E: "; +// for (int i = 0; i < 9; ++i) std::cout << E[i] << " "; +// std::cout << std::endl; + +// std::cout << "Calling stress_FD..." << std::endl; +// second_deriv_fwd_FD(E, mu, lambda, dW_dF_fwd_fd); +// std::cout << "Finite difference stress computed:" << std::endl; +// for (int i = 0; i < 81; ++i) { +// std::cout << "dW_dF_fd[" << i << "] = " << dW_dF_fwd_fd[i] << std::endl; +// } + +// std::cout << "Calling hand_code_deriv (fwd_rev)..." << std::endl; +// second_deriv_fwd_rev(E, mu, lambda, dW_dF_fwd_rev); +// std::cout << "Hand coded stress computed (fwd_rev):" << std::endl; +// for (int i = 0; i < 81; ++i) { +// std::cout << "dW_dF_hand[" << i << "] = " << dW_dF_fwd_rev[i] << std::endl; +// } + +// for (int i = 0; i < 81; ++i) { +// std::cout << "Comparing index " << i << ": FD = " << dW_dF_fwd_fd[i] +// << ", Hand derivative = " << dW_dF_fwd_rev[i] << std::endl; +// EXPECT_NEAR(dW_dF_fwd_fd[i], dW_dF_fwd_rev[i], 1e-6) << "Mismatch at index " << i; +// } +// } + +// TEST(StrainEnergyTest, StressFiniteDifferenceVsReverseAutodiff) { +// double E[9] = {0.1, 0.05, 0.02, 0.05, 0.2, 0.01, 0.02, 0.01, 0.15}; +// double mu = 1.0; +// double lambda = 1.0; +// double dW_dF_fwd_fd[81]; +// double dW_dF_rev_rev[81]; + +// std::cout << "E: "; +// for (int i = 0; i < 9; ++i) std::cout << E[i] << " "; +// std::cout << std::endl; + +// std::cout << "Calling stress_FD..." << std::endl; +// second_deriv_fwd_FD(E, mu, lambda, dW_dF_fwd_fd); +// std::cout << "Finite difference stress computed:" << std::endl; +// for (int i = 0; i < 81; ++i) { +// std::cout << "dW_dF_fd[" << i << "] = " << dW_dF_fwd_fd[i] << std::endl; +// } + +// std::cout << "Calling hand_code_deriv (rev_rev)..." << std::endl; +// second_deriv_rev_rev(E, mu, lambda, dW_dF_rev_rev); +// std::cout << "Hand coded stress computed (rev_rev):" << std::endl; +// for (int i = 0; i < 81; ++i) { +// std::cout << "dW_dF_hand[" << i << "] = " << dW_dF_rev_rev[i] << std::endl; +// } + +// for (int i = 0; i < 81; ++i) { +// std::cout << "Comparing index " << i << ": FD = " << dW_dF_fwd_fd[i] +// << ", Hand derivative = " << dW_dF_rev_rev[i] << std::endl; +// EXPECT_NEAR(dW_dF_fwd_fd[i], dW_dF_rev_rev[i], 1e-6) << "Mismatch at index " << i; +// } +// } + +// TEST(StrainEnergyTest, StressFiniteDifferenceVsFwdRevAutodiff) { +// double E[9] = {0.1, 0.05, 0.02, 0.05, 0.2, 0.01, 0.02, 0.01, 0.15}; +// double mu = 1.0; +// double lambda = 1.0; +// double dW_dF_fwd_fd[81]; +// double dW_dF_fwd_rev[81]; + +// std::cout << "E: "; +// for (int i = 0; i < 9; ++i) std::cout << E[i] << " "; +// std::cout << std::endl; + +// std::cout << "Calling stress_FD..." << std::endl; +// second_deriv_fwd_FD(E, mu, lambda, dW_dF_fwd_fd); +// std::cout << "Finite difference stress computed:" << std::endl; +// for (int i = 0; i < 81; ++i) { +// std::cout << "dW_dF_fd[" << i << "] = " << dW_dF_fwd_fd[i] << std::endl; +// } + +// std::cout << "Calling hand_code_deriv (fwd_rev)..." << std::endl; +// second_deriv_fwd_rev(E, mu, lambda, dW_dF_fwd_rev); +// std::cout << "Hand coded stress computed (fwd_rev):" << std::endl; +// for (int i = 0; i < 81; ++i) { +// std::cout << "dW_dF_hand[" << i << "] = " << dW_dF_fwd_rev[i] << std::endl; +// } + +// for (int i = 0; i < 81; ++i) { +// std::cout << "Comparing index " << i << ": FD = " << dW_dF_fwd_fd[i] +// << ", Hand derivative = " << dW_dF_fwd_rev[i] << std::endl; +// EXPECT_NEAR(dW_dF_fwd_fd[i], dW_dF_fwd_rev[i], 1e-6) << "Mismatch at index " << i; +// } +// } + + +// TEST(StrainEnergyTest, StressFiniteDifferenceVsFwdRevAutodiff) { +// double E[9] = {0.12, -0.03, 0.01, -0.03, 0.08, 0.02, 0.01, 0.02, 0.11}; +// double mu = 1.0; +// double lambda = 1.0; +// double dW_dF_fwd_fd[81]; +// double dW_dF_fwd_rev[81]; + +// std::cout << "E: "; +// for (int i = 0; i < 9; ++i) std::cout << E[i] << " "; +// std::cout << std::endl; + +// std::cout << "Calling stress_FD..." << std::endl; +// second_deriv_fwd_FD(E, mu, lambda, dW_dF_fwd_fd); +// std::cout << "Finite difference stress computed:" << std::endl; +// for (int i = 0; i < 81; ++i) { +// std::cout << "dW_dF_fd[" << i << "] = " << dW_dF_fwd_fd[i] << std::endl; +// } + +// std::cout << "Calling hand_code_deriv (fwd_rev)..." << std::endl; +// second_deriv_hand_fwd(E, mu, lambda, dW_dF_fwd_rev); +// std::cout << "Hand coded stress computed (fwd_rev):" << std::endl; +// for (int i = 0; i < 81; ++i) { +// std::cout << "dW_dF_hand[" << i << "] = " << dW_dF_fwd_rev[i] << std::endl; +// } + +// for (int i = 0; i < 81; ++i) { +// std::cout << "Comparing index " << i << ": FD = " << dW_dF_fwd_fd[i] +// << ", Hand derivative = " << dW_dF_fwd_rev[i] << std::endl; +// EXPECT_NEAR(dW_dF_fwd_fd[i], dW_dF_fwd_rev[i], 1e-6) << "Mismatch at index " << i; +// } +// } + +TEST(StrainEnergyTest, StressFiniteDifferenceVsFwdRevAutodiff) { + double F[9] = {1.01, -0.03, 0.01, -0.03, 1.05, 0.02, 0.01, 0.02, 1.0}; + double E[9] = {0}; + double mu = 1.0; + double lambda = 1.0; + double* S = nullptr; + calc_E_from_F(F, E); + hand_code_deriv(E, mu, lambda, S); + std::cout << "{ "; + for (int i = 0; i < 9; ++i) { + std::cout << ", " << S[i]; + } + std::cout << " }" << std::endl; + std::cout << "E: { "; + for (int i = 0; i < 9; ++i) { + std::cout << ", " << E[i]; + } + std::cout << "}" << std::endl; + + double dW_dF_hand_fd[81]; + double dW_dF_hand_fwd[81]; + + std::cout << "E: "; + for (int i = 0; i < 9; ++i) std::cout << E[i] << " "; + std::cout << std::endl; + + std::cout << "Calling stress_FD..." << std::endl; + second_deriv_hand_FD(E, mu, lambda, dW_dF_hand_fd); + std::cout << "Finite difference stress computed:" << std::endl; + for (int i = 0; i < 81; ++i) { + std::cout << "dW_dF_hand_fd[" << i << "] = " << dW_dF_hand_fd[i] << std::endl; + } + + std::cout << "Calling hand_code_deriv (fwd_rev)..." << std::endl; + second_deriv_hand_fwd(E, mu, lambda, dW_dF_hand_fwd); + std::cout << "Hand coded stress computed (fwd_rev):" << std::endl; + for (int i = 0; i < 81; ++i) { + std::cout << "dW_dF_hand_fwd[" << i << "] = " << dW_dF_hand_fwd[i] << std::endl; + } + + for (int i = 0; i < 81; ++i) { + std::cout << "Comparing index " << i << ": FD = " << dW_dF_hand_fd[i] + << ", Hand derivative = " << dW_dF_hand_fwd[i] << std::endl; + EXPECT_NEAR(dW_dF_hand_fd[i], dW_dF_hand_fd[i], 1e-6) << "Mismatch at index " << i; + } +} \ No newline at end of file diff --git a/src/tests/CMakeLists.txt b/src/tests/CMakeLists.txt index de2941b2..e7bb5b97 100644 --- a/src/tests/CMakeLists.txt +++ b/src/tests/CMakeLists.txt @@ -174,6 +174,7 @@ if( TRIBOL_USE_ENZYME ) tribol_enzyme_nodal_normal.cpp tribol_enzyme_mortar_assembled.cpp tribol_enzyme_poly_intersect.cpp + enzyme_neohookean.cpp ) set(combined_test_depends tribol gtest) diff --git a/src/tests/enzyme_neohookean.cpp b/src/tests/enzyme_neohookean.cpp new file mode 100644 index 00000000..1309c8c1 --- /dev/null +++ b/src/tests/enzyme_neohookean.cpp @@ -0,0 +1,624 @@ +#include +#include +#include +#include "gtest/gtest.h" +#include +#include "axom/core.hpp" +#include "axom/slic/interface/slic.hpp" +#include +#include +#include "tribol/common/Enzyme.hpp" + +void multiply3x3( const double F[9], const double F_T[9], double C[9] ) +{ + for ( int i = 0; i < 3; ++i ) { // row of A + for ( int j = 0; j < 3; ++j ) { // column of B + C[i * 3 + j] = 0.0; + for ( int k = 0; k < 3; ++k ) { + C[i * 3 + j] += F_T[i * 3 + k] * F[k * 3 + j]; + } + } + } +} + +void calc_E_from_F( const double F[9], double E[9] ) +{ + double I[9] = { 1.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 1.0 }; + double F_T[9] = { 0.0 }; + F_T[0] = F[0]; + F_T[1] = F[3]; + F_T[2] = F[6]; + F_T[3] = F[1]; + F_T[4] = F[4]; + F_T[5] = F[7]; + F_T[6] = F[2]; + F_T[7] = F[5]; + F_T[8] = F[8]; + double C[9] = { 0.0 }; + multiply3x3( F, F_T, C ); + for ( int i = 0; i < 9; ++i ) { + E[i] = 0.5 * ( C[i] - I[i] ); + } +} + +// Calculates right cauchy stress tensor +void calc_cauchy_stress_tensor( const double E[9], double C[9] ) +{ + double I[9] = { 1, 0, 0, 0, 1, 0, 0, 0, 1 }; + for ( int i = 0; i < 9; ++i ) { + C[i] = 2 * E[i]; + C[i] += I[i]; + } +} + +// Calculates Trace +double calc_trace( const double C[9] ) +{ + double Tr_C = C[0] + C[4] + C[8]; + return Tr_C; +} + +double calc_determinant( const double C[9] ) +{ + double J = C[0] * ( C[4] * C[8] - C[5] * C[7] ) - C[1] * ( C[3] * C[8] - C[5] * C[6] ) + + C[2] * ( C[3] * C[7] - C[4] * C[6] ); + J = std::sqrt( J ); + return J; +} + +bool invert3x3( const double F[9], double Finv[9] ) +{ + // Compute the determinant + double det = F[0] * ( F[4] * F[8] - F[5] * F[7] ) - F[1] * ( F[3] * F[8] - F[5] * F[6] ) + + F[2] * ( F[3] * F[7] - F[4] * F[6] ); + + if ( std::abs( det ) < std::numeric_limits::epsilon() ) return false; // Singular matrix + + double invDet = 1.0 / det; + + // Compute the inverse using the formula for the inverse of a 3x3 matrix + Finv[0] = ( F[4] * F[8] - F[5] * F[7] ) * invDet; + Finv[1] = -( F[1] * F[8] - F[2] * F[7] ) * invDet; + Finv[2] = ( F[1] * F[5] - F[2] * F[4] ) * invDet; + + Finv[3] = -( F[3] * F[8] - F[5] * F[6] ) * invDet; + Finv[4] = ( F[0] * F[8] - F[2] * F[6] ) * invDet; + Finv[5] = -( F[0] * F[5] - F[2] * F[3] ) * invDet; + + Finv[6] = ( F[3] * F[7] - F[4] * F[6] ) * invDet; + Finv[7] = -( F[0] * F[7] - F[1] * F[6] ) * invDet; + Finv[8] = ( F[0] * F[4] - F[1] * F[3] ) * invDet; + + return true; +} + +// build strain energy equation +void strain_energy( double* E, double mu, double lambda, double* W ) +{ + double C[9] = { 0.0 }; + calc_cauchy_stress_tensor( E, C ); + double Tr_C = calc_trace( C ); + double J = calc_determinant( C ); + *W = mu / 2.0 * ( Tr_C - 3.0 ) - mu * log( J ) + lambda / 2.0 * pow( ( log( J ) ), 2.0 ); +} + +// calc stress using enzyme fwddiff +void stress( double* E, double mu, double lambda, double* dW_dE ) +{ + double W = 0.0; + for ( int i = 0; i < 9; ++i ) { + double dE[9] = { 0.0 }; + dE[i] = 1.0; + double dmu = 0.0; + double dlambda = 0.0; + double dw = 0.0; + __enzyme_fwddiff( (void*)strain_energy, E, dE, mu, dmu, lambda, dlambda, &W, &dw ); + dW_dE[i] = dw; + } +} + +// calc stress using enzyme autodiff +void stress_reverse( double* E, double mu, double lambda, double* dW_dE ) +{ + double dE[9] = { 0.0 }; + double W = 0.0; + double dW = 1.0; + __enzyme_autodiff( strain_energy, enzyme_dup, E, dE, enzyme_const, mu, enzyme_const, lambda, enzyme_dup, &W, + &dW ); + + for ( int i = 0; i < 9; ++i ) { + dW_dE[i] = dE[i]; + } +} + +// calc stress using finite Difference +void stress_FD( double* E, double mu, double lambda, double* dW_dE, double h = 1e-7 ) +{ + double E_plus[9] = { 0.0 }; + double E_minus[9] = { 0.0 }; + double W_plus; + double W_minus; + for ( int i = 0; i < 9; ++i ) { + for ( int j = 0; j < 9; ++j ) { + E_plus[j] = E[j]; + E_minus[j] = E[j]; + } + E_plus[i] = E[i] + h; + E_minus[i] = E[i] - h; + strain_energy( E_plus, mu, lambda, &W_plus ); + strain_energy( E_minus, mu, lambda, &W_minus ); + dW_dE[i] = ( W_plus - W_minus ) / ( 2 * h ); + } +} + +void hand_code_deriv( double* E, double mu, double lambda, double* S ) +{ + double C[9] = { 0.0 }; + double Cinv[9]; + double I[9] = { 1.0, 0.0, 0.0, 0, 1.0, 0.0, 0.0, 0.0, 1.0 }; + calc_cauchy_stress_tensor( E, C ); + double J = calc_determinant( C ); + invert3x3( C, Cinv ); + double first_term[9]; + for ( int i = 0; i < 9; ++i ) { + first_term[i] = lambda * std::log( J ) * Cinv[i]; + } + double second_term[9]; + for ( int i = 0; i < 9; ++i ) { + second_term[i] = I[i] - Cinv[i]; + second_term[i] *= mu; + } + for ( int i = 0; i < 9; ++i ) { + S[i] = first_term[i] + second_term[i]; + } +} + +void second_deriv_fwd_fwd( double* E, double mu, double lambda, double* d2W_d2E ) +{ + double dW[9] = { 0.0 }; + double d2w[9] = { 0.0 }; + + for ( int i = 0; i < 9; ++i ) { + double d2E[9] = { 0.0 }; + d2E[i] = 1.0; + double d2mu = 0.0; + double d2lambda = 0.0; + __enzyme_fwddiff( (void*)stress, E, d2E, mu, d2mu, lambda, d2lambda, &dW, &d2w ); + for ( int j = 0; j < 9; ++j ) { + d2W_d2E[9 * i + j] = d2w[j]; + } + } +} + +void second_deriv_rev_fwd( double* E, double mu, double lambda, double* d2W_d2E ) +{ + double dW[9] = { 0.0 }; + double d2w[9] = { 0.0 }; + + for ( int i = 0; i < 9; ++i ) { + double d2E[9] = { 0.0 }; + d2E[i] = 1.0; + double d2mu = 0.0; + double d2lambda = 0.0; + __enzyme_fwddiff( (void*)stress_reverse, E, d2E, mu, d2mu, lambda, d2lambda, &dW, &d2w ); + for ( int j = 0; j < 9; ++j ) { + d2W_d2E[9 * i + j] = d2w[j]; + } + } +} + +void second_deriv_rev_rev( double* E, double mu, double lambda, double* d2W_d2E ) +{ + for ( int i = 0; i < 9; ++i ) { + double d2E[81] = { 0.0 }; + double W[9] = { 0.0 }; + double d2W[9] = { 0.0 }; + d2W[i] = 1.0; + __enzyme_autodiff( stress_reverse, enzyme_dup, E, d2E, enzyme_const, mu, enzyme_const, lambda, enzyme_dup, &W, + &d2W ); + for ( int j = 0; j < 9; ++j ) { + d2W_d2E[9 * j + i] = d2E[j]; + } + } +} + +void second_deriv_fwd_rev( double* E, double mu, double lambda, double* d2W_d2E ) +{ + for ( int i = 0; i < 9; ++i ) { + double d2E[81] = { 0.0 }; + double W[9] = { 0.0 }; + double d2W[9] = { 0.0 }; + d2W[i] = 1.0; + __enzyme_autodiff( stress, enzyme_dup, E, d2E, enzyme_const, mu, enzyme_const, lambda, enzyme_dup, &W, &d2W ); + for ( int j = 0; j < 9; ++j ) { + d2W_d2E[9 * i + j] = d2E[j]; + } + } +} + +void second_deriv_fwd_FD( double* E, double mu, double lambda, double* d2W_d2E, double h = 1e-7 ) +{ + double E_plus[9] = { 0.0 }; + double E_minus[9] = { 0.0 }; + double dW_plus[9] = { 0.0 }; + double dW_minus[9] = { 0.0 }; + for ( int i = 0; i < 9; ++i ) { + for ( int j = 0; j < 9; ++j ) { + E_plus[j] = E[j]; + E_minus[j] = E[j]; + } + E_plus[i] = E[i] + h; + E_minus[i] = E[i] - h; + stress( E_plus, mu, lambda, dW_plus ); + stress( E_minus, mu, lambda, dW_minus ); + for ( int j = 0; j < 9; ++j ) { + d2W_d2E[9 * i + j] = ( dW_plus[j] - dW_minus[j] ) / ( 2 * h ); + } + } +} + +void second_deriv_hand_fwd( double* E, double mu, double lambda, double* d2W_d2E ) +{ + double dW[9] = { 0.0 }; + double d2w[9] = { 0.0 }; + + for ( int i = 0; i < 9; ++i ) { + double d2E[9] = { 0.0 }; + d2E[i] = 1.0; + double d2mu = 0.0; + double d2lambda = 0.0; + __enzyme_fwddiff( (void*)hand_code_deriv, E, d2E, mu, d2mu, lambda, d2lambda, &dW, &d2w ); + for ( int j = 0; j < 9; ++j ) { + d2W_d2E[9 * i + j] = d2w[j]; + } + } +} + +void second_deriv_hand_FD( double* E, double mu, double lambda, double* d2W_d2E, double h = 1e-7 ) +{ + double E_plus[9] = { 0.0 }; + double E_minus[9] = { 0.0 }; + double dw_minus[9] = { 0.0 }; + double dw_plus[9] = { 0.0 }; + for ( int i = 0; i < 9; ++i ) { + for ( int j = 0; j < 9; ++j ) { + E_plus[j] = E[j]; + E_minus[j] = E[j]; + } + E_plus[i] = E[i] + h; + E_minus[i] = E[i] - h; + hand_code_deriv( E_plus, mu, lambda, dw_plus ); + hand_code_deriv( E_minus, mu, lambda, dw_minus ); + for ( int j = 0; j < 9; ++j ) { + d2W_d2E[9 * i + j] = ( dw_plus[j] - dw_minus[j] ) / ( 2 * h ); + } + } +} + +class EnzymeNeohookeanTest : public testing::Test { + protected: + double F_all[27] = { 1.0, 0.5, 0.0, 0.0, 1.2, 0.1, 0.0, 0.0, 1.0, 1.1, 0.7, 0.1, 0.2, 1.4, + 0.2, 0.02, 0.04, 1.05, 1.05, 0.2, 0.04, 0.0, 1.1, 0.14, 0.07, 0.0, 1.09 }; + + double mu = 1.6; + double lambda = 1.2; + + double E_all[27] = { 0.0 }; + + double tol = 1e-6; + + void SetUp() override + { + for ( int i = 0; i < 3; ++i ) { + calc_E_from_F( F_all + i * 9, E_all + i * 9 ); + } + } +}; + +TEST_F( EnzymeNeohookeanTest, finite_difference_vs_exact ) +{ + for ( int i = 0; i < 3; ++i ) { + double S_FD[9] = { 0.0 }; + double S_exact[9] = { 0.0 }; + stress_FD( E_all + i * 9, mu, lambda, S_FD ); + hand_code_deriv( E_all + i * 9, mu, lambda, S_exact ); + for ( int j = 0; j < 9; ++j ) { + EXPECT_NEAR( S_FD[j], S_exact[j], tol ); + } + } +} + +// TEST(StrainEnergyTest, StressFiniteDifferenceVsAutodiff) { +// double F[9] = {1.0, 0.5, 0.0, 0.0, 1.2, 0.1, 0.0, 0.0, 1.0}; +// double E[9] = {0.0}; +// double mu = 1.0; +// double lambda = 1.0; +// double dW_dF_fd[9]; +// double dW_dF_ad[9]; +// double J = calc_determinant(F); +// calc_green_lagrange(F, E); +// std::cout << "Calling stress_FD..." << std::endl; +// stress_FD(E, J, mu, lambda, dW_dF_fd); +// std::cout << "Finite difference stress computed:" << std::endl; +// for (int i = 0; i < 9; ++i) { +// std::cout << "dW_dF_fd[" << i << "] = " << dW_dF_fd[i] << std::endl; +// } + +// std::cout << "Calling stress..." << std::endl; +// stress(E, J, mu, lambda, dW_dF_ad); +// std::cout << "Autodiff stress computed:" << std::endl; +// for (int i = 0; i < 9; ++i) { +// std::cout << "dW_dF_ad[" << i << "] = " << dW_dF_ad[i] << std::endl; +// } + +// for (int i = 0; i < 9; ++i) { +// std::cout << "Comparing index " << i << ": FD = " << dW_dF_fd[i] +// << ", AD = " << dW_dF_ad[i] << std::endl; +// EXPECT_NEAR(dW_dF_fd[i], dW_dF_ad[i], 1e-6) << "Mismatch at index " << i; +// } +// } + +// TEST(StrainEnergyTest, StressFiniteDifferenceVsReverseAutodiff) { +// double E[9] = {1.0, 0.5, 0.0, 0.0, 1.2, 0.1, 0.0, 0.0, 1.0}; +// double mu = 1.0; +// double lambda = 1.0; +// double dW_dF_fd[9]; +// double dW_dF_rev[9]; +// std::cout << "E: "; +// for (int i = 0; i < 9; ++i) std::cout << E[i] << " "; +// std::cout << std::endl; +// std::cout << "Calling stress_FD..." << std::endl; +// stress_FD(E, mu, lambda, dW_dF_fd); +// std::cout << "Finite difference stress computed:" << std::endl; +// for (int i = 0; i < 9; ++i) { +// std::cout << "dW_dF_fd[" << i << "] = " << dW_dF_fd[i] << std::endl; +// } + +// std::cout << "Calling stress_reverse..." << std::endl; +// stress_reverse(E, mu, lambda, dW_dF_rev); +// std::cout << "Reverse-mode autodiff stress computed:" << std::endl; +// for (int i = 0; i < 9; ++i) { +// std::cout << "dW_dF_rev[" << i << "] = " << dW_dF_rev[i] << std::endl; +// } + +// for (int i = 0; i < 9; ++i) { +// std::cout << "Comparing index " << i << ": FD = " << dW_dF_fd[i] +// << ", Reverse AD = " << dW_dF_rev[i] << std::endl; +// EXPECT_NEAR(dW_dF_fd[i], dW_dF_rev[i], 1e-6) << "Mismatch at index " << i; +// } +// } + +// TEST(StrainEnergyTest, StressFiniteDifferenceVsReverseAutodiff) { +// double E[9] = {0.1, 0.05, 0.02, 0.05, 0.2, 0.01, 0.02, 0.01, 0.15}; +// double mu = 1.0; +// double lambda = 1.0; +// double dW_dF_fd[9]; +// double dW_dF_hand[9]; +// std::cout << "E: "; +// for (int i = 0; i < 9; ++i) std::cout << E[i] << " "; +// std::cout << std::endl; +// std::cout << "Calling stress_FD..." << std::endl; +// stress_FD(E, mu, lambda, dW_dF_fd); +// std::cout << "Finite difference stress computed:" << std::endl; +// for (int i = 0; i < 9; ++i) { +// std::cout << "dW_dF_fd[" << i << "] = " << dW_dF_fd[i] << std::endl; +// } + +// std::cout << "Calling hand_code_deriv..." << std::endl; +// hand_code_deriv(E, mu, lambda, dW_dF_hand); +// std::cout << "Hand coded stress computed:" << std::endl; +// for (int i = 0; i < 9; ++i) { +// std::cout << "dW_dF_hand[" << i << "] = " << dW_dF_hand[i] << std::endl; +// } + +// for (int i = 0; i < 9; ++i) { +// std::cout << "Comparing index " << i << ": FD = " << dW_dF_fd[i] +// << ", Hand derivative = " << dW_dF_hand[i] << std::endl; +// EXPECT_NEAR(dW_dF_fd[i], dW_dF_hand[i], 1e-6) << "Mismatch at index " << i; +// } +// } + +// TEST(StrainEnergyTest, StressFiniteDifferenceVsReverseAutodiff) { +// double E[9] = {0.1, 0.05, 0.02, 0.05, 0.2, 0.01, 0.02, 0.01, 0.15}; +// double mu = 1.0; +// double lambda = 1.0; +// double dW_dF_fwd_fd[81]; +// double dW_dF_fwd_fwd[81]; +// std::cout << "E: "; +// for (int i = 0; i < 81; ++i) std::cout << E[i] << " "; +// std::cout << std::endl; +// std::cout << "Calling stress_FD..." << std::endl; +// second_deriv_fwd_FD(E, mu, lambda, dW_dF_fwd_fd); +// std::cout << "Finite difference stress computed:" << std::endl; +// for (int i = 0; i < 81; ++i) { +// std::cout << "dW_dF_fd[" << i << "] = " << dW_dF_fwd_fd[i] << std::endl; +// } + +// std::cout << "Calling hand_code_deriv..." << std::endl; +// second_deriv_fwd_fwd(E, mu, lambda, dW_dF_fwd_fwd); +// std::cout << "Hand coded stress computed:" << std::endl; +// for (int i = 0; i < 81; ++i) { +// std::cout << "dW_dF_hand[" << i << "] = " << dW_dF_fwd_fwd[i] << std::endl; +// } + +// for (int i = 0; i < 81; ++i) { +// std::cout << "Comparing index " << i << ": FD = " << dW_dF_fwd_fd[i] +// << ", Hand derivative = " << dW_dF_fwd_fwd[i] << std::endl; +// EXPECT_NEAR(dW_dF_fwd_fd[i], dW_dF_fwd_fwd[i], 1e-6) << "Mismatch at index " << i; +// } +// } + +// TEST(StrainEnergyTest, StressFiniteDifferenceVsReverseAutodiff) { +// double E[9] = {0.1, 0.05, 0.02, 0.05, 0.2, 0.01, 0.02, 0.01, 0.15}; +// double mu = 1.0; +// double lambda = 1.0; +// double dW_dF_fwd_fd[81]; +// double dW_dF_fwd_rev[81]; + +// std::cout << "E: "; +// for (int i = 0; i < 9; ++i) std::cout << E[i] << " "; +// std::cout << std::endl; + +// std::cout << "Calling stress_FD..." << std::endl; +// second_deriv_fwd_FD(E, mu, lambda, dW_dF_fwd_fd); +// std::cout << "Finite difference stress computed:" << std::endl; +// for (int i = 0; i < 81; ++i) { +// std::cout << "dW_dF_fd[" << i << "] = " << dW_dF_fwd_fd[i] << std::endl; +// } + +// std::cout << "Calling hand_code_deriv (fwd_rev)..." << std::endl; +// second_deriv_fwd_rev(E, mu, lambda, dW_dF_fwd_rev); +// std::cout << "Hand coded stress computed (fwd_rev):" << std::endl; +// for (int i = 0; i < 81; ++i) { +// std::cout << "dW_dF_hand[" << i << "] = " << dW_dF_fwd_rev[i] << std::endl; +// } + +// for (int i = 0; i < 81; ++i) { +// std::cout << "Comparing index " << i << ": FD = " << dW_dF_fwd_fd[i] +// << ", Hand derivative = " << dW_dF_fwd_rev[i] << std::endl; +// EXPECT_NEAR(dW_dF_fwd_fd[i], dW_dF_fwd_rev[i], 1e-6) << "Mismatch at index " << i; +// } +// } + +// TEST(StrainEnergyTest, StressFiniteDifferenceVsReverseAutodiff) { +// double E[9] = {0.1, 0.05, 0.02, 0.05, 0.2, 0.01, 0.02, 0.01, 0.15}; +// double mu = 1.0; +// double lambda = 1.0; +// double dW_dF_fwd_fd[81]; +// double dW_dF_rev_rev[81]; + +// std::cout << "E: "; +// for (int i = 0; i < 9; ++i) std::cout << E[i] << " "; +// std::cout << std::endl; + +// std::cout << "Calling stress_FD..." << std::endl; +// second_deriv_fwd_FD(E, mu, lambda, dW_dF_fwd_fd); +// std::cout << "Finite difference stress computed:" << std::endl; +// for (int i = 0; i < 81; ++i) { +// std::cout << "dW_dF_fd[" << i << "] = " << dW_dF_fwd_fd[i] << std::endl; +// } + +// std::cout << "Calling hand_code_deriv (rev_rev)..." << std::endl; +// second_deriv_rev_rev(E, mu, lambda, dW_dF_rev_rev); +// std::cout << "Hand coded stress computed (rev_rev):" << std::endl; +// for (int i = 0; i < 81; ++i) { +// std::cout << "dW_dF_hand[" << i << "] = " << dW_dF_rev_rev[i] << std::endl; +// } + +// for (int i = 0; i < 81; ++i) { +// std::cout << "Comparing index " << i << ": FD = " << dW_dF_fwd_fd[i] +// << ", Hand derivative = " << dW_dF_rev_rev[i] << std::endl; +// EXPECT_NEAR(dW_dF_fwd_fd[i], dW_dF_rev_rev[i], 1e-6) << "Mismatch at index " << i; +// } +// } + +// TEST(StrainEnergyTest, StressFiniteDifferenceVsFwdRevAutodiff) { +// double E[9] = {0.1, 0.05, 0.02, 0.05, 0.2, 0.01, 0.02, 0.01, 0.15}; +// double mu = 1.0; +// double lambda = 1.0; +// double dW_dF_fwd_fd[81]; +// double dW_dF_fwd_rev[81]; + +// std::cout << "E: "; +// for (int i = 0; i < 9; ++i) std::cout << E[i] << " "; +// std::cout << std::endl; + +// std::cout << "Calling stress_FD..." << std::endl; +// second_deriv_fwd_FD(E, mu, lambda, dW_dF_fwd_fd); +// std::cout << "Finite difference stress computed:" << std::endl; +// for (int i = 0; i < 81; ++i) { +// std::cout << "dW_dF_fd[" << i << "] = " << dW_dF_fwd_fd[i] << std::endl; +// } + +// std::cout << "Calling hand_code_deriv (fwd_rev)..." << std::endl; +// second_deriv_fwd_rev(E, mu, lambda, dW_dF_fwd_rev); +// std::cout << "Hand coded stress computed (fwd_rev):" << std::endl; +// for (int i = 0; i < 81; ++i) { +// std::cout << "dW_dF_hand[" << i << "] = " << dW_dF_fwd_rev[i] << std::endl; +// } + +// for (int i = 0; i < 81; ++i) { +// std::cout << "Comparing index " << i << ": FD = " << dW_dF_fwd_fd[i] +// << ", Hand derivative = " << dW_dF_fwd_rev[i] << std::endl; +// EXPECT_NEAR(dW_dF_fwd_fd[i], dW_dF_fwd_rev[i], 1e-6) << "Mismatch at index " << i; +// } +// } + +// TEST(StrainEnergyTest, StressFiniteDifferenceVsFwdRevAutodiff) { +// double E[9] = {0.12, -0.03, 0.01, -0.03, 0.08, 0.02, 0.01, 0.02, 0.11}; +// double mu = 1.0; +// double lambda = 1.0; +// double dW_dF_fwd_fd[81]; +// double dW_dF_fwd_rev[81]; + +// std::cout << "E: "; +// for (int i = 0; i < 9; ++i) std::cout << E[i] << " "; +// std::cout << std::endl; + +// std::cout << "Calling stress_FD..." << std::endl; +// second_deriv_fwd_FD(E, mu, lambda, dW_dF_fwd_fd); +// std::cout << "Finite difference stress computed:" << std::endl; +// for (int i = 0; i < 81; ++i) { +// std::cout << "dW_dF_fd[" << i << "] = " << dW_dF_fwd_fd[i] << std::endl; +// } + +// std::cout << "Calling hand_code_deriv (fwd_rev)..." << std::endl; +// second_deriv_hand_fwd(E, mu, lambda, dW_dF_fwd_rev); +// std::cout << "Hand coded stress computed (fwd_rev):" << std::endl; +// for (int i = 0; i < 81; ++i) { +// std::cout << "dW_dF_hand[" << i << "] = " << dW_dF_fwd_rev[i] << std::endl; +// } + +// for (int i = 0; i < 81; ++i) { +// std::cout << "Comparing index " << i << ": FD = " << dW_dF_fwd_fd[i] +// << ", Hand derivative = " << dW_dF_fwd_rev[i] << std::endl; +// EXPECT_NEAR(dW_dF_fwd_fd[i], dW_dF_fwd_rev[i], 1e-6) << "Mismatch at index " << i; +// } +// } + +// TEST(StrainEnergyTest, StressFiniteDifferenceVsFwdRevAutodiff) { +// double F[9] = {1.01, -0.03, 0.01, -0.03, 1.05, 0.02, 0.01, 0.02, 1.0}; +// double E[9] = {0}; +// double mu = 1.0; +// double lambda = 1.0; +// double* S = nullptr; +// calc_E_from_F(F, E); +// hand_code_deriv(E, mu, lambda, S); +// std::cout << "{ "; +// for (int i = 0; i < 9; ++i) { +// std::cout << ", " << S[i]; +// } +// std::cout << " }" << std::endl; +// std::cout << "E: { "; +// for (int i = 0; i < 9; ++i) { +// std::cout << ", " << E[i]; +// } +// std::cout << "}" << std::endl; + +// double dW_dF_hand_fd[81]; +// double dW_dF_hand_fwd[81]; + +// std::cout << "E: "; +// for (int i = 0; i < 9; ++i) std::cout << E[i] << " "; +// std::cout << std::endl; + +// std::cout << "Calling stress_FD..." << std::endl; +// second_deriv_hand_FD(E, mu, lambda, dW_dF_hand_fd); +// std::cout << "Finite difference stress computed:" << std::endl; +// for (int i = 0; i < 81; ++i) { +// std::cout << "dW_dF_hand_fd[" << i << "] = " << dW_dF_hand_fd[i] << std::endl; +// } + +// std::cout << "Calling hand_code_deriv (fwd_rev)..." << std::endl; +// second_deriv_hand_fwd(E, mu, lambda, dW_dF_hand_fwd); +// std::cout << "Hand coded stress computed (fwd_rev):" << std::endl; +// for (int i = 0; i < 81; ++i) { +// std::cout << "dW_dF_hand_fwd[" << i << "] = " << dW_dF_hand_fwd[i] << std::endl; +// } + +// for (int i = 0; i < 81; ++i) { +// std::cout << "Comparing index " << i << ": FD = " << dW_dF_hand_fd[i] +// << ", Hand derivative = " << dW_dF_hand_fwd[i] << std::endl; +// EXPECT_NEAR(dW_dF_hand_fd[i], dW_dF_hand_fd[i], 1e-6) << "Mismatch at index " << i; +// } +// } \ No newline at end of file From 86a588fa9505cbac65b617568c79b67a1a3ce4e4 Mon Sep 17 00:00:00 2001 From: Ryan Lutz Date: Tue, 17 Jun 2025 16:13:33 -0700 Subject: [PATCH 2/4] enzyme_neohookean test and new hpp file --- src/tests/enzyme_neohookean.cpp | 639 ++++---------------------------- src/tests/neohookean_common.hpp | 293 +++++++++++++++ 2 files changed, 368 insertions(+), 564 deletions(-) create mode 100644 src/tests/neohookean_common.hpp diff --git a/src/tests/enzyme_neohookean.cpp b/src/tests/enzyme_neohookean.cpp index 1309c8c1..558490fa 100644 --- a/src/tests/enzyme_neohookean.cpp +++ b/src/tests/enzyme_neohookean.cpp @@ -1,624 +1,135 @@ #include +#include +#include "neohookean_common.hpp" #include #include -#include "gtest/gtest.h" #include -#include "axom/core.hpp" -#include "axom/slic/interface/slic.hpp" #include #include #include "tribol/common/Enzyme.hpp" -void multiply3x3( const double F[9], const double F_T[9], double C[9] ) -{ - for ( int i = 0; i < 3; ++i ) { // row of A - for ( int j = 0; j < 3; ++j ) { // column of B - C[i * 3 + j] = 0.0; - for ( int k = 0; k < 3; ++k ) { - C[i * 3 + j] += F_T[i * 3 + k] * F[k * 3 + j]; - } - } - } -} - -void calc_E_from_F( const double F[9], double E[9] ) -{ - double I[9] = { 1.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 1.0 }; - double F_T[9] = { 0.0 }; - F_T[0] = F[0]; - F_T[1] = F[3]; - F_T[2] = F[6]; - F_T[3] = F[1]; - F_T[4] = F[4]; - F_T[5] = F[7]; - F_T[6] = F[2]; - F_T[7] = F[5]; - F_T[8] = F[8]; - double C[9] = { 0.0 }; - multiply3x3( F, F_T, C ); - for ( int i = 0; i < 9; ++i ) { - E[i] = 0.5 * ( C[i] - I[i] ); - } -} - -// Calculates right cauchy stress tensor -void calc_cauchy_stress_tensor( const double E[9], double C[9] ) -{ - double I[9] = { 1, 0, 0, 0, 1, 0, 0, 0, 1 }; - for ( int i = 0; i < 9; ++i ) { - C[i] = 2 * E[i]; - C[i] += I[i]; - } -} - -// Calculates Trace -double calc_trace( const double C[9] ) -{ - double Tr_C = C[0] + C[4] + C[8]; - return Tr_C; -} - -double calc_determinant( const double C[9] ) -{ - double J = C[0] * ( C[4] * C[8] - C[5] * C[7] ) - C[1] * ( C[3] * C[8] - C[5] * C[6] ) + - C[2] * ( C[3] * C[7] - C[4] * C[6] ); - J = std::sqrt( J ); - return J; -} - -bool invert3x3( const double F[9], double Finv[9] ) -{ - // Compute the determinant - double det = F[0] * ( F[4] * F[8] - F[5] * F[7] ) - F[1] * ( F[3] * F[8] - F[5] * F[6] ) + - F[2] * ( F[3] * F[7] - F[4] * F[6] ); - - if ( std::abs( det ) < std::numeric_limits::epsilon() ) return false; // Singular matrix - - double invDet = 1.0 / det; - - // Compute the inverse using the formula for the inverse of a 3x3 matrix - Finv[0] = ( F[4] * F[8] - F[5] * F[7] ) * invDet; - Finv[1] = -( F[1] * F[8] - F[2] * F[7] ) * invDet; - Finv[2] = ( F[1] * F[5] - F[2] * F[4] ) * invDet; - - Finv[3] = -( F[3] * F[8] - F[5] * F[6] ) * invDet; - Finv[4] = ( F[0] * F[8] - F[2] * F[6] ) * invDet; - Finv[5] = -( F[0] * F[5] - F[2] * F[3] ) * invDet; - Finv[6] = ( F[3] * F[7] - F[4] * F[6] ) * invDet; - Finv[7] = -( F[0] * F[7] - F[1] * F[6] ) * invDet; - Finv[8] = ( F[0] * F[4] - F[1] * F[3] ) * invDet; - - return true; -} - -// build strain energy equation -void strain_energy( double* E, double mu, double lambda, double* W ) -{ - double C[9] = { 0.0 }; - calc_cauchy_stress_tensor( E, C ); - double Tr_C = calc_trace( C ); - double J = calc_determinant( C ); - *W = mu / 2.0 * ( Tr_C - 3.0 ) - mu * log( J ) + lambda / 2.0 * pow( ( log( J ) ), 2.0 ); -} +class EnzymeNeohookeanTest : public testing::Test { + protected: + double F_all[27] = { 1.0, 0.5, 0.0, 0.0, 1.2, 0.1, 0.0, 0.0, 1.0, 1.1, 0.7, 0.1, 0.2, 1.4, + 0.2, 0.02, 0.04, 1.05, 1.05, 0.2, 0.04, 0.0, 1.1, 0.14, 0.07, 0.0, 1.09 }; -// calc stress using enzyme fwddiff -void stress( double* E, double mu, double lambda, double* dW_dE ) -{ - double W = 0.0; - for ( int i = 0; i < 9; ++i ) { - double dE[9] = { 0.0 }; - dE[i] = 1.0; - double dmu = 0.0; - double dlambda = 0.0; - double dw = 0.0; - __enzyme_fwddiff( (void*)strain_energy, E, dE, mu, dmu, lambda, dlambda, &W, &dw ); - dW_dE[i] = dw; - } -} + double mu = 1.6; + double lambda = 1.2; -// calc stress using enzyme autodiff -void stress_reverse( double* E, double mu, double lambda, double* dW_dE ) -{ - double dE[9] = { 0.0 }; - double W = 0.0; - double dW = 1.0; - __enzyme_autodiff( strain_energy, enzyme_dup, E, dE, enzyme_const, mu, enzyme_const, lambda, enzyme_dup, &W, - &dW ); + double E_all[27] = { 0.0 }; - for ( int i = 0; i < 9; ++i ) { - dW_dE[i] = dE[i]; - } -} + double tol = 1e-6; -// calc stress using finite Difference -void stress_FD( double* E, double mu, double lambda, double* dW_dE, double h = 1e-7 ) -{ - double E_plus[9] = { 0.0 }; - double E_minus[9] = { 0.0 }; - double W_plus; - double W_minus; - for ( int i = 0; i < 9; ++i ) { - for ( int j = 0; j < 9; ++j ) { - E_plus[j] = E[j]; - E_minus[j] = E[j]; + void SetUp() override + { + for ( int i = 0; i < 3; ++i ) { + calc_E_from_F( F_all + i * 9, E_all + i * 9 ); } - E_plus[i] = E[i] + h; - E_minus[i] = E[i] - h; - strain_energy( E_plus, mu, lambda, &W_plus ); - strain_energy( E_minus, mu, lambda, &W_minus ); - dW_dE[i] = ( W_plus - W_minus ) / ( 2 * h ); - } -} - -void hand_code_deriv( double* E, double mu, double lambda, double* S ) -{ - double C[9] = { 0.0 }; - double Cinv[9]; - double I[9] = { 1.0, 0.0, 0.0, 0, 1.0, 0.0, 0.0, 0.0, 1.0 }; - calc_cauchy_stress_tensor( E, C ); - double J = calc_determinant( C ); - invert3x3( C, Cinv ); - double first_term[9]; - for ( int i = 0; i < 9; ++i ) { - first_term[i] = lambda * std::log( J ) * Cinv[i]; } - double second_term[9]; - for ( int i = 0; i < 9; ++i ) { - second_term[i] = I[i] - Cinv[i]; - second_term[i] *= mu; - } - for ( int i = 0; i < 9; ++i ) { - S[i] = first_term[i] + second_term[i]; - } -} +}; -void second_deriv_fwd_fwd( double* E, double mu, double lambda, double* d2W_d2E ) +TEST_F( EnzymeNeohookeanTest, finite_difference_vs_exact ) { - double dW[9] = { 0.0 }; - double d2w[9] = { 0.0 }; - - for ( int i = 0; i < 9; ++i ) { - double d2E[9] = { 0.0 }; - d2E[i] = 1.0; - double d2mu = 0.0; - double d2lambda = 0.0; - __enzyme_fwddiff( (void*)stress, E, d2E, mu, d2mu, lambda, d2lambda, &dW, &d2w ); + for ( int i = 0; i < 3; ++i ) { + double S_FD[9] = { 0.0 }; + double S_exact[9] = { 0.0 }; + stress_FD( E_all + i * 9, mu, lambda, S_FD ); + hand_code_deriv( E_all + i * 9, mu, lambda, S_exact ); for ( int j = 0; j < 9; ++j ) { - d2W_d2E[9 * i + j] = d2w[j]; + EXPECT_NEAR( S_FD[j], S_exact[j], tol ); } } } -void second_deriv_rev_fwd( double* E, double mu, double lambda, double* d2W_d2E ) +TEST_F( EnzymeNeohookeanTest, finite_difference_vs_fwd_diff ) { - double dW[9] = { 0.0 }; - double d2w[9] = { 0.0 }; - - for ( int i = 0; i < 9; ++i ) { - double d2E[9] = { 0.0 }; - d2E[i] = 1.0; - double d2mu = 0.0; - double d2lambda = 0.0; - __enzyme_fwddiff( (void*)stress_reverse, E, d2E, mu, d2mu, lambda, d2lambda, &dW, &d2w ); + for ( int i = 0; i < 3; ++i ) { + double S_FD[9] = { 0.0 }; + double S_fwddiff[9] = { 0.0 }; + stress_FD( E_all + i * 9, mu, lambda, S_FD ); + stress( E_all + i * 9, mu, lambda, S_fwddiff ); for ( int j = 0; j < 9; ++j ) { - d2W_d2E[9 * i + j] = d2w[j]; + EXPECT_NEAR( S_FD[j], S_fwddiff[j], tol ); } } } -void second_deriv_rev_rev( double* E, double mu, double lambda, double* d2W_d2E ) +TEST_F( EnzymeNeohookeanTest, finite_difference_vs_autodiff ) { - for ( int i = 0; i < 9; ++i ) { - double d2E[81] = { 0.0 }; - double W[9] = { 0.0 }; - double d2W[9] = { 0.0 }; - d2W[i] = 1.0; - __enzyme_autodiff( stress_reverse, enzyme_dup, E, d2E, enzyme_const, mu, enzyme_const, lambda, enzyme_dup, &W, - &d2W ); + for ( int i = 0; i < 3; ++i ) { + double S_FD[9] = { 0.0 }; + double S_autodiff[9] = { 0.0 }; + stress_FD( E_all + i * 9, mu, lambda, S_FD ); + stress_reverse( E_all + i * 9, mu, lambda, S_autodiff ); for ( int j = 0; j < 9; ++j ) { - d2W_d2E[9 * j + i] = d2E[j]; + EXPECT_NEAR( S_FD[j], S_autodiff[j], tol ); } } } -void second_deriv_fwd_rev( double* E, double mu, double lambda, double* d2W_d2E ) +TEST_F( EnzymeNeohookeanTest, 2nd_deriv_finite_difference_vs_fwd_fwd ) { - for ( int i = 0; i < 9; ++i ) { - double d2E[81] = { 0.0 }; - double W[9] = { 0.0 }; - double d2W[9] = { 0.0 }; - d2W[i] = 1.0; - __enzyme_autodiff( stress, enzyme_dup, E, d2E, enzyme_const, mu, enzyme_const, lambda, enzyme_dup, &W, &d2W ); + for ( int i = 0; i < 3; ++i ) { + double dW2_FD[81] = { 0.0 }; + double dW2_fwd_fwd[81] = { 0.0 }; + second_deriv_fwd_FD( E_all + i * 9, mu, lambda, dW2_FD ); + second_deriv_fwd_fwd( E_all + i * 9, mu, lambda, dW2_fwd_fwd ); for ( int j = 0; j < 9; ++j ) { - d2W_d2E[9 * i + j] = d2E[j]; + EXPECT_NEAR( dW2_FD[j], dW2_fwd_fwd[j], tol ); } } } -void second_deriv_fwd_FD( double* E, double mu, double lambda, double* d2W_d2E, double h = 1e-7 ) +TEST_F( EnzymeNeohookeanTest, 2nd_deriv_finite_difference_vs_fwd_rev ) { - double E_plus[9] = { 0.0 }; - double E_minus[9] = { 0.0 }; - double dW_plus[9] = { 0.0 }; - double dW_minus[9] = { 0.0 }; - for ( int i = 0; i < 9; ++i ) { - for ( int j = 0; j < 9; ++j ) { - E_plus[j] = E[j]; - E_minus[j] = E[j]; - } - E_plus[i] = E[i] + h; - E_minus[i] = E[i] - h; - stress( E_plus, mu, lambda, dW_plus ); - stress( E_minus, mu, lambda, dW_minus ); - for ( int j = 0; j < 9; ++j ) { - d2W_d2E[9 * i + j] = ( dW_plus[j] - dW_minus[j] ) / ( 2 * h ); + for ( int i = 0; i < 3; ++i ) { + double dW2_FD[81] = { 0.0 }; + double dW2_fwd_rev[81] = { 0.0 }; + second_deriv_fwd_FD( E_all + i * 9, mu, lambda, dW2_FD ); + second_deriv_fwd_rev( E_all + i * 9, mu, lambda, dW2_fwd_rev ); + for ( int j = 0; j < 81; ++j ) { + EXPECT_NEAR( dW2_FD[j], dW2_fwd_rev[j], tol ); } } } -void second_deriv_hand_fwd( double* E, double mu, double lambda, double* d2W_d2E ) +TEST_F( EnzymeNeohookeanTest, 2nd_deriv_finite_difference_vs_rev_rev ) { - double dW[9] = { 0.0 }; - double d2w[9] = { 0.0 }; - - for ( int i = 0; i < 9; ++i ) { - double d2E[9] = { 0.0 }; - d2E[i] = 1.0; - double d2mu = 0.0; - double d2lambda = 0.0; - __enzyme_fwddiff( (void*)hand_code_deriv, E, d2E, mu, d2mu, lambda, d2lambda, &dW, &d2w ); - for ( int j = 0; j < 9; ++j ) { - d2W_d2E[9 * i + j] = d2w[j]; + for ( int i = 0; i < 3; ++i ) { + double dW2_FD[81] = { 0.0 }; + double dW2_rev_rev[81] = { 0.0 }; + second_deriv_fwd_FD( E_all + i * 9, mu, lambda, dW2_FD ); + second_deriv_rev_rev( E_all + i * 9, mu, lambda, dW2_rev_rev ); + for ( int j = 0; j < 81; ++j ) { + EXPECT_NEAR( dW2_FD[j], dW2_rev_rev[j], tol ); } } } -void second_deriv_hand_FD( double* E, double mu, double lambda, double* d2W_d2E, double h = 1e-7 ) +TEST_F( EnzymeNeohookeanTest, 2nd_deriv_finite_difference_vs_rev_fwd ) { - double E_plus[9] = { 0.0 }; - double E_minus[9] = { 0.0 }; - double dw_minus[9] = { 0.0 }; - double dw_plus[9] = { 0.0 }; - for ( int i = 0; i < 9; ++i ) { - for ( int j = 0; j < 9; ++j ) { - E_plus[j] = E[j]; - E_minus[j] = E[j]; - } - E_plus[i] = E[i] + h; - E_minus[i] = E[i] - h; - hand_code_deriv( E_plus, mu, lambda, dw_plus ); - hand_code_deriv( E_minus, mu, lambda, dw_minus ); - for ( int j = 0; j < 9; ++j ) { - d2W_d2E[9 * i + j] = ( dw_plus[j] - dw_minus[j] ) / ( 2 * h ); + for ( int i = 0; i < 3; ++i ) { + double dW2_FD[81] = { 0.0 }; + double dW2_rev_fwd[81] = { 0.0 }; + second_deriv_fwd_FD( E_all + i * 9, mu, lambda, dW2_FD ); + second_deriv_rev_rev( E_all + i * 9, mu, lambda, dW2_rev_fwd ); + for ( int j = 0; j < 81; ++j ) { + EXPECT_NEAR( dW2_FD[j], dW2_rev_fwd[j], tol ); } } } -class EnzymeNeohookeanTest : public testing::Test { - protected: - double F_all[27] = { 1.0, 0.5, 0.0, 0.0, 1.2, 0.1, 0.0, 0.0, 1.0, 1.1, 0.7, 0.1, 0.2, 1.4, - 0.2, 0.02, 0.04, 1.05, 1.05, 0.2, 0.04, 0.0, 1.1, 0.14, 0.07, 0.0, 1.09 }; - - double mu = 1.6; - double lambda = 1.2; - - double E_all[27] = { 0.0 }; - - double tol = 1e-6; - - void SetUp() override - { - for ( int i = 0; i < 3; ++i ) { - calc_E_from_F( F_all + i * 9, E_all + i * 9 ); - } - } -}; - -TEST_F( EnzymeNeohookeanTest, finite_difference_vs_exact ) +TEST_F( EnzymeNeohookeanTest, 2nd_deriv_finite_difference_vs_hand_fwd ) { for ( int i = 0; i < 3; ++i ) { - double S_FD[9] = { 0.0 }; - double S_exact[9] = { 0.0 }; - stress_FD( E_all + i * 9, mu, lambda, S_FD ); - hand_code_deriv( E_all + i * 9, mu, lambda, S_exact ); - for ( int j = 0; j < 9; ++j ) { - EXPECT_NEAR( S_FD[j], S_exact[j], tol ); - } + double dW2_FD[81] = { 0.0 }; + double dW2_hand_fwd[81] = { 0.0 }; + second_deriv_fwd_FD( E_all + i * 9, mu, lambda, dW2_FD ); + second_deriv_hand_fwd( E_all + i * 9, mu, lambda, dW2_hand_fwd ); + for(int j = 0; j < 9; ++j) { + EXPECT_NEAR( dW2_FD[j], dW2_hand_fwd[j], tol ); + } } } -// TEST(StrainEnergyTest, StressFiniteDifferenceVsAutodiff) { -// double F[9] = {1.0, 0.5, 0.0, 0.0, 1.2, 0.1, 0.0, 0.0, 1.0}; -// double E[9] = {0.0}; -// double mu = 1.0; -// double lambda = 1.0; -// double dW_dF_fd[9]; -// double dW_dF_ad[9]; -// double J = calc_determinant(F); -// calc_green_lagrange(F, E); -// std::cout << "Calling stress_FD..." << std::endl; -// stress_FD(E, J, mu, lambda, dW_dF_fd); -// std::cout << "Finite difference stress computed:" << std::endl; -// for (int i = 0; i < 9; ++i) { -// std::cout << "dW_dF_fd[" << i << "] = " << dW_dF_fd[i] << std::endl; -// } - -// std::cout << "Calling stress..." << std::endl; -// stress(E, J, mu, lambda, dW_dF_ad); -// std::cout << "Autodiff stress computed:" << std::endl; -// for (int i = 0; i < 9; ++i) { -// std::cout << "dW_dF_ad[" << i << "] = " << dW_dF_ad[i] << std::endl; -// } - -// for (int i = 0; i < 9; ++i) { -// std::cout << "Comparing index " << i << ": FD = " << dW_dF_fd[i] -// << ", AD = " << dW_dF_ad[i] << std::endl; -// EXPECT_NEAR(dW_dF_fd[i], dW_dF_ad[i], 1e-6) << "Mismatch at index " << i; -// } -// } - -// TEST(StrainEnergyTest, StressFiniteDifferenceVsReverseAutodiff) { -// double E[9] = {1.0, 0.5, 0.0, 0.0, 1.2, 0.1, 0.0, 0.0, 1.0}; -// double mu = 1.0; -// double lambda = 1.0; -// double dW_dF_fd[9]; -// double dW_dF_rev[9]; -// std::cout << "E: "; -// for (int i = 0; i < 9; ++i) std::cout << E[i] << " "; -// std::cout << std::endl; -// std::cout << "Calling stress_FD..." << std::endl; -// stress_FD(E, mu, lambda, dW_dF_fd); -// std::cout << "Finite difference stress computed:" << std::endl; -// for (int i = 0; i < 9; ++i) { -// std::cout << "dW_dF_fd[" << i << "] = " << dW_dF_fd[i] << std::endl; -// } - -// std::cout << "Calling stress_reverse..." << std::endl; -// stress_reverse(E, mu, lambda, dW_dF_rev); -// std::cout << "Reverse-mode autodiff stress computed:" << std::endl; -// for (int i = 0; i < 9; ++i) { -// std::cout << "dW_dF_rev[" << i << "] = " << dW_dF_rev[i] << std::endl; -// } - -// for (int i = 0; i < 9; ++i) { -// std::cout << "Comparing index " << i << ": FD = " << dW_dF_fd[i] -// << ", Reverse AD = " << dW_dF_rev[i] << std::endl; -// EXPECT_NEAR(dW_dF_fd[i], dW_dF_rev[i], 1e-6) << "Mismatch at index " << i; -// } -// } - -// TEST(StrainEnergyTest, StressFiniteDifferenceVsReverseAutodiff) { -// double E[9] = {0.1, 0.05, 0.02, 0.05, 0.2, 0.01, 0.02, 0.01, 0.15}; -// double mu = 1.0; -// double lambda = 1.0; -// double dW_dF_fd[9]; -// double dW_dF_hand[9]; -// std::cout << "E: "; -// for (int i = 0; i < 9; ++i) std::cout << E[i] << " "; -// std::cout << std::endl; -// std::cout << "Calling stress_FD..." << std::endl; -// stress_FD(E, mu, lambda, dW_dF_fd); -// std::cout << "Finite difference stress computed:" << std::endl; -// for (int i = 0; i < 9; ++i) { -// std::cout << "dW_dF_fd[" << i << "] = " << dW_dF_fd[i] << std::endl; -// } - -// std::cout << "Calling hand_code_deriv..." << std::endl; -// hand_code_deriv(E, mu, lambda, dW_dF_hand); -// std::cout << "Hand coded stress computed:" << std::endl; -// for (int i = 0; i < 9; ++i) { -// std::cout << "dW_dF_hand[" << i << "] = " << dW_dF_hand[i] << std::endl; -// } - -// for (int i = 0; i < 9; ++i) { -// std::cout << "Comparing index " << i << ": FD = " << dW_dF_fd[i] -// << ", Hand derivative = " << dW_dF_hand[i] << std::endl; -// EXPECT_NEAR(dW_dF_fd[i], dW_dF_hand[i], 1e-6) << "Mismatch at index " << i; -// } -// } - -// TEST(StrainEnergyTest, StressFiniteDifferenceVsReverseAutodiff) { -// double E[9] = {0.1, 0.05, 0.02, 0.05, 0.2, 0.01, 0.02, 0.01, 0.15}; -// double mu = 1.0; -// double lambda = 1.0; -// double dW_dF_fwd_fd[81]; -// double dW_dF_fwd_fwd[81]; -// std::cout << "E: "; -// for (int i = 0; i < 81; ++i) std::cout << E[i] << " "; -// std::cout << std::endl; -// std::cout << "Calling stress_FD..." << std::endl; -// second_deriv_fwd_FD(E, mu, lambda, dW_dF_fwd_fd); -// std::cout << "Finite difference stress computed:" << std::endl; -// for (int i = 0; i < 81; ++i) { -// std::cout << "dW_dF_fd[" << i << "] = " << dW_dF_fwd_fd[i] << std::endl; -// } - -// std::cout << "Calling hand_code_deriv..." << std::endl; -// second_deriv_fwd_fwd(E, mu, lambda, dW_dF_fwd_fwd); -// std::cout << "Hand coded stress computed:" << std::endl; -// for (int i = 0; i < 81; ++i) { -// std::cout << "dW_dF_hand[" << i << "] = " << dW_dF_fwd_fwd[i] << std::endl; -// } - -// for (int i = 0; i < 81; ++i) { -// std::cout << "Comparing index " << i << ": FD = " << dW_dF_fwd_fd[i] -// << ", Hand derivative = " << dW_dF_fwd_fwd[i] << std::endl; -// EXPECT_NEAR(dW_dF_fwd_fd[i], dW_dF_fwd_fwd[i], 1e-6) << "Mismatch at index " << i; -// } -// } - -// TEST(StrainEnergyTest, StressFiniteDifferenceVsReverseAutodiff) { -// double E[9] = {0.1, 0.05, 0.02, 0.05, 0.2, 0.01, 0.02, 0.01, 0.15}; -// double mu = 1.0; -// double lambda = 1.0; -// double dW_dF_fwd_fd[81]; -// double dW_dF_fwd_rev[81]; - -// std::cout << "E: "; -// for (int i = 0; i < 9; ++i) std::cout << E[i] << " "; -// std::cout << std::endl; - -// std::cout << "Calling stress_FD..." << std::endl; -// second_deriv_fwd_FD(E, mu, lambda, dW_dF_fwd_fd); -// std::cout << "Finite difference stress computed:" << std::endl; -// for (int i = 0; i < 81; ++i) { -// std::cout << "dW_dF_fd[" << i << "] = " << dW_dF_fwd_fd[i] << std::endl; -// } - -// std::cout << "Calling hand_code_deriv (fwd_rev)..." << std::endl; -// second_deriv_fwd_rev(E, mu, lambda, dW_dF_fwd_rev); -// std::cout << "Hand coded stress computed (fwd_rev):" << std::endl; -// for (int i = 0; i < 81; ++i) { -// std::cout << "dW_dF_hand[" << i << "] = " << dW_dF_fwd_rev[i] << std::endl; -// } - -// for (int i = 0; i < 81; ++i) { -// std::cout << "Comparing index " << i << ": FD = " << dW_dF_fwd_fd[i] -// << ", Hand derivative = " << dW_dF_fwd_rev[i] << std::endl; -// EXPECT_NEAR(dW_dF_fwd_fd[i], dW_dF_fwd_rev[i], 1e-6) << "Mismatch at index " << i; -// } -// } - -// TEST(StrainEnergyTest, StressFiniteDifferenceVsReverseAutodiff) { -// double E[9] = {0.1, 0.05, 0.02, 0.05, 0.2, 0.01, 0.02, 0.01, 0.15}; -// double mu = 1.0; -// double lambda = 1.0; -// double dW_dF_fwd_fd[81]; -// double dW_dF_rev_rev[81]; - -// std::cout << "E: "; -// for (int i = 0; i < 9; ++i) std::cout << E[i] << " "; -// std::cout << std::endl; - -// std::cout << "Calling stress_FD..." << std::endl; -// second_deriv_fwd_FD(E, mu, lambda, dW_dF_fwd_fd); -// std::cout << "Finite difference stress computed:" << std::endl; -// for (int i = 0; i < 81; ++i) { -// std::cout << "dW_dF_fd[" << i << "] = " << dW_dF_fwd_fd[i] << std::endl; -// } - -// std::cout << "Calling hand_code_deriv (rev_rev)..." << std::endl; -// second_deriv_rev_rev(E, mu, lambda, dW_dF_rev_rev); -// std::cout << "Hand coded stress computed (rev_rev):" << std::endl; -// for (int i = 0; i < 81; ++i) { -// std::cout << "dW_dF_hand[" << i << "] = " << dW_dF_rev_rev[i] << std::endl; -// } - -// for (int i = 0; i < 81; ++i) { -// std::cout << "Comparing index " << i << ": FD = " << dW_dF_fwd_fd[i] -// << ", Hand derivative = " << dW_dF_rev_rev[i] << std::endl; -// EXPECT_NEAR(dW_dF_fwd_fd[i], dW_dF_rev_rev[i], 1e-6) << "Mismatch at index " << i; -// } -// } - -// TEST(StrainEnergyTest, StressFiniteDifferenceVsFwdRevAutodiff) { -// double E[9] = {0.1, 0.05, 0.02, 0.05, 0.2, 0.01, 0.02, 0.01, 0.15}; -// double mu = 1.0; -// double lambda = 1.0; -// double dW_dF_fwd_fd[81]; -// double dW_dF_fwd_rev[81]; - -// std::cout << "E: "; -// for (int i = 0; i < 9; ++i) std::cout << E[i] << " "; -// std::cout << std::endl; - -// std::cout << "Calling stress_FD..." << std::endl; -// second_deriv_fwd_FD(E, mu, lambda, dW_dF_fwd_fd); -// std::cout << "Finite difference stress computed:" << std::endl; -// for (int i = 0; i < 81; ++i) { -// std::cout << "dW_dF_fd[" << i << "] = " << dW_dF_fwd_fd[i] << std::endl; -// } - -// std::cout << "Calling hand_code_deriv (fwd_rev)..." << std::endl; -// second_deriv_fwd_rev(E, mu, lambda, dW_dF_fwd_rev); -// std::cout << "Hand coded stress computed (fwd_rev):" << std::endl; -// for (int i = 0; i < 81; ++i) { -// std::cout << "dW_dF_hand[" << i << "] = " << dW_dF_fwd_rev[i] << std::endl; -// } - -// for (int i = 0; i < 81; ++i) { -// std::cout << "Comparing index " << i << ": FD = " << dW_dF_fwd_fd[i] -// << ", Hand derivative = " << dW_dF_fwd_rev[i] << std::endl; -// EXPECT_NEAR(dW_dF_fwd_fd[i], dW_dF_fwd_rev[i], 1e-6) << "Mismatch at index " << i; -// } -// } - -// TEST(StrainEnergyTest, StressFiniteDifferenceVsFwdRevAutodiff) { -// double E[9] = {0.12, -0.03, 0.01, -0.03, 0.08, 0.02, 0.01, 0.02, 0.11}; -// double mu = 1.0; -// double lambda = 1.0; -// double dW_dF_fwd_fd[81]; -// double dW_dF_fwd_rev[81]; - -// std::cout << "E: "; -// for (int i = 0; i < 9; ++i) std::cout << E[i] << " "; -// std::cout << std::endl; - -// std::cout << "Calling stress_FD..." << std::endl; -// second_deriv_fwd_FD(E, mu, lambda, dW_dF_fwd_fd); -// std::cout << "Finite difference stress computed:" << std::endl; -// for (int i = 0; i < 81; ++i) { -// std::cout << "dW_dF_fd[" << i << "] = " << dW_dF_fwd_fd[i] << std::endl; -// } - -// std::cout << "Calling hand_code_deriv (fwd_rev)..." << std::endl; -// second_deriv_hand_fwd(E, mu, lambda, dW_dF_fwd_rev); -// std::cout << "Hand coded stress computed (fwd_rev):" << std::endl; -// for (int i = 0; i < 81; ++i) { -// std::cout << "dW_dF_hand[" << i << "] = " << dW_dF_fwd_rev[i] << std::endl; -// } - -// for (int i = 0; i < 81; ++i) { -// std::cout << "Comparing index " << i << ": FD = " << dW_dF_fwd_fd[i] -// << ", Hand derivative = " << dW_dF_fwd_rev[i] << std::endl; -// EXPECT_NEAR(dW_dF_fwd_fd[i], dW_dF_fwd_rev[i], 1e-6) << "Mismatch at index " << i; -// } -// } - -// TEST(StrainEnergyTest, StressFiniteDifferenceVsFwdRevAutodiff) { -// double F[9] = {1.01, -0.03, 0.01, -0.03, 1.05, 0.02, 0.01, 0.02, 1.0}; -// double E[9] = {0}; -// double mu = 1.0; -// double lambda = 1.0; -// double* S = nullptr; -// calc_E_from_F(F, E); -// hand_code_deriv(E, mu, lambda, S); -// std::cout << "{ "; -// for (int i = 0; i < 9; ++i) { -// std::cout << ", " << S[i]; -// } -// std::cout << " }" << std::endl; -// std::cout << "E: { "; -// for (int i = 0; i < 9; ++i) { -// std::cout << ", " << E[i]; -// } -// std::cout << "}" << std::endl; - -// double dW_dF_hand_fd[81]; -// double dW_dF_hand_fwd[81]; - -// std::cout << "E: "; -// for (int i = 0; i < 9; ++i) std::cout << E[i] << " "; -// std::cout << std::endl; - -// std::cout << "Calling stress_FD..." << std::endl; -// second_deriv_hand_FD(E, mu, lambda, dW_dF_hand_fd); -// std::cout << "Finite difference stress computed:" << std::endl; -// for (int i = 0; i < 81; ++i) { -// std::cout << "dW_dF_hand_fd[" << i << "] = " << dW_dF_hand_fd[i] << std::endl; -// } - -// std::cout << "Calling hand_code_deriv (fwd_rev)..." << std::endl; -// second_deriv_hand_fwd(E, mu, lambda, dW_dF_hand_fwd); -// std::cout << "Hand coded stress computed (fwd_rev):" << std::endl; -// for (int i = 0; i < 81; ++i) { -// std::cout << "dW_dF_hand_fwd[" << i << "] = " << dW_dF_hand_fwd[i] << std::endl; -// } - -// for (int i = 0; i < 81; ++i) { -// std::cout << "Comparing index " << i << ": FD = " << dW_dF_hand_fd[i] -// << ", Hand derivative = " << dW_dF_hand_fwd[i] << std::endl; -// EXPECT_NEAR(dW_dF_hand_fd[i], dW_dF_hand_fd[i], 1e-6) << "Mismatch at index " << i; -// } -// } \ No newline at end of file diff --git a/src/tests/neohookean_common.hpp b/src/tests/neohookean_common.hpp new file mode 100644 index 00000000..32df6730 --- /dev/null +++ b/src/tests/neohookean_common.hpp @@ -0,0 +1,293 @@ +#include +#include +#include +#include +#include +#include +#include "tribol/common/Enzyme.hpp" + +void multiply3x3( const double F[9], const double F_T[9], double C[9] ) +{ + for ( int i = 0; i < 3; ++i ) { // row of A + for ( int j = 0; j < 3; ++j ) { // column of B + C[i * 3 + j] = 0.0; + for ( int k = 0; k < 3; ++k ) { + C[i * 3 + j] += F_T[i * 3 + k] * F[k * 3 + j]; + } + } + } +} + +void calc_E_from_F( const double F[9], double E[9] ) +{ + double I[9] = { 1.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 1.0 }; + double F_T[9] = { 0.0 }; + F_T[0] = F[0]; + F_T[1] = F[3]; + F_T[2] = F[6]; + F_T[3] = F[1]; + F_T[4] = F[4]; + F_T[5] = F[7]; + F_T[6] = F[2]; + F_T[7] = F[5]; + F_T[8] = F[8]; + double C[9] = { 0.0 }; + multiply3x3( F, F_T, C ); + for ( int i = 0; i < 9; ++i ) { + E[i] = 0.5 * ( C[i] - I[i] ); + } +} + +// Calculates right cauchy stress tensor +void calc_cauchy_stress_tensor( const double E[9], double C[9] ) +{ + double I[9] = { 1, 0, 0, 0, 1, 0, 0, 0, 1 }; + for ( int i = 0; i < 9; ++i ) { + C[i] = 2 * E[i]; + C[i] += I[i]; + } +} + +// Calculates Trace +double calc_trace( const double C[9] ) +{ + double Tr_C = C[0] + C[4] + C[8]; + return Tr_C; +} + +double calc_determinant( const double C[9] ) +{ + double J = C[0] * ( C[4] * C[8] - C[5] * C[7] ) - C[1] * ( C[3] * C[8] - C[5] * C[6] ) + + C[2] * ( C[3] * C[7] - C[4] * C[6] ); + J = std::sqrt( J ); + return J; +} + +bool invert3x3( const double F[9], double Finv[9] ) +{ + // Compute the determinant + double det = F[0] * ( F[4] * F[8] - F[5] * F[7] ) - F[1] * ( F[3] * F[8] - F[5] * F[6] ) + + F[2] * ( F[3] * F[7] - F[4] * F[6] ); + + if ( std::abs( det ) < std::numeric_limits::epsilon() ) return false; // Singular matrix + + double invDet = 1.0 / det; + + // Compute the inverse using the formula for the inverse of a 3x3 matrix + Finv[0] = ( F[4] * F[8] - F[5] * F[7] ) * invDet; + Finv[1] = -( F[1] * F[8] - F[2] * F[7] ) * invDet; + Finv[2] = ( F[1] * F[5] - F[2] * F[4] ) * invDet; + + Finv[3] = -( F[3] * F[8] - F[5] * F[6] ) * invDet; + Finv[4] = ( F[0] * F[8] - F[2] * F[6] ) * invDet; + Finv[5] = -( F[0] * F[5] - F[2] * F[3] ) * invDet; + + Finv[6] = ( F[3] * F[7] - F[4] * F[6] ) * invDet; + Finv[7] = -( F[0] * F[7] - F[1] * F[6] ) * invDet; + Finv[8] = ( F[0] * F[4] - F[1] * F[3] ) * invDet; + + return true; +} + +// build strain energy equation +void strain_energy( double* E, double mu, double lambda, double* W ) +{ + double C[9] = { 0.0 }; + calc_cauchy_stress_tensor( E, C ); + double Tr_C = calc_trace( C ); + double J = calc_determinant( C ); + *W = mu / 2.0 * ( Tr_C - 3.0 ) - mu * log( J ) + lambda / 2.0 * pow( ( log( J ) ), 2.0 ); +} + +// calc stress using enzyme fwddiff +void stress( double* E, double mu, double lambda, double* dW_dE ) +{ + double W = 0.0; + for ( int i = 0; i < 9; ++i ) { + double dE[9] = { 0.0 }; + dE[i] = 1.0; + double dmu = 0.0; + double dlambda = 0.0; + double dw = 0.0; + __enzyme_fwddiff( (void*)strain_energy, E, dE, mu, dmu, lambda, dlambda, &W, &dw ); + dW_dE[i] = dw; + } +} + +// calc stress using enzyme autodiff +void stress_reverse( double* E, double mu, double lambda, double* dW_dE ) +{ + double dE[9] = { 0.0 }; + double W = 0.0; + double dW = 1.0; + __enzyme_autodiff( strain_energy, enzyme_dup, E, dE, enzyme_const, mu, enzyme_const, lambda, enzyme_dup, &W, + &dW ); + + for ( int i = 0; i < 9; ++i ) { + dW_dE[i] = dE[i]; + } +} + +// calc stress using finite Difference +void stress_FD( double* E, double mu, double lambda, double* dW_dE, double h = 1e-7 ) +{ + double E_plus[9] = { 0.0 }; + double E_minus[9] = { 0.0 }; + double W_plus; + double W_minus; + for ( int i = 0; i < 9; ++i ) { + for ( int j = 0; j < 9; ++j ) { + E_plus[j] = E[j]; + E_minus[j] = E[j]; + } + E_plus[i] = E[i] + h; + E_minus[i] = E[i] - h; + strain_energy( E_plus, mu, lambda, &W_plus ); + strain_energy( E_minus, mu, lambda, &W_minus ); + dW_dE[i] = ( W_plus - W_minus ) / ( 2 * h ); + } +} + +void hand_code_deriv( double* E, double mu, double lambda, double* S ) +{ + double C[9] = { 0.0 }; + double Cinv[9]; + double I[9] = { 1.0, 0.0, 0.0, 0, 1.0, 0.0, 0.0, 0.0, 1.0 }; + calc_cauchy_stress_tensor( E, C ); + double J = calc_determinant( C ); + invert3x3( C, Cinv ); + double first_term[9]; + for ( int i = 0; i < 9; ++i ) { + first_term[i] = lambda * std::log( J ) * Cinv[i]; + } + double second_term[9]; + for ( int i = 0; i < 9; ++i ) { + second_term[i] = I[i] - Cinv[i]; + second_term[i] *= mu; + } + for ( int i = 0; i < 9; ++i ) { + S[i] = first_term[i] + second_term[i]; + } +} + +void second_deriv_fwd_fwd( double* E, double mu, double lambda, double* d2W_d2E ) +{ + double dW[9] = { 0.0 }; + double d2w[9] = { 0.0 }; + + for ( int i = 0; i < 9; ++i ) { + double d2E[9] = { 0.0 }; + d2E[i] = 1.0; + double d2mu = 0.0; + double d2lambda = 0.0; + __enzyme_fwddiff( (void*)stress, E, d2E, mu, d2mu, lambda, d2lambda, &dW, &d2w ); + for ( int j = 0; j < 9; ++j ) { + d2W_d2E[9 * i + j] = d2w[j]; + } + } +} + +void second_deriv_rev_fwd( double* E, double mu, double lambda, double* d2W_d2E ) +{ + double dW[9] = { 0.0 }; + double d2w[9] = { 0.0 }; + + for ( int i = 0; i < 9; ++i ) { + double d2E[9] = { 0.0 }; + d2E[i] = 1.0; + double d2mu = 0.0; + double d2lambda = 0.0; + __enzyme_fwddiff( (void*)stress_reverse, E, d2E, mu, d2mu, lambda, d2lambda, &dW, &d2w ); + for ( int j = 0; j < 9; ++j ) { + d2W_d2E[9 * i + j] = d2w[j]; + } + } +} + +void second_deriv_rev_rev( double* E, double mu, double lambda, double* d2W_d2E ) +{ + for ( int i = 0; i < 9; ++i ) { + double d2E[81] = { 0.0 }; + double W[9] = { 0.0 }; + double d2W[9] = { 0.0 }; + d2W[i] = 1.0; + __enzyme_autodiff( stress_reverse, enzyme_dup, E, d2E, enzyme_const, mu, enzyme_const, lambda, enzyme_dup, &W, + &d2W ); + for ( int j = 0; j < 9; ++j ) { + d2W_d2E[9 * j + i] = d2E[j]; + } + } +} + +void second_deriv_fwd_rev( double* E, double mu, double lambda, double* d2W_d2E ) +{ + for ( int i = 0; i < 9; ++i ) { + double d2E[81] = { 0.0 }; + double W[9] = { 0.0 }; + double d2W[9] = { 0.0 }; + d2W[i] = 1.0; + __enzyme_autodiff( stress, enzyme_dup, E, d2E, enzyme_const, mu, enzyme_const, lambda, enzyme_dup, &W, &d2W ); + for ( int j = 0; j < 9; ++j ) { + d2W_d2E[9 * i + j] = d2E[j]; + } + } +} + +void second_deriv_fwd_FD( double* E, double mu, double lambda, double* d2W_d2E, double h = 1e-7 ) +{ + double E_plus[9] = { 0.0 }; + double E_minus[9] = { 0.0 }; + double dW_plus[9] = { 0.0 }; + double dW_minus[9] = { 0.0 }; + for ( int i = 0; i < 9; ++i ) { + for ( int j = 0; j < 9; ++j ) { + E_plus[j] = E[j]; + E_minus[j] = E[j]; + } + E_plus[i] = E[i] + h; + E_minus[i] = E[i] - h; + stress( E_plus, mu, lambda, dW_plus ); + stress( E_minus, mu, lambda, dW_minus ); + for ( int j = 0; j < 9; ++j ) { + d2W_d2E[9 * i + j] = ( dW_plus[j] - dW_minus[j] ) / ( 2 * h ); + } + } +} + +void second_deriv_hand_fwd( double* E, double mu, double lambda, double* d2W_d2E ) +{ + double dW[9] = { 0.0 }; + double d2w[9] = { 0.0 }; + + for ( int i = 0; i < 9; ++i ) { + double d2E[9] = { 0.0 }; + d2E[i] = 1.0; + double d2mu = 0.0; + double d2lambda = 0.0; + __enzyme_fwddiff( (void*)hand_code_deriv, E, d2E, mu, d2mu, lambda, d2lambda, &dW, &d2w ); + for ( int j = 0; j < 9; ++j ) { + d2W_d2E[9 * i + j] = d2w[j]; + } + } +} + +void second_deriv_hand_FD( double* E, double mu, double lambda, double* d2W_d2E, double h = 1e-7 ) +{ + double E_plus[9] = { 0.0 }; + double E_minus[9] = { 0.0 }; + double dw_minus[9] = { 0.0 }; + double dw_plus[9] = { 0.0 }; + for ( int i = 0; i < 9; ++i ) { + for ( int j = 0; j < 9; ++j ) { + E_plus[j] = E[j]; + E_minus[j] = E[j]; + } + E_plus[i] = E[i] + h; + E_minus[i] = E[i] - h; + hand_code_deriv( E_plus, mu, lambda, dw_plus ); + hand_code_deriv( E_minus, mu, lambda, dw_minus ); + for ( int j = 0; j < 9; ++j ) { + d2W_d2E[9 * i + j] = ( dw_plus[j] - dw_minus[j] ) / ( 2 * h ); + } + } +} \ No newline at end of file From ad1561b8ab38cc4130ce36a2ebb6ce8dc554e4fc Mon Sep 17 00:00:00 2001 From: Ryan Lutz Date: Tue, 17 Jun 2025 16:37:20 -0700 Subject: [PATCH 3/4] updated hpp and finalized neohookean_timing --- src/examples/neohookean_timing.cpp | 740 +---------------------------- src/tests/enzyme_neohookean.cpp | 8 +- src/tests/neohookean_common.hpp | 110 +++++ 3 files changed, 137 insertions(+), 721 deletions(-) diff --git a/src/examples/neohookean_timing.cpp b/src/examples/neohookean_timing.cpp index a694de55..8aee35f7 100644 --- a/src/examples/neohookean_timing.cpp +++ b/src/examples/neohookean_timing.cpp @@ -1,7 +1,7 @@ #include -#include +#include "../tests/neohookean_common.hpp" +#include #include -#include "gtest/gtest.h" #include #include "axom/core.hpp" #include "axom/slic/interface/slic.hpp" @@ -9,719 +9,27 @@ #include #include "tribol/common/Enzyme.hpp" - -void multiply3x3(const double F[9], const double F_T[9], double C[9]) { - for (int i = 0; i < 3; ++i) { // row of A - for (int j = 0; j < 3; ++j) { // column of B - C[i * 3 + j] = 0.0; - for (int k = 0; k < 3; ++k) { - C[i * 3 + j] += F_T[i * 3 + k] * F[k * 3 + j]; - } - } - } -} - -void calc_E_from_F(const double F[9], double E[9]) { - double I[9] = {1.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 1.0}; - double F_T[9] = {0.0}; - F_T[0] = F[0]; - F_T[1] = F[3]; - F_T[2] = F[6]; - F_T[3] = F[1]; - F_T[4] = F[4]; - F_T[5] = F[7]; - F_T[6] = F[2]; - F_T[7] = F[5]; - F_T[8] = F[8]; - double C[9] = {0.0}; - multiply3x3(F, F_T, C); - for(int i = 0; i < 9; ++i) { - E[i] = 0.5 * (C[i] - I[i]); - } -} - -//Calculates right cauchy stress tensor -void calc_cauchy_stress_tensor(const double E[9], double C[9]){ - double I[9] = {1, 0, 0, 0, 1, 0, 0, 0, 1}; - for(int i = 0; i < 9; ++i) { - C[i] = 2 * E[i]; - C[i] += I[i]; - } -} - - -//Calculates Trace -double calc_trace(const double C[9]) { - double Tr_C = C[0] + C[4] + C[8]; - return Tr_C; -} - - - -double calc_determinant(const double C[9]) { - double J = C[0] * (C[4] * C[8] - C[5] * C[7]) - C[1] * (C[3] * C[8] - C[5] * C[6]) + C[2] * (C[3] * C[7] - C[4] * C[6]); - J = std::sqrt(J); - return J; -} - - -bool invert3x3(const double F[9], double Finv[9]) +int main() { - // Compute the determinant - double det = - F[0]*(F[4]*F[8] - F[5]*F[7]) - - F[1]*(F[3]*F[8] - F[5]*F[6]) + - F[2]*(F[3]*F[7] - F[4]*F[6]); - - if (std::abs(det) < std::numeric_limits::epsilon()) - return false; // Singular matrix - - double invDet = 1.0 / det; - - // Compute the inverse using the formula for the inverse of a 3x3 matrix - Finv[0] = (F[4]*F[8] - F[5]*F[7]) * invDet; - Finv[1] = -(F[1]*F[8] - F[2]*F[7]) * invDet; - Finv[2] = (F[1]*F[5] - F[2]*F[4]) * invDet; - - Finv[3] = -(F[3]*F[8] - F[5]*F[6]) * invDet; - Finv[4] = (F[0]*F[8] - F[2]*F[6]) * invDet; - Finv[5] = -(F[0]*F[5] - F[2]*F[3]) * invDet; - - Finv[6] = (F[3]*F[7] - F[4]*F[6]) * invDet; - Finv[7] = -(F[0]*F[7] - F[1]*F[6]) * invDet; - Finv[8] = (F[0]*F[4] - F[1]*F[3]) * invDet; - - return true; -} - - -//build strain energy equation -void strain_energy(double* E, double mu, double lambda, double* W) { - double C[9] = {0.0}; - calc_cauchy_stress_tensor(E, C); - double Tr_C = calc_trace(C); - double J = calc_determinant(C); - *W = mu/2.0 * (Tr_C - 3.0) - mu * log(J) + lambda/2.0 * pow((log(J)), 2.0); -} - -//calc stress using enzyme fwddiff -void stress(double* E, double mu, double lambda, double* dW_dE) { - double W = 0.0; - for(int i = 0; i < 9; ++i) { - double dE[9] = {0.0}; - dE[i] = 1.0; - double dmu = 0.0; - double dlambda = 0.0; - double dw = 0.0; - __enzyme_fwddiff( (void*) strain_energy, E, dE, mu, dmu, lambda, dlambda, &W, &dw); - dW_dE[i] = dw; - } - -} - -//calc stress using enzyme autodiff -void stress_reverse(double* E, double mu, double lambda, double* dW_dE) { - double dE[9] = {0.0}; - double W = 0.0; - double dW = 1.0; - __enzyme_autodiff( strain_energy, enzyme_dup, E, dE, enzyme_const, mu, enzyme_const, lambda, enzyme_dup, &W, &dW); - - for (int i = 0; i < 9; ++i) { - dW_dE[i] = dE[i]; - } -} - - -//calc stress using finite Difference -void stress_FD(double* E, double mu, double lambda, double* dW_dE, double h = 1e-7) { - double E_plus[9] = {0.0}; - double E_minus[9] = {0.0}; - double W_plus; - double W_minus; - for(int i = 0; i < 9; ++i) { - for(int j = 0; j < 9; ++j) { - E_plus[j] = E[j]; - E_minus[j] = E[j]; - } - E_plus[i] = E[i] + h; - E_minus[i] = E[i] - h; - strain_energy(E_plus, mu, lambda, &W_plus); - strain_energy(E_minus, mu, lambda, &W_minus); - dW_dE[i] = (W_plus - W_minus) / (2 * h); - - } - -} - - -void hand_code_deriv(double* E, double mu, double lambda, double* S) { - double C[9] = {0.0}; - double Cinv[9]; - double I[9] = {1.0, 0.0, 0.0, 0, 1.0, 0.0, 0.0, 0.0, 1.0}; - calc_cauchy_stress_tensor(E, C); - double J = calc_determinant(C); - invert3x3(C, Cinv); - double first_term[9]; - for(int i = 0; i < 9; ++i) { - first_term[i] = lambda * std::log(J) * Cinv[i]; - } - double second_term[9]; - for(int i = 0; i < 9; ++i) { - second_term[i] = I[i] - Cinv[i]; - second_term[i] *= mu; - } - for(int i = 0; i < 9; ++i) { - S[i] = first_term[i] + second_term[i]; - } - -} - -void second_deriv_fwd_fwd(double* E, double mu, double lambda, double* d2W_d2E) { - double dW[9] = {0.0}; - double d2w[9] = {0.0}; - - for(int i = 0; i < 9; ++i) { - double d2E[9] = {0.0}; - d2E[i] = 1.0; - double d2mu = 0.0; - double d2lambda = 0.0; - __enzyme_fwddiff ( (void*) stress, E, d2E, mu, d2mu, lambda, d2lambda, &dW, &d2w ); - for(int j = 0; j < 9; ++j) { - d2W_d2E[9 * i + j] = d2w[j]; - } - } -} - -void second_deriv_rev_fwd(double* E, double mu, double lambda, double* d2W_d2E) { - double dW[9] = {0.0}; - double d2w[9] = {0.0}; - - for(int i = 0; i < 9; ++i) { - double d2E[9] = {0.0}; - d2E[i] = 1.0; - double d2mu = 0.0; - double d2lambda = 0.0; - __enzyme_fwddiff ( (void*) stress_reverse, E, d2E, mu, d2mu, lambda, d2lambda, &dW, &d2w ); - for(int j = 0; j < 9; ++j) { - d2W_d2E[9 * i + j] = d2w[j]; - } - } -} - -void second_deriv_rev_rev(double* E, double mu, double lambda, double* d2W_d2E) { - for (int i = 0; i < 9; ++i) { - double d2E[81] = {0.0}; - double W[9] = {0.0}; - double d2W[9] = {0.0}; - d2W[i] = 1.0; - __enzyme_autodiff( stress_reverse, enzyme_dup, E, d2E, enzyme_const, mu, enzyme_const, lambda, enzyme_dup, &W, &d2W); - for(int j = 0; j < 9; ++j) { - d2W_d2E[9 * j + i] = d2E[j]; - } - } -} - -void second_deriv_fwd_rev(double* E, double mu, double lambda, double* d2W_d2E) { - for (int i = 0; i < 9; ++i) { - double d2E[81] = {0.0}; - double W[9] = {0.0}; - double d2W[9] = {0.0}; - d2W[i] = 1.0; - __enzyme_autodiff( stress, enzyme_dup, E, d2E, enzyme_const, mu, enzyme_const, lambda, enzyme_dup, &W, &d2W); - for(int j = 0; j < 9; ++j) { - d2W_d2E[9 * i + j] = d2E[j]; - } - } -} - - -void second_deriv_fwd_FD(double* E, double mu, double lambda, double* d2W_d2E, double h = 1e-7){ - double E_plus[9] = {0.0}; - double E_minus[9] = {0.0}; - double dW_plus[9] = {0.0}; - double dW_minus[9] = {0.0}; - for(int i = 0; i < 9; ++i) { - for(int j = 0; j < 9; ++j) { - E_plus[j] = E[j]; - E_minus[j] = E[j]; - } - E_plus[i] = E[i] + h; - E_minus[i] = E[i] - h; - stress(E_plus, mu, lambda, dW_plus); - stress(E_minus, mu, lambda, dW_minus); - for(int j = 0; j < 9; ++j) { - d2W_d2E[9 * i + j] = (dW_plus[j] - dW_minus[j]) / (2 * h); - } - } -} - -void second_deriv_hand_fwd(double *E, double mu, double lambda, double* d2W_d2E) { - double dW[9] = {0.0}; - double d2w[9] = {0.0}; - - - for(int i = 0; i < 9; ++i) { - double d2E[9] = {0.0}; - d2E[i] = 1.0; - double d2mu = 0.0; - double d2lambda = 0.0; - __enzyme_fwddiff ( (void*) hand_code_deriv, E, d2E, mu, d2mu, lambda, d2lambda, &dW, &d2w ); - for(int j = 0; j < 9; ++j) { - d2W_d2E[9 * i + j] = d2w[j]; - } - } -} - -void second_deriv_hand_FD(double* E, double mu, double lambda, double* d2W_d2E, double h = 1e-7) { - double E_plus[9] = {0.0}; - double E_minus[9] = {0.0}; - double dw_minus[9] = {0.0}; - double dw_plus[9] = {0.0}; - for(int i = 0; i < 9; ++i){ - for (int j = 0; j < 9; ++j) { - E_plus[j] = E[j]; - E_minus[j] = E[j]; - } - E_plus[i] = E[i] + h; - E_minus[i] = E[i] - h; - hand_code_deriv(E_plus, mu, lambda, dw_plus); - hand_code_deriv(E_minus, mu, lambda, dw_minus); - for(int j = 0; j < 9; ++j) { - d2W_d2E[9 * i + j] = (dw_plus[j] - dw_minus[j]) / (2 * h); - } - } -} - - - - - -void run_fwd_mode(double* E, double mu, double lambda, double* dw_df, int N) { - axom::utilities::Timer timer{ false }; - stress(E, mu, lambda, dw_df); - double Dw[9] = {0.0}; - timer.start(); - for(int i = 0; i < N; ++i) { - stress(E, mu, lambda, dw_df); - for(int j = 0; j < 9; ++j) { - Dw[j] += dw_df[j]; - } - } - timer.stop(); - std::cout << axom::fmt::format( "Time to calc fwd_diff: {0:f}ms", timer.elapsedTimeInMilliSec() ) << std::endl; + double E[9] = { 0.1, 0.05, 0.02, 0.05, 0.2, 0.01, 0.02, 0.01, 0.15 }; + double mu = 1.0; + double lambda = 1.0; + int N = 0; + std::cout << "Enter cycles: "; + std::cin >> N; + double dw_dE[9] = { 0.0 }; + double d2W_d2E[81] = { 0.0 }; + double h = 1e-7; + + run_fwd_mode( E, mu, lambda, dw_dE, N ); + run_bkwd_mode( E, mu, lambda, dw_dE, N ); + run_hand_derivative( E, mu, lambda, dw_dE, N ); + run_fwd_fwd( E, mu, lambda, d2W_d2E, N ); + run_rev_fwd( E, mu, lambda, d2W_d2E, N ); + run_fwd_rev( E, mu, lambda, d2W_d2E, N ); + run_rev_rev( E, mu, lambda, d2W_d2E, N ); + run_fwd_FD( E, mu, lambda, d2W_d2E, N, h ); + run_hand_fwd( E, mu, lambda, d2W_d2E, N ); + + return 0; } - -void run_bkwd_mode(double* E, double mu, double lambda, double* dw_dE, int N) { - axom::utilities::Timer timer{ false }; - - stress_reverse(E, mu, lambda, dw_dE); - double Dw[9] = {0.0}; - - timer.start(); - for(int i = 0; i < N; ++i) { - stress_reverse(E, mu, lambda, dw_dE); - for(int j = 0; j < 9; ++j) { - Dw[j] += dw_dE[j]; - } - } - timer.stop(); - std::cout << axom::fmt::format( "Time to calc backward_diff: {0:f}ms", timer.elapsedTimeInMilliSec() ) << std::endl; -} - -void run_hand_derivative(double* E, double mu, double lambda, double* S, int N) { - axom::utilities::Timer timer{ false }; - hand_code_deriv(E, mu, lambda, S); - double Dw[9] = {0.0}; - timer.start(); - for(int i = 0; i < N; ++i) { - hand_code_deriv(E, mu, lambda, S); - for(int j = 0; j < 9; ++j) { - Dw[j] += S[j]; - } - } - timer.stop(); - std::cout << axom::fmt::format( "Time to calc hand_diff: {0:f}ms", timer.elapsedTimeInMilliSec() ) << std::endl; -} - -void run_fwd_fwd(double* E, double mu, double lambda, double* S, int N) { - axom::utilities::Timer timer{ false }; - timer.start(); - for(int i = 0; i < N; ++i) { - second_deriv_fwd_fwd(E, mu, lambda, S); - } - timer.stop(); - std::cout << axom::fmt::format( "Time to calc fwd_fwd_diff: {0:f}ms", timer.elapsedTimeInMilliSec() ) << std::endl; -} - -void run_rev_fwd(double* E, double mu, double lambda, double* S, int N) { - axom::utilities::Timer timer{ false }; - timer.start(); - for(int i = 0; i < N; ++i) { - second_deriv_rev_fwd(E, mu, lambda, S); - } - timer.stop(); - std::cout << axom::fmt::format( "Time to calc rev_fwd_diff: {0:f}ms", timer.elapsedTimeInMilliSec() ) << std::endl; -} - -void run_fwd_rev(double* E, double mu, double lambda, double* S, int N) { - axom::utilities::Timer timer{ false }; - timer.start(); - for(int i = 0; i < N; ++i) { - second_deriv_fwd_rev(E, mu, lambda, S); - } - timer.stop(); - std::cout << axom::fmt::format( "Time to calc fwd_rev_diff: {0:f}ms", timer.elapsedTimeInMilliSec() ) << std::endl; -} - -void run_rev_rev(double* E, double mu, double lambda, double* S, int N) { - axom::utilities::Timer timer{ false }; - timer.start(); - for(int i = 0; i < N; ++i) { - second_deriv_rev_rev(E, mu, lambda, S); - } - timer.stop(); - std::cout << axom::fmt::format( "Time to calc rev_rev_diff: {0:f}ms", timer.elapsedTimeInMilliSec() ) << std::endl; -} - -void run_fwd_FD(double* E, double mu, double lambda, double* S, int N, double h) { - axom::utilities::Timer timer{ false }; - timer.start(); - for(int i = 0; i < N; ++i) { - second_deriv_fwd_FD(E, mu, lambda, S, h); - } - timer.stop(); - std::cout << axom::fmt::format( "Time to calc fwd_FD_diff: {0:f}ms", timer.elapsedTimeInMilliSec() ) << std::endl; -} - -void run_hand_fwd(double* E, double mu, double lambda, double* S, int N) { - axom::utilities::Timer timer{ false }; - timer.start(); - for(int i = 0; i < N; ++i) { - second_deriv_hand_fwd(E, mu, lambda, S); - } - timer.stop(); - std::cout << axom::fmt::format( "Time to calc hand_fwd_diff: {0:f}ms", timer.elapsedTimeInMilliSec() ) << std::endl; -} - - -int main() { - double E[9] = {0.1, 0.05, 0.02, 0.05, 0.2, 0.01, 0.02, 0.01, 0.15}; - double mu = 1.0; - double lambda = 1.0; - int N = 0; - std::cout << "Enter Number: "; - std::cin >> N; - double dw_dE[9] = {0.0}; - double d2W_d2E[81] = {0.0}; - double h = 1e-7; - - run_fwd_mode(E, mu, lambda, dw_dE, N); - run_bkwd_mode(E, mu, lambda, dw_dE, N); - run_hand_derivative(E, mu, lambda, dw_dE, N); - run_fwd_fwd(E, mu, lambda, d2W_d2E, N); - run_rev_fwd(E, mu, lambda, d2W_d2E, N); - run_fwd_rev(E, mu, lambda, d2W_d2E, N); - run_rev_rev(E, mu, lambda, d2W_d2E, N); - run_fwd_FD(E, mu, lambda, d2W_d2E, N, h); - run_hand_fwd(E, mu, lambda, d2W_d2E, N); - - - return 0; - -} - -// TEST(StrainEnergyTest, StressFiniteDifferenceVsAutodiff) { -// double F[9] = {1.0, 0.5, 0.0, 0.0, 1.2, 0.1, 0.0, 0.0, 1.0}; -// double E[9] = {0.0}; -// double mu = 1.0; -// double lambda = 1.0; -// double dW_dF_fd[9]; -// double dW_dF_ad[9]; -// double J = calc_determinant(F); -// calc_green_lagrange(F, E); -// std::cout << "Calling stress_FD..." << std::endl; -// stress_FD(E, J, mu, lambda, dW_dF_fd); -// std::cout << "Finite difference stress computed:" << std::endl; -// for (int i = 0; i < 9; ++i) { -// std::cout << "dW_dF_fd[" << i << "] = " << dW_dF_fd[i] << std::endl; -// } - -// std::cout << "Calling stress..." << std::endl; -// stress(E, J, mu, lambda, dW_dF_ad); -// std::cout << "Autodiff stress computed:" << std::endl; -// for (int i = 0; i < 9; ++i) { -// std::cout << "dW_dF_ad[" << i << "] = " << dW_dF_ad[i] << std::endl; -// } - -// for (int i = 0; i < 9; ++i) { -// std::cout << "Comparing index " << i << ": FD = " << dW_dF_fd[i] -// << ", AD = " << dW_dF_ad[i] << std::endl; -// EXPECT_NEAR(dW_dF_fd[i], dW_dF_ad[i], 1e-6) << "Mismatch at index " << i; -// } -// } - -// TEST(StrainEnergyTest, StressFiniteDifferenceVsReverseAutodiff) { -// double E[9] = {1.0, 0.5, 0.0, 0.0, 1.2, 0.1, 0.0, 0.0, 1.0}; -// double mu = 1.0; -// double lambda = 1.0; -// double dW_dF_fd[9]; -// double dW_dF_rev[9]; -// std::cout << "E: "; -// for (int i = 0; i < 9; ++i) std::cout << E[i] << " "; -// std::cout << std::endl; -// std::cout << "Calling stress_FD..." << std::endl; -// stress_FD(E, mu, lambda, dW_dF_fd); -// std::cout << "Finite difference stress computed:" << std::endl; -// for (int i = 0; i < 9; ++i) { -// std::cout << "dW_dF_fd[" << i << "] = " << dW_dF_fd[i] << std::endl; -// } - -// std::cout << "Calling stress_reverse..." << std::endl; -// stress_reverse(E, mu, lambda, dW_dF_rev); -// std::cout << "Reverse-mode autodiff stress computed:" << std::endl; -// for (int i = 0; i < 9; ++i) { -// std::cout << "dW_dF_rev[" << i << "] = " << dW_dF_rev[i] << std::endl; -// } - -// for (int i = 0; i < 9; ++i) { -// std::cout << "Comparing index " << i << ": FD = " << dW_dF_fd[i] -// << ", Reverse AD = " << dW_dF_rev[i] << std::endl; -// EXPECT_NEAR(dW_dF_fd[i], dW_dF_rev[i], 1e-6) << "Mismatch at index " << i; -// } -// } - -// TEST(StrainEnergyTest, StressFiniteDifferenceVsReverseAutodiff) { -// double E[9] = {0.1, 0.05, 0.02, 0.05, 0.2, 0.01, 0.02, 0.01, 0.15}; -// double mu = 1.0; -// double lambda = 1.0; -// double dW_dF_fd[9]; -// double dW_dF_hand[9]; -// std::cout << "E: "; -// for (int i = 0; i < 9; ++i) std::cout << E[i] << " "; -// std::cout << std::endl; -// std::cout << "Calling stress_FD..." << std::endl; -// stress_FD(E, mu, lambda, dW_dF_fd); -// std::cout << "Finite difference stress computed:" << std::endl; -// for (int i = 0; i < 9; ++i) { -// std::cout << "dW_dF_fd[" << i << "] = " << dW_dF_fd[i] << std::endl; -// } - -// std::cout << "Calling hand_code_deriv..." << std::endl; -// hand_code_deriv(E, mu, lambda, dW_dF_hand); -// std::cout << "Hand coded stress computed:" << std::endl; -// for (int i = 0; i < 9; ++i) { -// std::cout << "dW_dF_hand[" << i << "] = " << dW_dF_hand[i] << std::endl; -// } - -// for (int i = 0; i < 9; ++i) { -// std::cout << "Comparing index " << i << ": FD = " << dW_dF_fd[i] -// << ", Hand derivative = " << dW_dF_hand[i] << std::endl; -// EXPECT_NEAR(dW_dF_fd[i], dW_dF_hand[i], 1e-6) << "Mismatch at index " << i; -// } -// } - -// TEST(StrainEnergyTest, StressFiniteDifferenceVsReverseAutodiff) { -// double E[9] = {0.1, 0.05, 0.02, 0.05, 0.2, 0.01, 0.02, 0.01, 0.15}; -// double mu = 1.0; -// double lambda = 1.0; -// double dW_dF_fwd_fd[81]; -// double dW_dF_fwd_fwd[81]; -// std::cout << "E: "; -// for (int i = 0; i < 81; ++i) std::cout << E[i] << " "; -// std::cout << std::endl; -// std::cout << "Calling stress_FD..." << std::endl; -// second_deriv_fwd_FD(E, mu, lambda, dW_dF_fwd_fd); -// std::cout << "Finite difference stress computed:" << std::endl; -// for (int i = 0; i < 81; ++i) { -// std::cout << "dW_dF_fd[" << i << "] = " << dW_dF_fwd_fd[i] << std::endl; -// } - -// std::cout << "Calling hand_code_deriv..." << std::endl; -// second_deriv_fwd_fwd(E, mu, lambda, dW_dF_fwd_fwd); -// std::cout << "Hand coded stress computed:" << std::endl; -// for (int i = 0; i < 81; ++i) { -// std::cout << "dW_dF_hand[" << i << "] = " << dW_dF_fwd_fwd[i] << std::endl; -// } - -// for (int i = 0; i < 81; ++i) { -// std::cout << "Comparing index " << i << ": FD = " << dW_dF_fwd_fd[i] -// << ", Hand derivative = " << dW_dF_fwd_fwd[i] << std::endl; -// EXPECT_NEAR(dW_dF_fwd_fd[i], dW_dF_fwd_fwd[i], 1e-6) << "Mismatch at index " << i; -// } -// } - -// TEST(StrainEnergyTest, StressFiniteDifferenceVsReverseAutodiff) { -// double E[9] = {0.1, 0.05, 0.02, 0.05, 0.2, 0.01, 0.02, 0.01, 0.15}; -// double mu = 1.0; -// double lambda = 1.0; -// double dW_dF_fwd_fd[81]; -// double dW_dF_fwd_rev[81]; - -// std::cout << "E: "; -// for (int i = 0; i < 9; ++i) std::cout << E[i] << " "; -// std::cout << std::endl; - -// std::cout << "Calling stress_FD..." << std::endl; -// second_deriv_fwd_FD(E, mu, lambda, dW_dF_fwd_fd); -// std::cout << "Finite difference stress computed:" << std::endl; -// for (int i = 0; i < 81; ++i) { -// std::cout << "dW_dF_fd[" << i << "] = " << dW_dF_fwd_fd[i] << std::endl; -// } - -// std::cout << "Calling hand_code_deriv (fwd_rev)..." << std::endl; -// second_deriv_fwd_rev(E, mu, lambda, dW_dF_fwd_rev); -// std::cout << "Hand coded stress computed (fwd_rev):" << std::endl; -// for (int i = 0; i < 81; ++i) { -// std::cout << "dW_dF_hand[" << i << "] = " << dW_dF_fwd_rev[i] << std::endl; -// } - -// for (int i = 0; i < 81; ++i) { -// std::cout << "Comparing index " << i << ": FD = " << dW_dF_fwd_fd[i] -// << ", Hand derivative = " << dW_dF_fwd_rev[i] << std::endl; -// EXPECT_NEAR(dW_dF_fwd_fd[i], dW_dF_fwd_rev[i], 1e-6) << "Mismatch at index " << i; -// } -// } - -// TEST(StrainEnergyTest, StressFiniteDifferenceVsReverseAutodiff) { -// double E[9] = {0.1, 0.05, 0.02, 0.05, 0.2, 0.01, 0.02, 0.01, 0.15}; -// double mu = 1.0; -// double lambda = 1.0; -// double dW_dF_fwd_fd[81]; -// double dW_dF_rev_rev[81]; - -// std::cout << "E: "; -// for (int i = 0; i < 9; ++i) std::cout << E[i] << " "; -// std::cout << std::endl; - -// std::cout << "Calling stress_FD..." << std::endl; -// second_deriv_fwd_FD(E, mu, lambda, dW_dF_fwd_fd); -// std::cout << "Finite difference stress computed:" << std::endl; -// for (int i = 0; i < 81; ++i) { -// std::cout << "dW_dF_fd[" << i << "] = " << dW_dF_fwd_fd[i] << std::endl; -// } - -// std::cout << "Calling hand_code_deriv (rev_rev)..." << std::endl; -// second_deriv_rev_rev(E, mu, lambda, dW_dF_rev_rev); -// std::cout << "Hand coded stress computed (rev_rev):" << std::endl; -// for (int i = 0; i < 81; ++i) { -// std::cout << "dW_dF_hand[" << i << "] = " << dW_dF_rev_rev[i] << std::endl; -// } - -// for (int i = 0; i < 81; ++i) { -// std::cout << "Comparing index " << i << ": FD = " << dW_dF_fwd_fd[i] -// << ", Hand derivative = " << dW_dF_rev_rev[i] << std::endl; -// EXPECT_NEAR(dW_dF_fwd_fd[i], dW_dF_rev_rev[i], 1e-6) << "Mismatch at index " << i; -// } -// } - -// TEST(StrainEnergyTest, StressFiniteDifferenceVsFwdRevAutodiff) { -// double E[9] = {0.1, 0.05, 0.02, 0.05, 0.2, 0.01, 0.02, 0.01, 0.15}; -// double mu = 1.0; -// double lambda = 1.0; -// double dW_dF_fwd_fd[81]; -// double dW_dF_fwd_rev[81]; - -// std::cout << "E: "; -// for (int i = 0; i < 9; ++i) std::cout << E[i] << " "; -// std::cout << std::endl; - -// std::cout << "Calling stress_FD..." << std::endl; -// second_deriv_fwd_FD(E, mu, lambda, dW_dF_fwd_fd); -// std::cout << "Finite difference stress computed:" << std::endl; -// for (int i = 0; i < 81; ++i) { -// std::cout << "dW_dF_fd[" << i << "] = " << dW_dF_fwd_fd[i] << std::endl; -// } - -// std::cout << "Calling hand_code_deriv (fwd_rev)..." << std::endl; -// second_deriv_fwd_rev(E, mu, lambda, dW_dF_fwd_rev); -// std::cout << "Hand coded stress computed (fwd_rev):" << std::endl; -// for (int i = 0; i < 81; ++i) { -// std::cout << "dW_dF_hand[" << i << "] = " << dW_dF_fwd_rev[i] << std::endl; -// } - -// for (int i = 0; i < 81; ++i) { -// std::cout << "Comparing index " << i << ": FD = " << dW_dF_fwd_fd[i] -// << ", Hand derivative = " << dW_dF_fwd_rev[i] << std::endl; -// EXPECT_NEAR(dW_dF_fwd_fd[i], dW_dF_fwd_rev[i], 1e-6) << "Mismatch at index " << i; -// } -// } - - -// TEST(StrainEnergyTest, StressFiniteDifferenceVsFwdRevAutodiff) { -// double E[9] = {0.12, -0.03, 0.01, -0.03, 0.08, 0.02, 0.01, 0.02, 0.11}; -// double mu = 1.0; -// double lambda = 1.0; -// double dW_dF_fwd_fd[81]; -// double dW_dF_fwd_rev[81]; - -// std::cout << "E: "; -// for (int i = 0; i < 9; ++i) std::cout << E[i] << " "; -// std::cout << std::endl; - -// std::cout << "Calling stress_FD..." << std::endl; -// second_deriv_fwd_FD(E, mu, lambda, dW_dF_fwd_fd); -// std::cout << "Finite difference stress computed:" << std::endl; -// for (int i = 0; i < 81; ++i) { -// std::cout << "dW_dF_fd[" << i << "] = " << dW_dF_fwd_fd[i] << std::endl; -// } - -// std::cout << "Calling hand_code_deriv (fwd_rev)..." << std::endl; -// second_deriv_hand_fwd(E, mu, lambda, dW_dF_fwd_rev); -// std::cout << "Hand coded stress computed (fwd_rev):" << std::endl; -// for (int i = 0; i < 81; ++i) { -// std::cout << "dW_dF_hand[" << i << "] = " << dW_dF_fwd_rev[i] << std::endl; -// } - -// for (int i = 0; i < 81; ++i) { -// std::cout << "Comparing index " << i << ": FD = " << dW_dF_fwd_fd[i] -// << ", Hand derivative = " << dW_dF_fwd_rev[i] << std::endl; -// EXPECT_NEAR(dW_dF_fwd_fd[i], dW_dF_fwd_rev[i], 1e-6) << "Mismatch at index " << i; -// } -// } - -TEST(StrainEnergyTest, StressFiniteDifferenceVsFwdRevAutodiff) { - double F[9] = {1.01, -0.03, 0.01, -0.03, 1.05, 0.02, 0.01, 0.02, 1.0}; - double E[9] = {0}; - double mu = 1.0; - double lambda = 1.0; - double* S = nullptr; - calc_E_from_F(F, E); - hand_code_deriv(E, mu, lambda, S); - std::cout << "{ "; - for (int i = 0; i < 9; ++i) { - std::cout << ", " << S[i]; - } - std::cout << " }" << std::endl; - std::cout << "E: { "; - for (int i = 0; i < 9; ++i) { - std::cout << ", " << E[i]; - } - std::cout << "}" << std::endl; - - double dW_dF_hand_fd[81]; - double dW_dF_hand_fwd[81]; - - std::cout << "E: "; - for (int i = 0; i < 9; ++i) std::cout << E[i] << " "; - std::cout << std::endl; - - std::cout << "Calling stress_FD..." << std::endl; - second_deriv_hand_FD(E, mu, lambda, dW_dF_hand_fd); - std::cout << "Finite difference stress computed:" << std::endl; - for (int i = 0; i < 81; ++i) { - std::cout << "dW_dF_hand_fd[" << i << "] = " << dW_dF_hand_fd[i] << std::endl; - } - - std::cout << "Calling hand_code_deriv (fwd_rev)..." << std::endl; - second_deriv_hand_fwd(E, mu, lambda, dW_dF_hand_fwd); - std::cout << "Hand coded stress computed (fwd_rev):" << std::endl; - for (int i = 0; i < 81; ++i) { - std::cout << "dW_dF_hand_fwd[" << i << "] = " << dW_dF_hand_fwd[i] << std::endl; - } - - for (int i = 0; i < 81; ++i) { - std::cout << "Comparing index " << i << ": FD = " << dW_dF_hand_fd[i] - << ", Hand derivative = " << dW_dF_hand_fwd[i] << std::endl; - EXPECT_NEAR(dW_dF_hand_fd[i], dW_dF_hand_fd[i], 1e-6) << "Mismatch at index " << i; - } -} \ No newline at end of file diff --git a/src/tests/enzyme_neohookean.cpp b/src/tests/enzyme_neohookean.cpp index 558490fa..7c85b3ce 100644 --- a/src/tests/enzyme_neohookean.cpp +++ b/src/tests/enzyme_neohookean.cpp @@ -8,7 +8,6 @@ #include #include "tribol/common/Enzyme.hpp" - class EnzymeNeohookeanTest : public testing::Test { protected: double F_all[27] = { 1.0, 0.5, 0.0, 0.0, 1.2, 0.1, 0.0, 0.0, 1.0, 1.1, 0.7, 0.1, 0.2, 1.4, @@ -127,9 +126,8 @@ TEST_F( EnzymeNeohookeanTest, 2nd_deriv_finite_difference_vs_hand_fwd ) double dW2_hand_fwd[81] = { 0.0 }; second_deriv_fwd_FD( E_all + i * 9, mu, lambda, dW2_FD ); second_deriv_hand_fwd( E_all + i * 9, mu, lambda, dW2_hand_fwd ); - for(int j = 0; j < 9; ++j) { - EXPECT_NEAR( dW2_FD[j], dW2_hand_fwd[j], tol ); - } + for ( int j = 0; j < 9; ++j ) { + EXPECT_NEAR( dW2_FD[j], dW2_hand_fwd[j], tol ); + } } } - diff --git a/src/tests/neohookean_common.hpp b/src/tests/neohookean_common.hpp index 32df6730..766a0370 100644 --- a/src/tests/neohookean_common.hpp +++ b/src/tests/neohookean_common.hpp @@ -5,6 +5,8 @@ #include #include #include "tribol/common/Enzyme.hpp" +#include "axom/core.hpp" +#include "axom/slic/interface/slic.hpp" void multiply3x3( const double F[9], const double F_T[9], double C[9] ) { @@ -290,4 +292,112 @@ void second_deriv_hand_FD( double* E, double mu, double lambda, double* d2W_d2E, d2W_d2E[9 * i + j] = ( dw_plus[j] - dw_minus[j] ) / ( 2 * h ); } } +} + + +void run_fwd_mode(double* E, double mu, double lambda, double* dw_df, int N) { + axom::utilities::Timer timer{ false }; + stress(E, mu, lambda, dw_df); + double Dw[9] = {0.0}; + timer.start(); + for(int i = 0; i < N; ++i) { + stress(E, mu, lambda, dw_df); + for(int j = 0; j < 9; ++j) { + Dw[j] += dw_df[j]; + } + } + timer.stop(); + std::cout << axom::fmt::format( "Time to calc fwd_diff: {0:f}ms", timer.elapsedTimeInMilliSec() ) << std::endl; +} + +void run_bkwd_mode(double* E, double mu, double lambda, double* dw_dE, int N) { + axom::utilities::Timer timer{ false }; + + stress_reverse(E, mu, lambda, dw_dE); + double Dw[9] = {0.0}; + + timer.start(); + for(int i = 0; i < N; ++i) { + stress_reverse(E, mu, lambda, dw_dE); + for(int j = 0; j < 9; ++j) { + Dw[j] += dw_dE[j]; + } + } + timer.stop(); + std::cout << axom::fmt::format( "Time to calc backward_diff: {0:f}ms", timer.elapsedTimeInMilliSec() ) << std::endl; +} + +void run_hand_derivative(double* E, double mu, double lambda, double* S, int N) { + axom::utilities::Timer timer{ false }; + hand_code_deriv(E, mu, lambda, S); + double Dw[9] = {0.0}; + timer.start(); + for(int i = 0; i < N; ++i) { + hand_code_deriv(E, mu, lambda, S); + for(int j = 0; j < 9; ++j) { + Dw[j] += S[j]; + } + } + timer.stop(); + std::cout << axom::fmt::format( "Time to calc hand_diff: {0:f}ms", timer.elapsedTimeInMilliSec() ) << std::endl; +} + +void run_fwd_fwd(double* E, double mu, double lambda, double* S, int N) { + axom::utilities::Timer timer{ false }; + timer.start(); + for(int i = 0; i < N; ++i) { + second_deriv_fwd_fwd(E, mu, lambda, S); + } + timer.stop(); + std::cout << axom::fmt::format( "Time to calc fwd_fwd_diff: {0:f}ms", timer.elapsedTimeInMilliSec() ) << std::endl; +} + +void run_rev_fwd(double* E, double mu, double lambda, double* S, int N) { + axom::utilities::Timer timer{ false }; + timer.start(); + for(int i = 0; i < N; ++i) { + second_deriv_rev_fwd(E, mu, lambda, S); + } + timer.stop(); + std::cout << axom::fmt::format( "Time to calc rev_fwd_diff: {0:f}ms", timer.elapsedTimeInMilliSec() ) << std::endl; +} + +void run_fwd_rev(double* E, double mu, double lambda, double* S, int N) { + axom::utilities::Timer timer{ false }; + timer.start(); + for(int i = 0; i < N; ++i) { + second_deriv_fwd_rev(E, mu, lambda, S); + } + timer.stop(); + std::cout << axom::fmt::format( "Time to calc fwd_rev_diff: {0:f}ms", timer.elapsedTimeInMilliSec() ) << std::endl; +} + +void run_rev_rev(double* E, double mu, double lambda, double* S, int N) { + axom::utilities::Timer timer{ false }; + timer.start(); + for(int i = 0; i < N; ++i) { + second_deriv_rev_rev(E, mu, lambda, S); + } + timer.stop(); + std::cout << axom::fmt::format( "Time to calc rev_rev_diff: {0:f}ms", timer.elapsedTimeInMilliSec() ) << std::endl; +} + +void run_fwd_FD(double* E, double mu, double lambda, double* S, int N, double h) { + axom::utilities::Timer timer{ false }; + timer.start(); + for(int i = 0; i < N; ++i) { + second_deriv_fwd_FD(E, mu, lambda, S, h); + } + timer.stop(); + std::cout << axom::fmt::format( "Time to calc fwd_FD_diff: {0:f}ms", timer.elapsedTimeInMilliSec() ) << std::endl; +} + +void run_hand_fwd(double* E, double mu, double lambda, double* S, int N) { + axom::utilities::Timer timer{ false }; + timer.start(); + for(int i = 0; i < N; ++i) { + second_deriv_hand_fwd(E, mu, lambda, S); + } + timer.stop(); + std::cout << axom::fmt::format( "Time to calc hand_fwd_diff: {0:f}ms", timer.elapsedTimeInMilliSec() ) << std::endl; } \ No newline at end of file From 0f62d2cda11814e6a2fbc71ba9c669eaaaad6af7 Mon Sep 17 00:00:00 2001 From: Ryan Lutz Date: Wed, 18 Jun 2025 16:41:27 -0700 Subject: [PATCH 4/4] cleaned up code/ added lambda functions to reduce --- src/examples/neohookean_timing.cpp | 100 ++++++++++++++++++++---- src/tests/enzyme_neohookean.cpp | 117 +++++++++++++---------------- src/tests/neohookean_common.hpp | 106 -------------------------- 3 files changed, 137 insertions(+), 186 deletions(-) diff --git a/src/examples/neohookean_timing.cpp b/src/examples/neohookean_timing.cpp index 8aee35f7..244ba132 100644 --- a/src/examples/neohookean_timing.cpp +++ b/src/examples/neohookean_timing.cpp @@ -2,13 +2,49 @@ #include "../tests/neohookean_common.hpp" #include #include -#include #include "axom/core.hpp" #include "axom/slic/interface/slic.hpp" -#include -#include #include "tribol/common/Enzyme.hpp" +template +void do_dWdE_timing( double* E, double mu, double lambda, int N, TangentFn&& tangent_fn) +{ + double S[9] = {0.0}; + axom::utilities::Timer timer{ false }; + tangent_fn(E, mu, lambda, S); + + + timer.start(); + for(int i = 0; i < N; ++i) { + tangent_fn(E, mu, lambda, S); + // for(int j = 0; j < 9; ++j) { + // Dw[j] += S[j]; + // } + } + timer.stop(); + std::cout << axom::fmt::format( "{:f}ms", timer.elapsedTimeInMilliSec() ) << std::endl; +} + +template +void do_dSdE_timing( double* E, double mu, double lambda, int N, TangentFn&& tangent_fn) +{ + double S[81] = {0.0}; + axom::utilities::Timer timer{ false }; + tangent_fn(E, mu, lambda, S); + + + timer.start(); + for(int i = 0; i < N; ++i) { + tangent_fn(E, mu, lambda, S); + // for(int j = 0; j < 9; ++j) { + // Dw[j] += S[j]; + // } + } + timer.stop(); + std::cout << axom::fmt::format( "{:f}ms", timer.elapsedTimeInMilliSec() ) << std::endl; +} + + int main() { double E[9] = { 0.1, 0.05, 0.02, 0.05, 0.2, 0.01, 0.02, 0.01, 0.15 }; @@ -17,19 +53,51 @@ int main() int N = 0; std::cout << "Enter cycles: "; std::cin >> N; - double dw_dE[9] = { 0.0 }; - double d2W_d2E[81] = { 0.0 }; - double h = 1e-7; - - run_fwd_mode( E, mu, lambda, dw_dE, N ); - run_bkwd_mode( E, mu, lambda, dw_dE, N ); - run_hand_derivative( E, mu, lambda, dw_dE, N ); - run_fwd_fwd( E, mu, lambda, d2W_d2E, N ); - run_rev_fwd( E, mu, lambda, d2W_d2E, N ); - run_fwd_rev( E, mu, lambda, d2W_d2E, N ); - run_rev_rev( E, mu, lambda, d2W_d2E, N ); - run_fwd_FD( E, mu, lambda, d2W_d2E, N, h ); - run_hand_fwd( E, mu, lambda, d2W_d2E, N ); + std::cout << "Time to calc stress with fwddiff: "; + + do_dWdE_timing(E, mu, lambda, N, [](double * E, double mu, double lambda, double* S) { + stress(E, mu, lambda, S); + }); + + std::cout << "Time to calc stress with autodiff: "; + do_dWdE_timing(E, mu, lambda, N, [](double * E, double mu, double lambda, double* S) { + stress_reverse(E, mu, lambda, S); + }); + + std::cout << "Time to calc stress with hand coded deriv: "; + do_dWdE_timing(E, mu, lambda, N, [](double * E, double mu, double lambda, double* S) { + hand_code_deriv(E, mu, lambda, S); + }); + + std::cout << "Time to calc 2nd deriv with fwd_fwd: "; + do_dSdE_timing(E, mu, lambda, N, [](double * E, double mu, double lambda, double* S) { + second_deriv_fwd_fwd(E, mu, lambda, S); + }); + + std::cout << "Time to calc 2nd deriv with rev_fwd: "; + do_dSdE_timing(E, mu, lambda, N, [](double * E, double mu, double lambda, double* S) { + second_deriv_rev_fwd(E, mu, lambda, S); + }); + + std::cout << "Time to calc 2nd deriv with fwd_rev: "; + do_dSdE_timing(E, mu, lambda, N, [](double * E, double mu, double lambda, double* S) { + second_deriv_fwd_rev(E, mu, lambda, S); + }); + + std::cout << "Time to calc 2nd deriv with rev_rev: "; + do_dSdE_timing(E, mu, lambda, N, [](double * E, double mu, double lambda, double* S) { + second_deriv_rev_rev(E, mu, lambda, S); + }); + + std::cout << "Time to calc 2nd deriv with fwd_FD: "; + do_dSdE_timing(E, mu, lambda, N, [](double * E, double mu, double lambda, double* S) { + second_deriv_fwd_FD(E, mu, lambda, S); + }); + + std::cout << "Time to calc 2nd deriv with hand_fwd: "; + do_dSdE_timing(E, mu, lambda, N, [](double * E, double mu, double lambda, double* S) { + second_deriv_hand_fwd(E, mu, lambda, S); + }); return 0; } diff --git a/src/tests/enzyme_neohookean.cpp b/src/tests/enzyme_neohookean.cpp index 7c85b3ce..cdc6f6ae 100644 --- a/src/tests/enzyme_neohookean.cpp +++ b/src/tests/enzyme_neohookean.cpp @@ -3,16 +3,14 @@ #include "neohookean_common.hpp" #include #include -#include -#include -#include -#include "tribol/common/Enzyme.hpp" + class EnzymeNeohookeanTest : public testing::Test { protected: +//clang-format off double F_all[27] = { 1.0, 0.5, 0.0, 0.0, 1.2, 0.1, 0.0, 0.0, 1.0, 1.1, 0.7, 0.1, 0.2, 1.4, 0.2, 0.02, 0.04, 1.05, 1.05, 0.2, 0.04, 0.0, 1.1, 0.14, 0.07, 0.0, 1.09 }; - +//clang-format on double mu = 1.6; double lambda = 1.2; @@ -28,95 +26,86 @@ class EnzymeNeohookeanTest : public testing::Test { } }; -TEST_F( EnzymeNeohookeanTest, finite_difference_vs_exact ) + +template +void test_stress( double* E_all, double mu, double lambda, double tol, StressFn&& stress_test_fn) { for ( int i = 0; i < 3; ++i ) { double S_FD[9] = { 0.0 }; - double S_exact[9] = { 0.0 }; + double S_test[9] = { 0.0 }; stress_FD( E_all + i * 9, mu, lambda, S_FD ); - hand_code_deriv( E_all + i * 9, mu, lambda, S_exact ); + stress_test_fn( E_all + i * 9, mu, lambda, S_test ); for ( int j = 0; j < 9; ++j ) { - EXPECT_NEAR( S_FD[j], S_exact[j], tol ); + EXPECT_NEAR( S_FD[j], S_test[j], tol ); } } } -TEST_F( EnzymeNeohookeanTest, finite_difference_vs_fwd_diff ) + +template +void test_stress_2nd_deriv( double* E_all, double mu, double lambda, double tol, StressFn&& deriv2_test_fn) { for ( int i = 0; i < 3; ++i ) { - double S_FD[9] = { 0.0 }; - double S_fwddiff[9] = { 0.0 }; - stress_FD( E_all + i * 9, mu, lambda, S_FD ); - stress( E_all + i * 9, mu, lambda, S_fwddiff ); - for ( int j = 0; j < 9; ++j ) { - EXPECT_NEAR( S_FD[j], S_fwddiff[j], tol ); + double dW2_FD[81] = { 0.0 }; + double dW2_test[81] = { 0.0 }; + second_deriv_fwd_FD( E_all + i * 9, mu, lambda, dW2_FD ); + deriv2_test_fn( E_all + i * 9, mu, lambda, dW2_test ); + for ( int j = 0; j < 81; ++j ) { + EXPECT_NEAR( dW2_FD[j], dW2_test[j], tol ); } } } + + +TEST_F( EnzymeNeohookeanTest, finite_difference_vs_exact ) +{ + test_stress(E_all, mu, lambda, tol, [](double* E, double mu, double lambda, double* S) { + hand_code_deriv(E, mu, lambda, S); + }); +} + +TEST_F( EnzymeNeohookeanTest, finite_difference_vs_fwd_diff ) +{ + test_stress(E_all, mu, lambda, tol, [](double* E, double mu, double lambda, double* S) { + stress(E, mu, lambda, S); + }); +} + + TEST_F( EnzymeNeohookeanTest, finite_difference_vs_autodiff ) { - for ( int i = 0; i < 3; ++i ) { - double S_FD[9] = { 0.0 }; - double S_autodiff[9] = { 0.0 }; - stress_FD( E_all + i * 9, mu, lambda, S_FD ); - stress_reverse( E_all + i * 9, mu, lambda, S_autodiff ); - for ( int j = 0; j < 9; ++j ) { - EXPECT_NEAR( S_FD[j], S_autodiff[j], tol ); - } - } + test_stress(E_all, mu, lambda, tol, [](double* E, double mu, double lambda, double* S) { + stress_reverse(E, mu, lambda, S); + }); } TEST_F( EnzymeNeohookeanTest, 2nd_deriv_finite_difference_vs_fwd_fwd ) { - for ( int i = 0; i < 3; ++i ) { - double dW2_FD[81] = { 0.0 }; - double dW2_fwd_fwd[81] = { 0.0 }; - second_deriv_fwd_FD( E_all + i * 9, mu, lambda, dW2_FD ); - second_deriv_fwd_fwd( E_all + i * 9, mu, lambda, dW2_fwd_fwd ); - for ( int j = 0; j < 9; ++j ) { - EXPECT_NEAR( dW2_FD[j], dW2_fwd_fwd[j], tol ); - } - } + test_stress_2nd_deriv(E_all, mu, lambda, tol, [](double* E, double mu, double lambda, double* S) { + second_deriv_fwd_fwd(E, mu, lambda, S); + }); } TEST_F( EnzymeNeohookeanTest, 2nd_deriv_finite_difference_vs_fwd_rev ) -{ - for ( int i = 0; i < 3; ++i ) { - double dW2_FD[81] = { 0.0 }; - double dW2_fwd_rev[81] = { 0.0 }; - second_deriv_fwd_FD( E_all + i * 9, mu, lambda, dW2_FD ); - second_deriv_fwd_rev( E_all + i * 9, mu, lambda, dW2_fwd_rev ); - for ( int j = 0; j < 81; ++j ) { - EXPECT_NEAR( dW2_FD[j], dW2_fwd_rev[j], tol ); - } - } + { + test_stress_2nd_deriv(E_all, mu, lambda, tol, [](double* E, double mu, double lambda, double* S) { + second_deriv_fwd_rev(E, mu, lambda, S); + }); } TEST_F( EnzymeNeohookeanTest, 2nd_deriv_finite_difference_vs_rev_rev ) -{ - for ( int i = 0; i < 3; ++i ) { - double dW2_FD[81] = { 0.0 }; - double dW2_rev_rev[81] = { 0.0 }; - second_deriv_fwd_FD( E_all + i * 9, mu, lambda, dW2_FD ); - second_deriv_rev_rev( E_all + i * 9, mu, lambda, dW2_rev_rev ); - for ( int j = 0; j < 81; ++j ) { - EXPECT_NEAR( dW2_FD[j], dW2_rev_rev[j], tol ); - } - } + { + test_stress_2nd_deriv(E_all, mu, lambda, tol, [](double* E, double mu, double lambda, double* S) { + second_deriv_rev_rev(E, mu, lambda, S); + }); } TEST_F( EnzymeNeohookeanTest, 2nd_deriv_finite_difference_vs_rev_fwd ) -{ - for ( int i = 0; i < 3; ++i ) { - double dW2_FD[81] = { 0.0 }; - double dW2_rev_fwd[81] = { 0.0 }; - second_deriv_fwd_FD( E_all + i * 9, mu, lambda, dW2_FD ); - second_deriv_rev_rev( E_all + i * 9, mu, lambda, dW2_rev_fwd ); - for ( int j = 0; j < 81; ++j ) { - EXPECT_NEAR( dW2_FD[j], dW2_rev_fwd[j], tol ); - } - } + { + test_stress_2nd_deriv(E_all, mu, lambda, tol, [](double* E, double mu, double lambda, double* S) { + second_deriv_rev_fwd(E, mu, lambda, S); + }); } TEST_F( EnzymeNeohookeanTest, 2nd_deriv_finite_difference_vs_hand_fwd ) diff --git a/src/tests/neohookean_common.hpp b/src/tests/neohookean_common.hpp index 766a0370..959f379d 100644 --- a/src/tests/neohookean_common.hpp +++ b/src/tests/neohookean_common.hpp @@ -295,109 +295,3 @@ void second_deriv_hand_FD( double* E, double mu, double lambda, double* d2W_d2E, } -void run_fwd_mode(double* E, double mu, double lambda, double* dw_df, int N) { - axom::utilities::Timer timer{ false }; - stress(E, mu, lambda, dw_df); - double Dw[9] = {0.0}; - timer.start(); - for(int i = 0; i < N; ++i) { - stress(E, mu, lambda, dw_df); - for(int j = 0; j < 9; ++j) { - Dw[j] += dw_df[j]; - } - } - timer.stop(); - std::cout << axom::fmt::format( "Time to calc fwd_diff: {0:f}ms", timer.elapsedTimeInMilliSec() ) << std::endl; -} - -void run_bkwd_mode(double* E, double mu, double lambda, double* dw_dE, int N) { - axom::utilities::Timer timer{ false }; - - stress_reverse(E, mu, lambda, dw_dE); - double Dw[9] = {0.0}; - - timer.start(); - for(int i = 0; i < N; ++i) { - stress_reverse(E, mu, lambda, dw_dE); - for(int j = 0; j < 9; ++j) { - Dw[j] += dw_dE[j]; - } - } - timer.stop(); - std::cout << axom::fmt::format( "Time to calc backward_diff: {0:f}ms", timer.elapsedTimeInMilliSec() ) << std::endl; -} - -void run_hand_derivative(double* E, double mu, double lambda, double* S, int N) { - axom::utilities::Timer timer{ false }; - hand_code_deriv(E, mu, lambda, S); - double Dw[9] = {0.0}; - timer.start(); - for(int i = 0; i < N; ++i) { - hand_code_deriv(E, mu, lambda, S); - for(int j = 0; j < 9; ++j) { - Dw[j] += S[j]; - } - } - timer.stop(); - std::cout << axom::fmt::format( "Time to calc hand_diff: {0:f}ms", timer.elapsedTimeInMilliSec() ) << std::endl; -} - -void run_fwd_fwd(double* E, double mu, double lambda, double* S, int N) { - axom::utilities::Timer timer{ false }; - timer.start(); - for(int i = 0; i < N; ++i) { - second_deriv_fwd_fwd(E, mu, lambda, S); - } - timer.stop(); - std::cout << axom::fmt::format( "Time to calc fwd_fwd_diff: {0:f}ms", timer.elapsedTimeInMilliSec() ) << std::endl; -} - -void run_rev_fwd(double* E, double mu, double lambda, double* S, int N) { - axom::utilities::Timer timer{ false }; - timer.start(); - for(int i = 0; i < N; ++i) { - second_deriv_rev_fwd(E, mu, lambda, S); - } - timer.stop(); - std::cout << axom::fmt::format( "Time to calc rev_fwd_diff: {0:f}ms", timer.elapsedTimeInMilliSec() ) << std::endl; -} - -void run_fwd_rev(double* E, double mu, double lambda, double* S, int N) { - axom::utilities::Timer timer{ false }; - timer.start(); - for(int i = 0; i < N; ++i) { - second_deriv_fwd_rev(E, mu, lambda, S); - } - timer.stop(); - std::cout << axom::fmt::format( "Time to calc fwd_rev_diff: {0:f}ms", timer.elapsedTimeInMilliSec() ) << std::endl; -} - -void run_rev_rev(double* E, double mu, double lambda, double* S, int N) { - axom::utilities::Timer timer{ false }; - timer.start(); - for(int i = 0; i < N; ++i) { - second_deriv_rev_rev(E, mu, lambda, S); - } - timer.stop(); - std::cout << axom::fmt::format( "Time to calc rev_rev_diff: {0:f}ms", timer.elapsedTimeInMilliSec() ) << std::endl; -} - -void run_fwd_FD(double* E, double mu, double lambda, double* S, int N, double h) { - axom::utilities::Timer timer{ false }; - timer.start(); - for(int i = 0; i < N; ++i) { - second_deriv_fwd_FD(E, mu, lambda, S, h); - } - timer.stop(); - std::cout << axom::fmt::format( "Time to calc fwd_FD_diff: {0:f}ms", timer.elapsedTimeInMilliSec() ) << std::endl; -} - -void run_hand_fwd(double* E, double mu, double lambda, double* S, int N) { - axom::utilities::Timer timer{ false }; - timer.start(); - for(int i = 0; i < N; ++i) { - second_deriv_hand_fwd(E, mu, lambda, S); - } - timer.stop(); - std::cout << axom::fmt::format( "Time to calc hand_fwd_diff: {0:f}ms", timer.elapsedTimeInMilliSec() ) << std::endl; -} \ No newline at end of file