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..244ba132 --- /dev/null +++ b/src/examples/neohookean_timing.cpp @@ -0,0 +1,103 @@ +#include +#include "../tests/neohookean_common.hpp" +#include +#include +#include "axom/core.hpp" +#include "axom/slic/interface/slic.hpp" +#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 }; + double mu = 1.0; + double lambda = 1.0; + int N = 0; + std::cout << "Enter cycles: "; + std::cin >> 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/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..cdc6f6ae --- /dev/null +++ b/src/tests/enzyme_neohookean.cpp @@ -0,0 +1,122 @@ +#include +#include +#include "neohookean_common.hpp" +#include +#include + + +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; + + 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 ); + } + } +}; + + +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_test[9] = { 0.0 }; + stress_FD( E_all + i * 9, mu, lambda, S_FD ); + stress_test_fn( E_all + i * 9, mu, lambda, S_test ); + for ( int j = 0; j < 9; ++j ) { + EXPECT_NEAR( S_FD[j], S_test[j], tol ); + } + } +} + + +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 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 ) +{ + 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 ) +{ + 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 ) + { + 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 ) + { + 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 ) + { + 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 ) +{ + for ( int i = 0; i < 3; ++i ) { + 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 ); + } + } +} diff --git a/src/tests/neohookean_common.hpp b/src/tests/neohookean_common.hpp new file mode 100644 index 00000000..959f379d --- /dev/null +++ b/src/tests/neohookean_common.hpp @@ -0,0 +1,297 @@ +#include +#include +#include +#include +#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] ) +{ + 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 ); + } + } +} + +